问如何用龙格-库塔法求解Lorenz 96模型？

EN

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