1
  1. 您好，我必须编写一个 python 函数来使用 Runge-Kutta 2cond 等级求解 Lorenz 微分方程

微分方程

sigma=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)
4

3 回答 3

6

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

您应该大幅减小步长 ( h)，或者只使用由 提供的自适应龙格库塔方法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.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

在此处输入图像描述

于 2020-06-14T19:46:28.693 回答
1

您使用的步长为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为这种数值方法提供数值最佳结果。

于 2020-06-15T08:02:47.503 回答
0

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

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

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

于 2020-06-14T19:56:33.313 回答