ODE system initial condition results in infinity/NaN error. Initial conditions cannot be changed. How to evade?

Hi all. I currently have my ODE system code 90% complete using the runge-kutta 4th order(R4K) iteration method but I've hit a bit of a snag due to both my initial conditions, and the diferential invloved.

For one of my differentials I'm calculating the change in velocity, where u_p(0) = 0

fu_p=@(t,D_p,Tp,u_p) -3.*CD(t)*rhoG*(u_p^2)/(4.*rhoP*D_p)+((rhoP-rhoG)/rhoP)*g;

Now the issue is, u_p is related to 2 other non ODE's in the system:

Re =@(t) (rhoG*D_p*u_p)/muG;
CD =@(t) 24./Re(t);

Clearly we can see that when u_p(0) =0, Re(0) =0 which then leads to CD(0)= inf

This completely breaks the differential as It results in NaN values at every time step, which also breaks the other differentials in the system(although for now they aren't related to the root cause).

I have tried incorporating a loop where it would use a different function for u_p at the initial time t=0, and use another when t>0 but it didn't seem to work as the function is above the R4K loop section.

Is there any way of "bypassing" this in a quick line, or am I just missing something simple.

Answer 1

[Alan Stevens]

Simply write your differential function as:

fu_p = @(D_p,u_p) -18*muG*u_p/(rhoP*D_p^2) + g*(rhoP-rhoG)/rhoP;

avoiding any explicit reference to Re or CD.

(Incidentally, although you had Tp as an input argument to fu_p, it isn't clear where that is used, so I've left it out of the expression above!

If, however, rho_p is a function ofTp, then you should probably replace rhoP in the function by rhoP(Tp).

[abdu bouridane]

Thanks! Much appreciated indeed! I simplified that CD equation for ease of troubleshooting but the method would be exactly the same, just rearrange for Re, and sub it into the function.

I do have other issues that I've noticed, but that comes later!