Applying Forward Euler Method to a Three-Box Model System of ODEs

Question

I am trying to model a system of coupled ODEs which represent a three-box ocean model of phosphorous concentration (y) in the low-latitude surface ocean (Box 1), high-latitude deep ocean (Box 2), and deep ocean (Box 3). The ODEs are given below:

dy1dt = (Q/V1)*(y3-y1) - y1/tau                                        # dP/dt in Box 1
dy2dt = (Q/V2)*(y1-y2) + (qh/V2)*(y3-y2) - (fh/V2)                     # dP/dt in Box 2
dy3dt = (Q/V3)*(y2-y3) + (qh/V3)*(y2-y3) + (V1*y1)/(V3*tau) + (fh/V3)  # dP/dt in Box 3

The constants and box volumes are given by:

### Define Constants
tau = 86400              # s
VT  = 1.37e18            # m3 
Q   = 25e6               # m3/s 
qh  = 38e6               # m3/s
fh  = 0.0022             # mol/m3
avp = 0.00215            # mol/m3
    
### Calculate Surface Area of Ocean
r = 6.4e6                # m
earth = 4*np.pi*(r**2)   # m2
ocean = .70*earth        # m2   
    
### Calculate Volumes for Each Box
V1 = .85*100*ocean       # m3
V2 = .15*250*ocean       # m3
V3 = VT-V1-V2            # m3

This can be put into matrix form y = Ay + f, where y = [y1, y2, y3]. I have provided the matrices and initial conditions below:

A = np.array([[(-Q/V1)-(1/tau),0,(Q/V1)],
                  [(Q/V2),(-Q/V2)-(qh/V2),(qh/V2)],
                  [(V1/(V3*tau)),((Q+qh)/V3),((-Q-qh)/V3)]])
f = np.array([[0],[-fh/V2],[fh/V3]])

y1 = y2 = y3 = 0.00215 # mol/m3

I am having trouble adapting the Forward Euler method to apply to a system on linear ODEs, rather than just one. This is what I have come up with so far (it runs with no issues but doesn't work if that makes sense; I think it has something t do with initial conditions?):

### Define a Function for the Euler-Forward Scheme
import numpy as np
def ForwardEuler(t0,y0,N,dt):
        N = 100
        dt = 0.1
        # Create empty 2D arrays for t and y
        t = np.zeros([N+1,3,3]) # steps, # variables, # solutions
        y = np.zeros([N+1,3,3])
        # Assign each ODE to a vector element
        y1 = y[0]
        y2 = y[1]
        y3 = y[2]
        # Set initial conditions for each solution
        t[0, 0] = t0[0] 
        y[0, 0] = y0[0]
        t[0, 1] = t0[1] 
        y[0, 1] = y0[1]
        t[0, 2] = t0[2] 
        y[0, 2] = y0[2]
        
        for i in trange(int(N)):
            t[i+1]  = t[i]  + t[i]*dt
            y1[i+1] = y1[i] + ((Q/V1)*(y3[i]-y1[i]) - (y1[i]/tau))*dt
            y2[i+1] = y2[i] + ((Q/V2)*(y1[i]-y2[i]) + (qh/V2)*(y3[i]-y2[i]) - (fh/V2))*dt
            y3[i+1] = y3[i] + ((Q/V3)*(y2[i]-y3[i]) + (qh/V3)*(y2[i]-y3[i]) + (V1*y1[i])/(V3*tau) + (fh/V3))*dt
        return t, y1, y2, y3

Any help on this is greatly appreciated. I have not found any resources online that go through the Euler Forward for a system of 3 ODEs, and am at a loss. I am happy to explain further if there are more questions.

Why is `t` a 2D array, or a list of 2D arrays? Why does `y` have the format `3x3`? Why do you define the matrix formulation and then do not use it? You will catch more errors, simplify testing, and avoid bad design decisions by separating out the ODE function and the Euler integration as their own procedures. This would also simplify switching to a higher-order method or using the `scipy.integrate` steppers and solvers. — Lutz Lehmann, Feb 11 '22 at 15:12
@LutzLehmann I was trying to adapt a script I found; I am really new to Python and don't have a great understanding of how Python works through matrix math or ODEs in general. What do you mean by separating out the ODEs and Euler integration? How does that work? — Savannah, Feb 11 '22 at 21:07
The structure would be similar to https://stackoverflow.com/a/66883561/3088138, with the extension to ODE systems as in https://stackoverflow.com/questions/47711186/rk4 or https://stackoverflow.com/questions/29481496, https://stackoverflow.com/questions/53645649/cannot-get-rk4 — Lutz Lehmann, Feb 11 '22 at 21:43
@Savannah According to the suggestions of Lutz Lehmann, I have proposed a solution. Please check it out and let me know if that solved your problem. — Muhammad Mohsin Khan, Feb 15 '22 at 08:27

Muhammad Mohsin Khan · Accepted Answer · 2022-02-24T08:05:46.790

As Lutz Lehmann pointed out, you need to design a simple ODE system. You could define the whole ODE system inside a function as follows:

import numpy as np

def fun(t, RHS):

    # get initial boundary condition values
    y1 = RHS[0]
    y2 = RHS[1]
    y3 = RHS[2]

    # calculte rate of respective variables
    dy1dt = (Q/V1)*(y3-y1) - y1/tau
    dy2dt = (Q/V2)*(y1-y2) + (qh/V2)*(y3-y2) - (fh/V2)
    dy3dt = (Q/V3)*(y2-y3) + (qh/V3)*(y2-y3) + (V1*y1)/(V3*tau) + (fh/V3)

    # Left-hand side of ODE
    LHS = np.zeros([3,])

    LHS[0] = dy1dt
    LHS[1] = dy2dt
    LHS[2] = dy3dt

    return LHS

In the above function, we get time t as an argument and the initial values of y1, y2, and y3 as a list in the variable RHS, which is then unpacked to get the respective variables. Afterward, the rate equation of each variable is defined. In the end, the calculated rates are returned also as a list in the variable LHS.

Now, we can define a simple Euler Forward method to solve this ODE system as follows:

def ForwardEuler(fun,t0,y0,N,dt):

    t = np.arange(t0,N+dt,dt)
    y = np.zeros([len(t), len(y0)])
    y[0] = y0

    for i in range(len(t)-1):
        y[i+1,:] = y[i,:] + fun(t[i],y[i,:]) * dt

    return t, y

Here, we create a time range from 0 to 100 with a step size of 0.1. Afterward, we create an array of zeros with the shape (len(t), len(y0)) which is in this case (1001,3). We need to do this because we want to solve fun for the range of t (1001) and the RHS variable of fun has a shape of (3,) ([y1, y2, y3]). So for each and every point in t, we will solve for the three variables of RHS, which will be returned as LHS.

In the end, we can solve this ODE system as follows:

dt = 0.1
N = 100
y0 = [0.00215, 0.00215, 0.00215]
t0 = 0

t,y = ForwardEuler(fun,t0,y0,N,dt)

Solution using scipy.integrate

As Lutz Lehmann also pointed out, you can use scipy.integrate for this purpose as well which is far easier. Here you can use the above defined fun and simply solve the ODE as follows:

import numpy as np
from scipy.integrate import odeint

dt = 0.1
N = 100
t = np.linspace(0,N,int(N/dt))
y0 = [0.00215, 0.00215, 0.00215]

res = odeint(fun, y0, t, tfirst=True)
print(res)

Applying Forward Euler Method to a Three-Box Model System of ODEs

1 Answers1

Solution using scipy.integrate

Linked