在没有包的情况下使用 Python 中的 Runge Kutta 4th Order 求解 Lorentz 模型

Question

我希望在没有包帮助的情况下在 Python 中求解 Lorentz 模型，但我的代码似乎没有达到我的预期。我不知道为什么我没有得到预期的结果和洛伦兹吸引子。我猜的主要问题与如何分别存储 x、y 和 z 的解的各种值有关。以下是我的 Lorentz 模型的 Runge-Kutta 45 代码，带有 3D 解图：

import numpy as np
import matplotlib.pyplot as plt
#from scipy.integrate import odeint

#a) Defining the Runge-Kutta45 method

def fx(x,y,z,t):
    dxdt=sigma*(y-z)
    return dxdt

def fy(x,y,z,t):
    dydt=x*(rho-z)-y
    return dydt

def fz(x,y,z,t):
    dzdt=x*y-beta*z
    return dzdt


def RungeKutta45(x,y,z,fx,fy,fz,t,h):
    k1x,k1y,k1z=h*fx(x,y,z,t),h*fy(x,y,z,t),h*fz(x,y,z,t)
    k2x,k2y,k2z=h*fx(x+k1x/2,y+k1y/2,z+k1z/2,t+h/2),h*fy(x+k1x/2,y+k1y/2,z+k1z/2,t+h/2),h*fz(x+k1x/2,y+k1y/2,z+k1z/2,t+h/2)
    k3x,k3y,k3z=h*fx(x+k2x/2,y+k2y/2,z+k2z/2,t+h/2),h*fy(x+k2x/2,y+k2y/2,z+k2z/2,t+h/2),h*fz(x+k2x/2,y+k2y/2,z+k2z/2,t+h/2)
    k4x,k4y,k4z=h*fx(x+k3x,y+k3y,z+k3z,t+h),h*fy(x+k3x,y+k3y,z+k3z,t+h),h*fz(x+k3x,y+k3y,z+k3z,t+h)
    return x+(k1x+2*k2x+2*k3x+k4x)/6,y+(k1y+2*k2y+2*k3y+k4y)/6,z+(k1z+2*k2z+2*k3z+k4z)/6



sigma=10.
beta=8./3.
rho=28.
tIn=0.
tFin=10.
h=0.05
totalSteps=int(np.floor((tFin-tIn)/h))

t=np.zeros(totalSteps)
x=np.zeros(totalSteps) 
y=np.zeros(totalSteps) 
z=np.zeros(totalSteps)

for i in range(1, totalSteps):
    x[i-1]=1. #Initial condition
    y[i-1]=1. #Initial condition
    z[i-1]=1. #Initial condition
    t[0]=0. #Starting value of t
    t[i]=t[i-1]+h
    x,y,z=RungeKutta45(x,y,z,fx,fy,fz,t[i-1],h)


#Plotting solution

from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

fig=plt.figure()
ax=fig.gca(projection='3d')
ax.plot(x,y,z,'r',label='Lorentz 3D Solution')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.legend()

Answer 1

我更改了集成步骤（顺便说一句，经典的 4 阶 Runge-Kutta，而不是自适应 RK54）以广泛使用列表和列表操作的 Python 核心概念，以减少定义计算的位置数量。那里没有要纠正的错误，但我认为算法本身更集中。

您在系统中出现了一个错误，将其更改为一个快速发散的系统。你有fx = sigma*(y-z)，而洛伦兹系统有fx = sigma*(y-x)。

接下来你的主循环有一些奇怪的任务。在每个循环中，您首先将先前的坐标设置为初始条件，然后用应用于完整数组的 RK 步骤替换完整数组。我完全替换了它，正确的解决方案没有小步骤。

import numpy as np
import matplotlib.pyplot as plt
#from scipy.integrate import odeint


def fx(x,y,z,t): return sigma*(y-x)
def fy(x,y,z,t): return x*(rho-z)-y
def fz(x,y,z,t): return x*y-beta*z

#a) Defining the classical Runge-Kutta 4th order method

def RungeKutta45(x,y,z,fx,fy,fz,t,h):
    k1x,k1y,k1z = ( h*f(x,y,z,t) for f in (fx,fy,fz) )
    xs, ys,zs,ts = ( r+0.5*kr for r,kr in zip((x,y,z,t),(k1x,k1y,k1z,h)) )
    k2x,k2y,k2z = ( h*f(xs,ys,zs,ts) for f in (fx,fy,fz) )
    xs, ys,zs,ts = ( r+0.5*kr for r,kr in zip((x,y,z,t),(k2x,k2y,k2z,h)) )
    k3x,k3y,k3z = ( h*f(xs,ys,zs,ts) for f in (fx,fy,fz) )
    xs, ys,zs,ts = ( r+kr for r,kr in zip((x,y,z,t),(k3x,k3y,k3z,h)) )
    k4x,k4y,k4z  =( h*f(xs,ys,zs,ts) for f in (fx,fy,fz) )
    return (r+(k1r+2*k2r+2*k3r+k4r)/6 for r,k1r,k2r,k3r,k4r in 
            zip((x,y,z),(k1x,k1y,k1z),(k2x,k2y,k2z),(k3x,k3y,k3z),(k4x,k4y,k4z)))



sigma=10.
beta=8./3.
rho=28.
tIn=0.
tFin=10.
h=0.01
totalSteps=int(np.floor((tFin-tIn)/h))

t = totalSteps * [0.0]
x = totalSteps * [0.0]
y = totalSteps * [0.0]
z = totalSteps * [0.0]

x[0],y[0],z[0],t[0] = 1., 1., 1., 0.  #Initial condition
for i in range(1, totalSteps):
    x[i],y[i],z[i] = RungeKutta45(x[i-1],y[i-1],z[i-1], fx,fy,fz, t[i-1], h)

使用tFin = 40 和h=0.01 我得到图像

看起来像洛伦兹吸引子的典型图像。