洛伦兹吸引子与Runge-Kutta蟒蛇的关系。

Question

1. 你好，我必须用Runge-Kutta 2cond grade编程一个python函数来解决Lorenz微分方程。

sigma=10，r=28，b=83。

初始条件为(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)，或者直接使用Runge Kutta提供的自适应Runge Kutta方法。scipy:

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

该求解器使用4阶和5阶Runge Kutta组合，并通过调整步长来控制它们之间的偏差。更多使用文档请看这里。https:/docs.scipy.orgdocscipyreferencegeneratedscipy.integration.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见 ODE解决方案的混沌和连续依赖性所以 h<0.05 是绝对需要的。h=0.002 是更好的和 h=2e-5 给出了这种数值方法的最佳结果。