Solving system of three differential equations using Runge-Kutta 4 in python

Question

I wrote code for Runge-Kutta 4 for solving system of 3 ODEs

I think that it does not work fine for because I solved the system with Euler's method and I had have the following results

But with the RK4's code attached I get a very different result:

import matplotlib.pyplot as plt    
umax = 0.53
Km = 0.12
Yxs = 0.4
S1f = 4
Yxp = 0.4
D=0.5      #
X1f=0.2
P1f=0.1

# Equations:


#du/dt=V(u,t)
def V(u,t):
    X1, S1, P1, vx, vs, vp = u
    return np.array([ vx,vs,vp,
            D*(X1f - X1)+(umax*(1-np.exp(-S1/Km)))*X1, 
            D*(S1f - S1)-(umax*(1-np.exp(-S1/Km))*X1)/Yxs, 
            D*(P1f - P1)+(umax*(1-np.exp(-S1/Km)))*X1*Yxp ])

def rk4(f, u0, t0, tf , n):
    t = np.linspace(t0, tf, n+1)
    u = np.array((n+1)*[u0])
    h = (t[1]-t[0])/n
    for i in range(n):
        k1 = h * f(u[i], t[i])    
        k2 = h * f(u[i] + 0.5 * k1, t[i] + 0.5*h)
        k3 = h * f(u[i] + 0.5 * k2, t[i] + 0.5*h)
        k4 = h * f(u[i] + k3, t[i] + h)
        u[i+1] = u[i] + (k1 + 2*(k2 + k3 ) + k4) / 6
    return u, t

u, t  = rk4(V, np.array([0., 0., 0. , 0., 1. , 0.]) , 0. , 40. , 4000)
x,y,z, vx,vs,vp  = u.T
# plt.plot(t, x, t,y)
plt.plot(t, x, t,y,t,z)
plt.grid('on')

plt.show()

I have reviewed the code several times but I can't find the reason why the result is very different from the previous one

It would be helpful for reproducing the difference if you could add the Euler code. (In my opinion I would keep the separation of text/image problem description with a code appendix, perhaps with the Euler code last.) If it is not the case already, you could adapt the Euler code to be similar in structure to the RK4 code, using the same `V` function. Then again, this might already solve the problem of different outcomes, so in the end the original Euler code is needed to find the critical differences. — Lutz Lehmann, Feb 05 '23 at 11:53

score 2 · Accepted Answer · answered Feb 05 '23 at 16:00

You should notice that the function V there will only be the variables X1, S1, P1 = u and the return vector as

np.array([D*(X1f - X1)+(umax*(1-np.exp(-S1/Km)))*X1, D*(S1f - S1)-(umax*(1-np.exp(-S1/Km))*X1)/Yxs, D*(P1f - P1)+(umax*(1-np.exp(-S1/Km)))*X1*Yxp ])

So your code would look like the following

import numpy as np
import matplotlib.pyplot as plt


umax = 0.53
Km = 0.12
Yxs = 0.4
S1f = 4
Yxp = 0.4
D=0.5      # Para teissier 0<D<0.514
X1f=0.2
P1f=0.1

# Equations:
#du/dt=V(u,t)
def V(u,t):
    X1, S1, P1 = u
    return np.array([D*(X1f - X1)+(umax*(1-np.exp(-S1/Km)))*X1,D*(S1f -  S1)-(umax*(1-np.exp(-S1/Km))*X1)/Yxs,D*(P1f - P1)+(umax*(1-np.exp(-S1/Km)))*X1*Yxp] )

def rk4(f, u0, t0, tf , n):
    t = np.linspace(t0, tf, n+1)
    u = np.array((n+1)*[u0])
    h = (t[1]-t[0])
    for i in range(n):
        k1 = h * f(u[i], t[i])    
        k2 = h * f(u[i] + 0.5 * k1, t[i] + 0.5*h)
        k3 = h * f(u[i] + 0.5 * k2, t[i] + 0.5*h)
        k4 = h * f(u[i] + k3, t[i] + h)
        u[i+1] = u[i] + (k1 + 2*(k2 + k3 ) + k4) / 6
    return u, t

u, t  = rk4(V, np.array([0., 0., 1. ]) , 0. , 40. , 4000)
x,y,z= u.T
fig = plt.figure()                                  # create figure
plt.plot(t, x, linewidth = 2, label = 'Biomasa')    # plot Y to t 
plt.plot(t, y, linewidth = 2, label = 'Sustrato')    # plot P to tplt.title('Title', fontsize = 12)    # add some title to your plot
plt.plot(t, z, linewidth = 2, label = 'Producto')    # plot P to tplt.title('Title', fontsize = 12)    # add some title to your plot
plt.xlabel('t (en segundos)', fontsize = 12)
plt.ylabel('Biomasa, Sustrato, Producto', fontsize = 12)
plt.xticks(fontsize = 12)
plt.yticks(fontsize = 12)
plt.grid(True)                        # show grid
plt.legend()
plt.show()

Thus obtaining the graphics you want

score 1 · Answer 2 · answered Feb 04 '23 at 07:55

1

You did too much in

    h = (t[1]-t[0])/n

The difference t[1]-t[0] is already the step size, what you did makes it the step size for n^2 steps in the same interval.

Are the initial values really the same? In the first graph it does not look like all variables start at zero.

answered Feb 04 '23 at 07:55

Lutz Lehmann

25,219
2
22
51

Yes, the initial values are the sames but removing the n from h =(t[1]-t[0])/n, I continue with the same problem – Juan Saknussem Feb 05 '23 at 05:32
2

I suspect something is wrong in your definition of `V`. I implemented a quick Euler integration and I get similar results as your RK4 approach after fixing the problem Lutz pointed out: https://i.stack.imgur.com/1Iqac.png – Pranav Hosangadi Feb 05 '23 at 06:42
I had to significantly increase the resolution of the Euler integration, but that was expected. Your RK4 integration looks alright to me (other than what this answer identifies) – Pranav Hosangadi Feb 05 '23 at 06:43
Yep, your RK4 is fine. Here's a solution for V = cos(t). https://i.stack.imgur.com/ZZGjc.png – Pranav Hosangadi Feb 05 '23 at 06:58
Indeed, I found that the problem lies in the definition of the function V and I deleted the variables vx, vs and vp and in turn the vector of initial values is of 3 components and not 6 as was proposed – Juan Saknussem Feb 05 '23 at 15:51

Solving system of three differential equations using Runge-Kutta 4 in python

2 Answers2