Hi every one, I am trying to solve the ode function, but up until now I still couldnot figure it out why my value always NaN.
Here below is my function,
function f=model2(~,Y,ro) rt=Y(1);Rt=Y(2);% Explicit equations
pw=1*0.001; pc=3.16*0.001; wc=pw/pc*((1/((0.5^3)*4/3*pi))-1); T=293; b1=1231; C=5*10^7; ER=5364; yg=0.25; yw=0.15; RH=0.8; cw=((RH-0.75)/0.25)^3;v=2; rwc3s=0.586; rwc3a=0.172; De293=((rwc3s*100)^2.024)*(3.2*10^-14); kr293=(8.05*(10^-10))*((rwc3s*100+rwc3a*100)^0.975); B293=0.3*10^-10; alfa=1-(rt/ro)^3;kr=kr293*exp(-ER*(1/T-1/293)); De=De293*(log(1/alfa))^1.5;B=B293*exp(-b1*(1/T-1/293)); kd=(B/alfa^1.5)+C*(Rt-ro)^4; L=((4*pi*(wc*(pc/pw)+1)/3)^(1/3))*ro*0.001; Rt1=Rt/L;z=ro/L;if Rt1<=z Sr=0;elseif (Rt1>=z)&&(Rt1<0.5) Sr=4*pi*(Rt1)^2;elseif (0.5<=Rt1)&&(Rt1<0.5*(2^0.5)) Sr=(4*pi*(Rt1)^2)-(12*pi*(1-(0.5/(Rt1))));elseif (Rt1>=0.5*(2^0.5)) && (Rt1<0.5*(3^0.5)) syms y x fun = @(y,x) 8*(Rt1)/(sqrt((Rt1)^2-(x.^2)-(y.^2))); ymin=sqrt((Rt1)^2-0.5); xmin=@(x) sqrt((Rt1)^2-0.25-x.^2); Sr=integral2(fun,ymin,0.5,xmin,0.5);else Sr=0;end cst=Sr/(4*pi*(Rt^2));%Differential equations
drtdt=-((pw*cst*cw)/((yw+yg)*pc*rt^2))*1/((1/(kd*rt^2))+(((1/rt)-(1/Rt))/De)+(1/(kr*rt^2))); dRtdt=(v-1)*(rt^2)*drtdt/Sr;f=[drtdt;dRtdt];end
And here is the command to run the function
clc; clear;tspan=[0 100];ro=4*0.001; y0=[(ro-0.0001) (ro+0.0001)];[t,y]=ode45(@model2,tspan,y0,[],ro);figure (2)plot(t,y(:,1),t,y(:,2));legend('rt','Rt');ylabel('mm');xlabel('t');axis([tspan 0 100]),grid;
Thank you for you attention! Looking forward to help my question.
Best Regard,
Kevin
Best Answer