如何使用 Runge-Kutta 方法求解 Lorenz 96 模型？

问题描述

我已经编写了以下代码来使用 Runge-Kutta 方法求解 Lorenz 96 模型，但结果并不可靠，如下图所示：

三个变量的正确关系如下：

如何修改代码才能正确解决问题？有关洛伦兹 96 的更多信息，请参阅 Lorenz 96 model- Wikipedia

this.cardindexShow

解决方法

我不太确定您将矢量化再次分解为单个组件的行的动机是什么。无论如何，这是错误的。你的时间循环应该像

一样简单

for i in range(len(time)-1):
    v[i+1] = RK4_step(lambda t,y:L96(y,t),time[i],v[i],time[i+1]-time[i])

方法步进器只是食谱代码

def RK4_step(f,t,y,h):
    k1 = h*f(t,y)
    k2 = h*f(t+0.5*h,y+0.5*k1)
    k3 = h*f(t+0.5*h,y+0.5*k2)
    k4 = h*f(t+h,y+k3)
    return y+(k1+2*k2+2*k3+k4)/6

因此对 3 个组件没有硬编码限制，代码组件是通用的，L96 函数可以与库中的任意数量的数值求解器一起使用或自行实现，步进器和时间循环用于RK4 对任何其他微分方程组都有效。

请注意，在维基百科示例中，有 N=5 组件，其中 3D 图是从其中的前 3 个构建的。对于 N=3，您确实收敛到一个看起来像直线的固定点，对于 N=4，它看起来好像有一个极限环，只有对于 N=5，解才开始看起来混乱。

如果一次计算多个轨迹，则拥有一个 3 维数组是有意义的。 （但只有使用固定步长的方法，或者如果保证它们保持靠近，否则大多数情况下步长控制将不适合大多数轨迹。）该适应的完整脚本读取

# Define Lorenz 96 function

# These are our constants
N = 4  # Number of variables
F = 8  # Forcing


def L96(t,v):
    """Lorenz 96 model with constant forcing"""
    # Setting up vector
    dv_dt = 0*v
    # Loops over indices (with operations and Python underflow indexing handling edge cases)
    for i in range(N):
        dv_dt[i] = (v[(i + 1) % N] - v[i - 2]) * v[i - 1] - v[i] + F
    return dv_dt
   

#define the given range t
t0=0
tn=100
h=0.005

#define number of steps (n)

time = np.arange(t0,tn,h)
v = np.zeros([len(time),N,3],float)
v[0] +=F
v[0,0] -= 0.005
v[0,1] += 0.005
v[0,2] += 0.01
for i in range(len(time)-1):
    v[i+1] = RK4_step(L96,time[i+1]-time[i])

# Plot 3d plot of first 3 coordinates
from mpl_toolkits import mplot3d

fig = plt.figure(figsize=(8,5))
ax = plt.axes(projection='3d',elev=40,azim=-50)
ax.plot3D(v[:,0],v[:,1,2,'g',lw=0.2)
ax.plot3D(v[:,1],'r',2],'b',lw=0.2)
ax.set_xlabel('$v_1$')
ax.set_ylabel('$v_2$')
ax.set_zlabel('$v_3$')
ax.set_title(f"N={N}")
plt.show()