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

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

%dT_p/dt
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
    t(i+1)=t(i)+h;
    %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));
    k1h_p=fh_p(t(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);
    
end

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

Abdu

Answer 1

[Alan Stevens]

Have you tried replacing

u_p(i+1)=u_p(i)+h/6*(k1u_p+ 2*k2u_p + 2*k3u_p +k4u_p);

with

u_p(i+1)= max(u_p(i)+h/6*(k1u_p+ 2*k2u_p + 2*k3u_p +k4u_p),0);

[abdu bouridane]

Hi,

Indeed I have, it just results in the velocity calculation tending close to zero as normal, but then rising again which is what I want to also avoid. I need the iteration for the velocity term to just stop at a value close enough to zero, if that is even possible.

[Alan Stevens]

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);
else
    u_p(i+1) = 0;
    flag = 0;
end

where you calculate u_p(i+1).