Solving Non-Linear ODEs with RK4 Method: Time Step Considerations

[RandomGuy88]

Hi,

I am having a problem implementing 4th Order Runge Kutta (RK4) to solve a non linear ODE.

dy/dt = B(t) - A(t)*y(t)^2 where the coefficients A(t) and B(t) are specific to the problem I am working on.

For the test case I am running I specify a sinusoidal signal and from this signal I calculate A(t) and B(t). While the magnitudes of these coefficients are functions of time they are not functions of the frequency of the sinusoidal function I am specify. I am sure you are all familiar with the algorithm but having it will help me explain.

y(i+1) = y(i) + h/6*(k1 + 2k2 + 2k3 + k4)

k1, k2, k3 and k4 are functions of A(t), B(t) and y(t). As I mentioned above the coefficients A and B are not functions of the frequency of the sinusoid so the values of k1, k2, k3 and k4 are also not functions of the frequency.

My problem is that the values of the k coefficients are quite large. So when I use a high frequency signal the large values of k are compensated by the very small value of h so that the second term on the right hand side h/6*(k1 ...) is the same order of magnitude as y(i). If the frequency is small then I would think I could use a larger timestep h but the k coefficients are still very large so h/6*(k1 ...) >> y(i) and the method is unstable.

Is there anything I can do other than just using a very small time step for the low frequency cases? Perhaps some way of scaling the coefficients? Or should I just use a different method? It just seems wrong that for a lower frequency I need a smaller time step.

Thanks for your help.

 

[the_wolfman]

Consider what happens if you take make A and B constants in time. This is the ultimate low frequency approximation.

For instance let both A and B =1

Then your equations has the solution \frac{\left|y+1 \right|}{\left|y-1\right|}=Ae^{2t}

Notice that the solution is divided into three distinct regions y>1, -1<y<1 and y<-1. Your solution should not cross from one region into another. Solutions in these three regions are distinct.


However if your time step is too large, than you RK4 method will calculating k2,k3,k4 from multiple regions. This will give you an erroneous result. I don't know if this is the root cause of your numerical instability... but it does place limits on how large your time step can be.

 

[AlephZero]

The "size" of the $k$ coefficients is not a stability issue in itself. If you solve the same problem in different units, the numerical value of coefficients will change but the number of time steps required for a stable solution will not.

I suggest you take constant values of $A$ and $B$ as in post #2, and make a phase-plane diagram (a plot of $y'$ against $y$) for different values of $A$ and/or $B$. You can plot that direct from your differential equation, and also try to reproduce the curves with your numerical solution.

The analysis in #2 might be the explanation of what is going on. The situation may be worse than that, because if $A$ and $B$ are time-dependent, the boundaries between the regions will also be functions of time, which may force the solution of the DE to go to infinity if it can't be a continuous function that crosses the boundary.

Assuming your DE represents the physics correctly, the real problem might be choosing the correct initial conditions for your solutions. You might consider integrating the solution backwards in time, if it's easier to choose physically sensible conditions at the final value of $t$.