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 as follow.
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.
Thank you.
Why do you expect any type of reasonable result for ridiculously large step sizes?
The most simple demonstration is y'=-y and the explicit Euler method. If you use step sizes smaller 1 you will get qualititively correct solutions. For step sizes smaller 0.1, you will start to get also quantitatively correct step sizes.
However, if you use a step size h=10, then the iteration
y[k+1]= y[k] - h*y[k] = -9*y[k]
will explode. The same also happens for higher order methods, sufficiently small step sizes give quantitatively correct results, medium step sizes can still give a qualitatively correct picture, large step sizes lead to errors that are very quickly larger than the solution values.