带 Runge-Kutta 的洛伦兹吸引子 python

Question

你好，我必须编写一个 python 函数来使用 Runge-Kutta 二级

西格玛=10，r=28 和 b=8/3

初始条件 (x,y,z)=(0,1,0)

这是我写的代码，但它抛出一个错误说 double_scalars 中遇到溢出，而且我看不出程序有什么问题

from pylab import *
def runge_4(r0,a,b,n,f1,f2,f3):
    def f(r,t):
        x=r[0]
        y=r[1]
        z=r[2]
        fx=f1(x,y,z,t)
        fy=f2(x,y,z,t)
        fz=f3(x,y,z,t)
        return array([fx,fy,fz],float)
    h=(b-a)/n
    lista_t=arange(a,b,h)
    print(lista_t)
    X,Y,Z=[],[],[]
    for t in lista_t:
        k1=h*f(r0,t)
        print("k1=",k1)
        k2=h*f(r0+0.5*k1,t+0.5*h)
        print("k2=",k2)
        k3=h*f(r0+0.5*k2,t+0.5*h)
        print("k3=",k3)
        k4=h*f(r0+k3,t+h)
        print("k4=",k4)
        r0+=(k1+2*k2+2*k3+k4)/float(6)
        print(r0)
        X.append(r0[0])
        Y.append(r0[1])
        Z.append(r0[2])
    return array([X,Y,Z])

def f1(x,y,z,t):
    return 10*(y-x)
def f2(x,y,z,t):
    return 28*x-y-x*z
def f3(x,y,z,t):
    return x*y-(8.0/3.0)*z
#and I run it
r0=[1,1,1]

runge_4(r0,1,50,20,f1,f2,f3)

Answer 1

以数值方式求解微分方程可能具有挑战性。如果您选择的步长过大，解决方案将累积大量错误，甚至会变得不稳定，就像您的情况一样。

要么大幅减小步长 (h)，要么只使用 scipy:

提供的自适应 Runge Kutta 方法

from numpy import array, linspace
from scipy.integrate import solve_ivp
import pylab
from mpl_toolkits import mplot3d

def func(t, r):
    x, y, z = r 
    fx = 10 * (y - x)
    fy = 28 * x - y - x * z
    fz = x * y - (8.0 / 3.0) * z
    return array([fx, fy, fz], float)

# and I run it
r0 = [0, 1, 0]
sol = solve_ivp(func, [0, 50], r0, t_eval=linspace(0, 50, 5000))

# and plot it
fig = pylab.figure()
ax = pylab.axes(projection="3d")
ax.plot3D(sol.y[0,:], sol.y[1,:], sol.y[2,:], 'blue')
pylab.show()

该求解器使用四阶和五阶龙格库塔组合，并通过调整步长来控制它们之间的偏差。在此处查看更多使用文档：https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

Answer 2

它可能与被零除或超出类型限制（浮点类型）有关。

要弄清楚它发生的时间和地点，您可以设置 numpy.seterr('raise')，它会引发异常，以便您调试并查看发生了什么。看来你的算法发散了。

在这里你可以了解如何使用numpy.seterr

Answer 3

您使用 h=2.5 的步长。

对于 RK4，给定 Lipschitz 常数 L 的有用步长在 L*h=1e-3 到 0.1 范围内，最多 L*h=2.5 可能会得到一些看起来不错的结果.在此之上，该方法变得混乱，与底层 ODE 的任何相似之处都丢失了。

Lorenz 系统的 Lipschitz 常数约为 L=50，参见 Chaos and continuous dependency of ODE solution，因此 h<0.05 是绝对必需的，h=0.002 更好，h=2e-5给出了该数值方法的最佳数值结果。

带 Runge-Kutta 的洛伦兹吸引子 python

Lorenz attractor with Runge-Kutta python

python

differential-equations

runge-kutta

lorenz-system