Help me to Solve ODE problem

Question

Kevin Isakayoga 2020-9-25

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/599674-help-me-to-solve-ode-problem

编辑： Cris LaPierre 2020-9-26

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Cris LaPierre 2020-9-26

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/599674-help-me-to-solve-ode-problem#answer_500878

编辑：Cris LaPierre 2020-9-26

Check your equations. Sr is always zero. This causes cst to always be zero, which then causes drdt to also be zero, causing dRdt to be NaN (drdt/Sr = 0/0 = NaN).

Because the output of the first loop is the initial conditions for the next, and you are getting NaN in the first loop, your problem perpetuates. I tried tweaking the initial conditions some, but it didn't help much.

Check why Sr is always 0.