Help me to solve ordinary differential equation

Hi,

I am currently working on some codes and I have to find numerical solutions to the following ODE's using Matlab (with graphs):

clc; clear;
mc=0.3;                                                 %g/cm3
wc=0.4;
pc=3.16;                                                %density of cement g/cm3
pw=1;                                                   %density of water g/cm3
T=293;
vc=1/((wc*pc)+1);                                       %the cement volume fraction (dalam bentuk 1 satuan)
ro=(3*vc/(4*pi))^1/3;                                   %jari2 partikel awal (for ex: 2, 4, 8, 16, 32) (dr micron jadi ke mm) tapi masi ratio
L=((4*pi*(wc*pc+1)/3)^(1/3))*ro;                        %nilainya 1   
b1=1231;                                                %K^-1
b2=7579;                                                %K^-1
C=5*10^7;                                               %/mm2h
ER=5364;                                                %K^-1
yg=0.25;                                                %the stoichiometric ratio by mass of water to cement
yw=0.15;                                                %the ratio of water adsorbed in gel pore to reacted cement
RH=0.88;                                                %asumsi
cw=((RH-0.75)/0.25)^3;
v=2;%, ratio of volume of cement gel including gel pore assumed to be 2
rwc3s=0.586;                                            %mass contents of C3S
rwc3a=0.172;                                            %mass contents of C3A
De293=((rwc3s*100)^2.024)*(3.2*10^-14);                 %mm2/h
kr293=(8.05*(10^-10))*((rwc3s*100+rwc3a*100)^0.975);    %mm/h
B293=0.3*10^-10;                                        %mm/h
kr=kr293*exp(-ER*(1/T-1/293));%mm/h
%De=De293*(log(1./alfa)).^1.5;
B=B293*exp(-b1*(1/T-1/293));%mm/h
%kd=(B/alfa^1.5)+C*(Rt-ro).^4; %(mm/h)
tspan = [0 100];   
Rt=0.7;% in the range between 0.5-0.85
syms alfa(t)
input1= diff(alfa,t)==((3*cw)/((yw+yg)*pc*ro^2))*1/((1/(((B/alfa^1.5)+C*(Rt-ro)^4)*ro*alfa^(2/3)))+(((alfa^(-1/3))-((2-alfa)^(-1/3)))/(De293*(log(1/alfa))^1.5))+(1/(kr*ro*alfa^(-2/3))));
cond = alfa(0)==0;
Solve_alfa(t)=ode45(input1,tspan,cond); %or should i use dsolve?

and i also tried to make another code the same but still got NaN value.

clc; clear;
tspan=[0 1000];
y0=[0;0];
 
[t,r]=ode45(@myod,tspan,y0);
plot(t,r);
lgd=legend('r(t)','R(t)');
lgd.FontSize=10;
 
