Help me to Solve ODE problem
1 回表示 (過去 30 日間)
古いコメントを表示
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 件のコメント
採用された回答
Cris LaPierre
2020 年 9 月 26 日
編集済み: 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.
0 件のコメント
その他の回答 (0 件)
参考
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!