Help me to Solve ODE problem
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
Kevin Isakayoga
am 25 Sep. 2020
Bearbeitet: Cris LaPierre
am 26 Sep. 2020
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 Kommentare
Akzeptierte Antwort
Cris LaPierre
am 26 Sep. 2020
Bearbeitet: Cris LaPierre
am 26 Sep. 2020
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 Kommentare
Weitere Antworten (0)
Siehe auch
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!