欧拉法在三箱模型系统中的应用
我试图模拟一个耦合的海洋系统，它代表了低纬度表层海洋(方框1)、高纬度深海(框2)和深海(框3)中磷浓度(y)的三盒海洋模型。以下是颂歌：
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我很难将正向欧拉方法应用于线性ODEs系统，而不仅仅是一个。到目前为止，这就是我想出来的(它运行时没有问题，但如果这有意义的话就不起作用了；我认为它与初始条件无关吗？)：
### 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在这方面的任何帮助都是非常感谢的。我没有在网上找到任何资源，通过欧拉转发系统的3首歌，并在一个不知所措。如果还有更多的问题，我很乐意进一步解释。
回答 1
Stack Overflow用户
发布于 2022-02-14 10:45:03
正如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转发方法来解决这个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)。我们需要这样做，因为我们希望为t (1001)的范围求解t，而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)https://stackoverflow.com/questions/71080688
复制相似问题
腾讯云开发者
