MATLAB Answers

How to stop ODE(R4K) function when the value is close to 0?

49 views (last 30 days)
abdu bouridane
abdu bouridane on 26 Feb 2021
Commented: Alan Stevens on 28 Feb 2021
Follow up from my previous post about my ODE system using the runge-kutta 4th order method.
The main body of the code looks like this:
(it could be far cleaner I'm sure, but I am a novice at MATLAB so I didn't want to change much in case It results in errors)
fD_p=@(t,D_p,Tp,u_p) (-2.*kc(t,D_p,u_p).*(C_s(t,Tp)-C_inf))./rhoP;%differential eqn for diameter change of droplet
fTp=@(t,D_p,Tp,u_p) (6.*h_m(t,D_p,u_p)*(Ta-Tp)-(6.*L(t,Tp).*kc(t,D_p,u_p).*(C_s(t,Tp)-C_inf)))./(rhoP*cpP.*D_p);% differential eqn for temperature change of droplet
%du_p/dt, velocity differential "u_p"
fu_p=@(t,D_p,Tp,u_p) -18*muG.*u_p./(rhoP.*D_p.^2)+g.*(1-(rhoG./rhoP)); % differential eqn for velocity of droplet
fh_p=@(t,u_p) u_p; % differential for vertical height travelled as a function of time and velocity
%step size
h=0.0001; % step size
tfinal=18.4; % final time
N=tfinal/h; % number of iterations
%update loop
for i=1:N
%update time
%update D_p,T_p and u_p
%Runge-Kutta 4th order code
k1D_p= fD_p(t(i), D_p(i), Tp(i),u_p(i));
k1Tp= fTp(t(i), D_p(i),Tp(i),u_p(i));
k1u_p= fu_p(t(i), D_p(i), Tp(i),u_p(i));
k2D_p= fD_p(t(i)+h/2, D_p(i)+h/2*k1D_p, Tp(i)+h/2*k1Tp, u_p(i)+h/2*k1u_p);
k2Tp= fTp( t(i)+h/2, D_p(i)+h/2*k1D_p, Tp(i)+h/2*k1Tp, u_p(i)+h/2*k1u_p);
k2u_p= fu_p(t(i)+h/2, D_p(i)+h/2*k1D_p, Tp(i)+h/2*k1Tp, u_p(i)+h/2*k1u_p);
k2h_p= fh_p(t(i)+h/2, u_p(i)+h/2*k1u_p);
k3D_p= fD_p(t(i)+h/2, D_p(i)+h/2*k2D_p, Tp(i)+h/2*k2Tp, u_p(i)+h/2*k2u_p);
k3Tp= fTp( t(i)+h/2, D_p(i)+h/2*k2D_p, Tp(i)+h/2*k2Tp, u_p(i)+h/2*k2u_p);
k3u_p= fu_p(t(i)+h/2, D_p(i)+h/2*k2D_p, Tp(i)+h/2*k2Tp, u_p(i)+h/2*k2u_p);
k3h_p= fh_p(t(i)+h/2, u_p(i)+h/2*k2u_p);
k4D_p= fD_p(t(i)+h,D_p(i)+h*k3D_p,Tp(i)+h*k3Tp, u_p(i)+h*k3u_p);
k4Tp= fTp(t(i)+h, D_p(i)+h*k3D_p,Tp(i)+h*k3Tp, u_p(i)+h*k3u_p);
k4u_p= fu_p(t(i)+h, D_p(i)+h*k3D_p,Tp(i)+h*k3Tp, u_p(i)+h*k3u_p);
k4h_p= fh_p(t(i)+h, u_p(i)+h*k3u_p);
%update D_p value
D_p(i+1)=D_p(i)+h/6*(k1D_p+ 2*k2D_p + 2*k3D_p +k4D_p);
%update Tp value
Tp(i+1)=Tp(i)+h/6*(k1Tp+ 2*k2Tp + 2*k3Tp +k4Tp);
%update u_p value
u_p(i+1)=u_p(i)+h/6*(k1u_p+ 2*k2u_p + 2*k3u_p +k4u_p);
%update h_p value
h_p(i+1)=h_p(i)+h/6*(k1h_p+ 2*k2h_p + 2*k3h_p +k4h_p);
The other vairables in the functions which call at time t,such as
h_m(t,D_p,u_p) % heat transfer coefficient, depends on the diameter and velocity of particle at time t.
..are calculated outside the R4K loop as they arent differentials, they depend on the parameters at each time step t.
Now on to the issue at hand. If you look at the velocity differential, eventually it starts to decrease and get close to 0 which is expected. The issue is, the iteration keeps going past 0, into the negatives, resulting in the other 3 differentials to break. Is there a way to stop the iteration for it,such that the differential becomes a constant at or as close to 0 as possible?
I cant assign a loop for each differential since they all depend on each other,so I'm a bit confused on what to do..
Many thanks

Answers (1)

Alan Stevens
Alan Stevens on 27 Feb 2021
Have you tried replacing
u_p(i+1)=u_p(i)+h/6*(k1u_p+ 2*k2u_p + 2*k3u_p +k4u_p);
u_p(i+1)= max(u_p(i)+h/6*(k1u_p+ 2*k2u_p + 2*k3u_p +k4u_p),0);
Alan Stevens
Alan Stevens on 28 Feb 2021
Try setting
flag = 1;
before your for loop, then have
if (u_p(i)>0) && (flag == 1)
u_p(i+1)=u_p(i)+h/6*(k1u_p+ 2*k2u_p + 2*k3u_p +k4u_p);
u_p(i+1) = 0;
flag = 0;
where you calculate u_p(i+1).

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!