Help me to Solve ODE problem
1 visualización (últimos 30 días)
Mostrar comentarios más antiguos
Kevin Isakayoga
el 25 de Sept. de 2020
Editada: Cris LaPierre
el 26 de Sept. de 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 comentarios
Respuesta aceptada
Cris LaPierre
el 26 de Sept. de 2020
Editada: Cris LaPierre
el 26 de Sept. de 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 comentarios
Más respuestas (0)
Ver también
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!