function drdt=myod(t,r) 
wc=0.4;
pc=3.16;                                                %density of cement g/cm3
T=293;
vc=1/((wc*pc)+1);                                       %the cement volume fraction (dalam bentuk 1 satuan)
ro=(3*vc/(4*pi))^1/3;                                   %jari2 partikel awal (for ex: 2, 4, 8, 16, 32) (dr micron jadi ke mm) tapi masi ratio
b1=1231;                                                %K^-1
C=5*10^7;                                               %/mm2h
ER=5364;                                                %K^-1
yg=0.25;                                                %the stoichiometric ratio by mass of water to cement
yw=0.15;                                                %the ratio of water adsorbed in gel pore to reacted cement
RH=0.88;                                                %asumsi
cw=((RH-0.75)/0.25)^3;
v=2;%, ratio of volume of cement gel including gel pore assumed to be 2
rwc3s=0.586;                                            %mass contents of C3S
rwc3a=0.172;                                            %mass contents of C3A
De293=((rwc3s*100)^2.024)*(3.2*10^-14);                 %mm2/h
kr293=(8.05*(10^-10))*((rwc3s*100+rwc3a*100)^0.975);    %mm/h
B293=0.3*10^-10;                                        %mm/h
kr=kr293*exp(-ER*(1/T-1/293));%mm/h
%De=De293*(log(1./alfa)).^1.5;
B=B293*exp(-b1*(1/T-1/293));%mm/h
%kd=(B/alfa^1.5)+C*(Rt-ro).^4; %(mm/h)
Rt=0.75;
    if Rt<=ro
    Sr=0;
    elseif (Rt>=ro)&&(Rt<0.5)
    Sr=4*pi*(Rt)^2;
    elseif (0.5<=Rt)&&(Rt<0.5*(2^0.5))
    Sr=(4*pi*(Rt)^2)-(6*pi*(1-(0.5/(Rt))));
    elseif (Rt>=0.5*(2^0.5)) && (Rt<0.5*(3^0.5))
    syms y x 
    fun = @(y,x) 8*(Rt)./(sqrt((Rt)^2-(x.^2)-(y.^2)));
    ymin=sqrt((Rt)^2-0.5);
    xmin=@(x) sqrt((Rt)^2-0.25-x.^2);
    Sr=integral2(fun,ymin,0.5,xmin,0.5);
    else
    Sr=0;
    end
pw=1;
drdt=zeros(2,1); 
drdt(1)=-((pw*(Sr/(4*pi*(r(2)^2)))*cw)/((yw+yg)*pc*r(1)^2))*1/((1/(((B/(1-(r(1)/ro)^3)^1.5)+C*(r(2)-ro)^4)*r(1)^2))+(((1/r(1))-(1/r(2)))/(De293*(log(1/((1-(r(1)/ro)^3))))^1.5))+(1/(kr*r(1)^2))); 
drdt(2)=(v-1)*(r(1)^2)*(-((pw*(Sr/(4*pi*(r(2)^2)))*cw)/((yw+yg)*pc*r(1)^2))*1/((1/(((B/(1-(r(1)/ro)^3)^1.5)+C*(r(2)-ro)^4)*r(1)^2))+(((1/r(1))-(1/r(2)))/(De293*(log(1/((1-(r(1)/ro)^3))))^1.5))+(1/(kr*r(1)^2))))/Sr;
end

I am new to Matlab and can't seem to figure out how to do these 2 codes. Would it be possible for someone to give me the Matlab codes/instructions I need to find the numerical solutions to these ODE's?

Thank you very much !

[Alan Stevens]

Are you sure you have your initial conditions right?

In your first code you have a term

log(1/alfa))^1.5

in input1. With alfa(0) = 0, this will cause problems. Similarly, in your second code you have

(1/(kr*r(1)^2)))

in drdt(1). With r(1) = 0, this will also cause problems.

[madhan ravi]

Are you sure you want everyone to see your thesis calculations a public forum? Because once you get the complete answer , you will want to delete your question , which makes the efforts of the answerer go to waste. P.S: This thread has been backed up.

[madhan ravi]

Ok, good :)

Answer 1

[Alan Stevens]

The following gets it working, though the output isn't very exciting!. However, you have a number of variables that aren't used.

tspan = [0 100];   

alfa0=0.0001;
[t,alfa]=ode45(@f,tspan,alfa0); %or should i use dsolve?

plot(t,alfa),grid

