如何用龙格-库塔法求解洛伦兹96模型？

Question

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

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

如何修改代码以正确解决问题？有关Lorenz 96的更多信息，请参见Lorenz 96 model- Wikipedia

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
   
###############################################################################

# Define Lorenz 96 function

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


def L96(v, t):
    """Lorenz 96 model with constant forcing"""
    # Setting up vector
    dv_dt = np.zeros(N)
    # 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_range=np.arange(t0, tn, h)

# preallocate space
v0=np.zeros((len(time_range)+1, 1, N))
t=np.zeros(len(time_range)+1)

v0[0][0][0] += 0.01  # Add small perturbation to the first variable
L96(v0[0][0],0)

# Solve Runge-Kutta
for i in range (len(time_range)):
    print(i)
    dv_dt=L96(v0[i][0],t[i])
    
    k1=dv_dt[0]
    l1=dv_dt[1]
    m1=dv_dt[2]
    
    k2=L96([v0[i][0][0]+0.5*k1*h,v0[i][0][1]+0.5*l1*h,v0[i][0][2]+0.5*m1*h],t[i]+h/2)
    l2=L96([v0[i][0][0]+0.5*k1*h,v0[i][0][1]+0.5*l1*h,v0[i][0][2]+0.5*m1*h],t[i]+h/2)
    m2=L96([v0[i][0][0]+0.5*k1*h,v0[i][0][1]+0.5*l1*h,v0[i][0][2]+0.5*m1*h],t[i]+h/2)
  
    k3=L96([v0[i][0][0]+0.5*k2[0]*h,v0[i][0][1]+0.5*l2[0]*h,v0[i][0][2]+0.5*m2[0]*h],t[i]+h/2)
    l3=L96([v0[i][0][0]+0.5*k2[1]*h,v0[i][0][1]+0.5*l2[1]*h,v0[i][0][2]+0.5*m2[1]*h],t[i]+h/2)
    m3=L96([v0[i][0][0]+0.5*k2[2]*h,v0[i][0][1]+0.5*l2[2]*h,v0[i][0][2]+0.5*m2[2]*h],t[i]+h/2)
  
    k4=L96([v0[i][0][0]+0.5*k3[0]*h,v0[i][0][1]+0.5*l3[0]*h,v0[i][0][2]+0.5*m3[0]*h],t[i]+h/2)
    l4=L96([v0[i][0][0]+0.5*k3[1]*h,v0[i][0][1]+0.5*l3[1]*h,v0[i][0][2]+0.5*m3[1]*h],t[i]+h/2)
    m4=L96([v0[i][0][0]+0.5*k3[2]*h,v0[i][0][1]+0.5*l3[2]*h,v0[i][0][2]+0.5*m3[2]*h],t[i]+h/2)
  
    v0[i+1][0][0] = v0[i][0][0] + h*(k1 +2*k2[0]  +2*k3[0]   +k4[0])/6      # final equations
    v0[i+1][0][1] = v0[i][0][1] + h*(l1  +2*k2[1]   +2*k3[1]    +k4[1])/6
    v0[i+1][0][2] = v0[i][0][2] + h*(m1+2*k2[2] +2*k3[2]  +k4[2])/6
    t[i+1]=time_range[i]

###############################################################################

v_array=np.array(v0)
v_array.shape

v1=v_array[:,0][:,0]
v2=v_array[:,0][:,1]
v3=v_array[:,0][:,2]


fig, (ax1, ax2) = plt.subplots(1, 2,constrained_layout=True,figsize=(10, 4))  
ax1.plot(v1,v3) 
ax1.set_title('v1 vs v2')    
ax2.plot(v2,v3)
ax2.set_title('v2 vs v3')    
    
# Plot 3d plot of v1,v2, and v3
from mpl_toolkits import mplot3d

fig = plt.figure(figsize=(8, 5))
ax = plt.axes(projection='3d', elev=40, azim=-50)
ax.plot3D(v1, v2, v3)
ax.set_xlabel('$v1$')
ax.set_ylabel('$v2$')
ax.set_zlabel('$v3$')
plt.show()

Answer 1

我真的不确定你的动机是什么，把矢量化再次拆分成单个组件的行。无论如何，这是错误的。你的时间循环应该尽可能简单

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])

其中，stepper方法只是烹饪书代码

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的步进器和时间循环适用于任何其他微分方程系统

请注意，在wikipedia示例中，有N=5个组件，其中3D绘图是由前3个组件构建的。对于N=3你确实会收敛到一个看起来像一条直线的固定点，对于N=4它似乎有一个极限环，只有对于N=5解才开始看起来混乱

如果同时计算多条轨迹，拥有一个三维阵列是有意义的。（但只有使用固定的步长方法，或者如果保证它们保持在一起，则步长控制在大多数情况下不适用于大多数轨迹。）该调整的完整脚本如下

# 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] -= 0.005
v[0,0,1] += 0.005
v[0,0,2] += 0.01
for i in range(len(time)-1):
    v[i+1] = RK4_step(L96,time[i],v[i],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,0], v[:,1,0], v[:,2,0],'g', lw=0.2)
ax.plot3D(v[:,0,1], v[:,1,1], v[:,2,1],'r', lw=0.2)
ax.plot3D(v[:,0,2], v[:,1,2], v[:,2,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()