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

Answer 1

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