function dalfadt = f(~,alfa)
%mc=0.3;                         %%%% unused            %g/cm3
wc=0.4;
pc=3.16;                                                %density of cement g/cm3
%pw=1;                           %%%% unused            %density of water g/cm3
T=293;
vc=1/((wc*pc)+1);                                       %the cement volume fraction (dalam bentuk 1 satuan)
ro=(3*vc/(4*pi))^1/3;                                   %jari2 partikel awal (for ex: 2, 4, 8, 16, 32) (dr micron jadi ke mm) tapi masi ratio
%L=((4*pi*(wc*pc+1)/3)^(1/3))*ro;    %%%% unused        %nilainya 1   
b1=1231;                                                %K^-1
%b2=7579;                            %%%% unused        %K^-1
C=5*10^7;                                               %/mm2h
ER=5364;                                                %K^-1
yg=0.25;                                                %the stoichiometric ratio by mass of water to cement
yw=0.15;                                                %the ratio of water adsorbed in gel pore to reacted cement
RH=0.88;                                                %asumsi
cw=((RH-0.75)/0.25)^3;
%v=2;%, %%%% unused  ratio of volume of cement gel including gel pore assumed to be 2
rwc3s=0.586;                                            %mass contents of C3S
rwc3a=0.172;                                            %mass contents of C3A
De293=((rwc3s*100)^2.024)*(3.2*10^-14);                 %mm2/h
kr293=(8.05*(10^-10))*((rwc3s*100+rwc3a*100)^0.975);    %mm/h
B293=0.3*10^-10;                                        %mm/h
kr=kr293*exp(-ER*(1/T-1/293));%mm/h
%De=De293*(log(1./alfa)).^1.5;
B=B293*exp(-b1*(1/T-1/293));%mm/h
%kd=(B/alfa^1.5)+C*(Rt-ro).^4; %(mm/h)
Rt=0.7;% in the range between 0.5-0.85

dalfadt=((3*cw)/((yw+yg)*pc*ro^2))*1/((1/(((B/alfa^1.5)+C*(Rt-ro)^4)*ro*alfa^(2/3)))+(((alfa^(-1/3))-((2-alfa)^(-1/3)))/(De293*(log(1/alfa))^1.5))+(1/(kr*ro*alfa^(-2/3))));
end

[Jamba wumbo]

Thank you very much Alan ! I will fix the rest of the code.

And I also have one last question, how to plot the parameter of Rt with alfa in the same plot? for example Rt 0.7 and 0.75?

Rt=0.7;% in the range between 0.5-0.85

and the time is always follows the unit of equation itself right? if my equation is in hour, so the time is also in hour?

Thankyou in advance.

[Alan Stevens]

Wel you could use the following

tspan = [0 100];   
Rt = [0.5 0.7 0.85];
for i = 1:numel(Rt)
alfa0=0.0001;
[t,alfa(:,i)]=ode45(@f,tspan,alfa0,[],Rt(i)); %or should i use dsolve?
end

plot(t,alfa),grid

function dalfadt = f(~,alfa,Rt)
%mc=0.3;                         %%%% unused            %g/cm3
wc=0.4;
pc=3.16;                                                %density of cement g/cm3
%pw=1;                           %%%% unused            %density of water g/cm3
T=293;
vc=1/((wc*pc)+1);                                       %the cement volume fraction (dalam bentuk 1 satuan)
ro=(3*vc/(4*pi))^1/3;                                   %jari2 partikel awal (for ex: 2, 4, 8, 16, 32) (dr micron jadi ke mm) tapi masi ratio
%L=((4*pi*(wc*pc+1)/3)^(1/3))*ro;    %%%% unused        %nilainya 1   
b1=1231;                                                %K^-1
%b2=7579;                            %%%% unused        %K^-1
C=5*10^7;                                               %/mm2h
ER=5364;                                                %K^-1
yg=0.25;                                                %the stoichiometric ratio by mass of water to cement
yw=0.15;                                                %the ratio of water adsorbed in gel pore to reacted cement
RH=0.88;                                                %asumsi
cw=((RH-0.75)/0.25)^3;
%v=2;%, %%%% unused  ratio of volume of cement gel including gel pore assumed to be 2
rwc3s=0.586;                                            %mass contents of C3S
rwc3a=0.172;                                            %mass contents of C3A
De293=((rwc3s*100)^2.024)*(3.2*10^-14);                 %mm2/h
kr293=(8.05*(10^-10))*((rwc3s*100+rwc3a*100)^0.975);    %mm/h
B293=0.3*10^-10;                                        %mm/h
kr=kr293*exp(-ER*(1/T-1/293));%mm/h
%De=De293*(log(1./alfa)).^1.5;
B=B293*exp(-b1*(1/T-1/293));%mm/h
%kd=(B/alfa^1.5)+C*(Rt-ro).^4; %(mm/h)
%Rt=0.7;% in the range between 0.5-0.85

