将正向欧拉方法应用于常微分方程的三盒模型系统

Question

我正在尝试建立耦合 ODE 系统的模型，该系统表示低纬度表层海洋（方框 1）、高纬度深海（方框 2）中磷浓度 (y) 的三方框海洋模型和深海（专栏 3）。 ODE 如下：

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

常量和盒子体积由下式给出：

### 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

这可以写成矩阵形式y = Ay + f, 其中y = [y1, y2, y3]。我在下面提供了矩阵和初始条件：

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

我在调整正向欧拉方法以应用于基于线性 ODE 的系统时遇到问题，而不仅仅是一个。这是我到目前为止的想法（它运行没有问题，但如果有意义的话就不起作用；我认为它与初始条件无关？）：

### 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

非常感谢对此的任何帮助。我还没有在网上找到任何通过 Euler Forward 获得 3 个 ODE 系统的资源，而且我很茫然。如果有更多问题，我很乐意进一步解释。

Answer 1

正如 Lutz Lehmann 所指出的，您需要设计一个简单的 ODE 系统。您可以在一个函数内定义整个 ODE 系统，如下所示：

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

在上面的函数中，我们将时间t作为参数，并将y1、y2和y3的初始值作为变量中的列表RHS，然后解包得到相应的变量。之后，定义了每个变量的速率方程。最后，计算出的费率也作为变量 LHS.

中的列表返回

现在，我们可以定义一个简单的 Euler Forward 方法来求解这个 ODE 系统，如下所示：

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

在这里，我们创建了一个从 0 到 100 的时间范围，步长为 0.1。之后，我们创建一个形状为 (len(t), len(y0)) 的零数组，在本例中为 (1001,3)。我们需要这样做是因为我们要解决 fun 范围 t (1001) 并且 fun 的 RHS 变量具有 (3,) 的形状([y1, y2, y3])。因此，对于 t 中的每个点，我们将求解 RHS 的三个变量，它们将返回为 LHS.

最后，我们可以如下求解这个ODE系统：

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

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

解决方案使用 scipy.integrate

正如 Lutz Lehmann 也指出的那样，您也可以使用 scipy.integrate 来达到此目的，这要容易得多。在这里你可以使用上面定义的 fun 并简单地解决 ODE 如下：

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

python

ode

runge-kutta

解决方案使用 scipy.integrate