且构网

分享程序员开发的那些事...
且构网 - 分享程序员编程开发的那些事

不带包装的Python中使用Runge Kutta 4th Order解决Lorentz模型

更新时间:2022-06-26 23:59:02

我更改了集成步骤(顺便说一句,经典的四阶Runge-Kutta,而不是自适应RK54),以广泛使用python核心概念的列表和列表操作以减少定义计算的位置数量.没有错误可以纠正,但是我认为算法本身更加集中.

I changed the integration step (btw., classical 4th order Runge-Kutta, not the adaptive RK54) to use the python core concept of lists and list operations extensively to reduce the number of places where the computation is defined. There were no errors there to correct, but I think the algorithm itself is more concentrated.

您的系统中出现错误,将其更改为快速发散的系统. Lorentz系统具有fx = sigma*(y-x)时,您拥有fx = sigma*(y-z).

You had an error in the system that changed it into a system that rapidly diverges. You had fx = sigma*(y-z) while the Lorentz system has fx = sigma*(y-x).

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

Next your main loop has some strange assignments. In every loop you first set the previous coordinates to the initial conditions and then replace the full arrays with the RK step applied to the full arrays. I replaced that completely, there are no small steps to a correct solution.

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 = 40h=0.01我得到了图像

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

looking like the typical image of the Lorentz attractor.