dalfadt=((3*cw)/((yw+yg)*pc*ro^2))*1/((1/(((B/alfa^1.5)+C*(Rt-ro)^4)*ro*alfa^(2/3)))+(((alfa^(-1/3))-((2-alfa)^(-1/3)))/(De293*(log(1/alfa))^1.5))+(1/(kr*ro*alfa^(-2/3))));
end

but the results are clearly insensitive to Rt in this range!

[Jamba wumbo]

I think this is the modified one by me, however I still stuck in the matrix size.

tspan =[0 100];   
ro = [0.004 0.008 0.016];
for i = 1:numel(ro)
alfa0=0.0001;
[t,alfa(:,i)]=ode45(@f,tspan,alfa0,[],ro(i)); %or should i use dsolve?
end
plot(t,alfa),grid
function dalfadt = f(~,alfa,ro)
pc=3.16*0.001;                                          %density of cement g/mm3
T=297;
%vc=1/((wc*pc)+1);                                      %the cement volume fraction (dalam bentuk 1 satuan)
%ro=2*0.001;                                             %jari2 partikel awal (for ex: 2, 4, 8, 16, 32) (dr micron jadi ke mm) tapi masi ratio
b1=1231;                                                %K^-1
C=5*10^7;                                               %/mm2h
ER=5364;                                                %K^-1
yg=0.25;                                                %the stoichiometric ratio by mass of water to cement
yw=0.15;                                                %the ratio of water adsorbed in gel pore to reacted cement
RH=0.88;                                                %asumsi
cw=((RH-0.75)/0.25)^3;
rwc3s=0.586;                                            %mass contents of C3S
rwc3a=0.172;                                            %mass contents of C3A
De293=((rwc3s*100)^2.024)*(3.2*10^-14);                 %mm2/h
kr293=(8.05*(10^-10))*((rwc3s*100+rwc3a*100)^0.975);    %mm/h
B293=0.3*10^-10;                                        %mm/h
kr=kr293*exp(-ER*(1/T-1/293));                          %mm/h
De=De293*(log(1/alfa))^1.5;
B=B293*exp(-b1*(1/T-1/293));                            %mm/h
kd=(B/alfa^1.5)+C*(alfa)^2;                             %mm/h
dalfadt=((3*cw)/((yw+yg)*pc*ro^2))*1/((1/((kd)*ro*alfa^(2/3)))+(((alfa^(-1/3))-((2-alfa)^(-1/3)))/(De))+(1/(kr*ro*alfa^(-2/3))));
end

Could you help me where should I change?

Thank you in advance, sorry to bother you a lot.

[Alan Stevens]

The problem arises from the fact that there were happened to be a different number of timesteps for each value of ro. Overcome this by forcing the number of timesteps to be the same each time. Begin the program with:

tspan = 0:100;   
ro = [0.004 0.008 0.016];
for i = 1:numel(ro)
alfa0=0.0001;
[t,alfa(:,i)]=ode45(@f,tspan,alfa0,[],ro(i)); 
lgnd(i,:) = num2str(ro(i));
end
plot(t,alfa),grid
legend(lgnd)

[Jamba wumbo]

And about ODE timestep, is it based on our equation unit? for example if the timesteps at 100 so the time will correspond to 100 hours, Is it like that?

Thank you! You helped me a lot. I truly thankful to you.

Best Regard,

Hadiyoga

[Alan Stevens]

"And about ODE timestep, is it based on our equation unit? for example if the timesteps at 100 so the time will correspond to 100 hours, Is it like that?"

Yes, if the variables that feed in to dalfadt use hours, then that is the case.