ODE System doesn't work when step size (h) is bigger than 1

Question

I am a beginner to Python. Currently I'm writing a code for developing a simple solver for non-linear ODE systems with initial value. The equations of the system are

The function of myu is evaluated first to get the value of myu, then used in dX/dt, dS/dt, and dDO/dt. At the next step, myu is evaluated again to get its new value based on new value of S and DO.

I am using General Linear Method (GLM), proposed by J. C. Butcher, as my method. This method use a transition matrix, which value and size depends on numerical method that we use. In this case, I use Runge Kutta Cash-Karp.

While you may find in the equation that D is also a function, here I set the value of D as a constant.

In initialization, the value of h is set first, to get the number of step. I create a vector named 'initValue', with 8 columns and 4 rows, consist of values of k for each equations (row 1 to 6), initial value for fourth order of the Runge Kutta (row 7. I set it to 0 since it just act as a 'predictor'), and initial value for each equations (row 8).

Transition matrix is created based on the GLM, which values inside it comes from the constants of stage equations (to find the value of k1 to k6) and step equations (to find the solutions) of Runge Kutta Cash-Karp.

In the looping, at the very first time, I simply add the initial values to an array named 'result'. At the first step, I simply multiple the transition matrix with vector 'initValue'. And at the next until final step, I initialize the 'initValue' based on result from previous step.

What I'm looking for is the solution which should look like this. My code works if h is less than 1. I compare my result with result from scipy.integrate.odeint. But when I set h bigger than 1, it show different result than the result it should be. For example, in the code, I set h = 100, which means that it will only display the initial value and final value (when time = 100). While X and S should going upward, and DO and Xr going downward, mine was the opposite of them. The result from odeint when h is set to bigger than 1 show the same result with the expected solution.

I need help to fix my code so it can display the expected solution for any value of h.

Answer 1

Increasing the timestep isn't always sensible either, because you are increasing the numerical error. Instead, you should be setting a requested value for the error, and letting the integrator make the stepsizes which are required to achieve that error estimate (adaptivity), or you need to balance between numerical error and runtime yourself. That said...

I need help to fix my code so it can display the expected solution for any value of h.

Sorry, that's just not possible. This is due to numerical stability issues. Explicit integrators (like explicit RK) have a maximum stepsize that is "stable", above which small perturbations are amplified causing the solver to diverge to infinity. This maximum stepsize is a property of the solver and the "stiffness" of the equation (related to the maximum and minimum eigenvalues of the equation).

If you don't want to be limited by dt at all, you need to use an integrator designed to handle stiffness. "Unlimited dt" only occurs in A-B-L-stable methods, which is an unreasonable goal. Things like the 2nd order Rosenbrock method, Midpoint method, or BDF2 display this though. You can't get higher than 2nd order with that property though, meaning error tends to be higher. So in order to do this well, you need to be adapting both the order and the timestep. This is one reason why I would recommend sticking to using some dedicated ODE solver software: it has a ton of bells and whistles to make this work.