Runge Kutta 常数对于 Lorenz 系统发散？

Question

我正在尝试使用四阶 Runge Kutta 方法求解 Lorenz system，其中

dx/dt=a*(y-x)
dy/dt=x(b-z)-y
dx/dt=x*y-c*z

由于这个系统不明确依赖于时间，它可能会在迭代中忽略那部分，所以我只需要 dX=F(x,y,z)

def func(x0):
    a=10
    b=38.63
    c=8/3
    fx=a*(x0[1]-x0[0])
    fy=x0[0]*(b-x0[2])-x0[1]
    fz=x0[0]*x0[1]-c*x0[2]
    return np.array([fx,fy,fz])

def kcontants(f,h,x0):
    k0=h*f(x0)
    k1=h*f(f(x0)+k0/2)
    k2=h*f(f(x0)+k1/2)
    k3=h*f(f(x0)+k2)
    #note returned K is a matrix
    return np.array([k0,k1,k2,k3])

x0=np.array([-8,8,27])
h=0.001

t=np.arange(0,50,h)
result=np.zeros([len(t),3])

for time in range(len(t)):
    if time==0:
        k=kcontants(func,h,x0)
        result[time]=func(x0)+(1/6)*(k[0]+2*k[1]+2*k[2]+k[3])
    else:
        k=kcontants(func,h,result[time-1])
        result[time]=result[time-1]+(1/6)*(k[0]+2*k[1]+2*k[2]+k[3])

结果应该是 Lorenz 吸引子，但是我的代码在第五次迭代时出现分歧，这是因为我在 kconstants 中创建的内容确实如此，但是我检查了一下，我很确定 runge kutta 实现不是错...（至少我认为）

编辑：

找到一个类似的，但无法弄清楚我做错了什么

Answer 1

您是否查看过不同的计算初始值？你选择的那些有意义吗？ IE。他们是身体上的吗？根据过去使用 rk 的经验，如果您选择愚蠢的起始参数，有时会得到非常混乱的结果。

Answer 2

您在k1、k2和k3的计算中额外调用了f(x0)。将函数 kcontants 更改为

def kcontants(f,h,x0):
    k0=h*f(x0)
    k1=h*f(x0 + k0/2)
    k2=h*f(x0 + k1/2)
    k3=h*f(x0 + k2)
    #note returned K is a matrix
    return np.array([k0,k1,k2,k3])

Runge Kutta 常数对于 Lorenz 系统发散？

Runge Kutta constants diverging for Lorenz system?

python

numerical-methods

runge-kutta

lorenz-system

编辑：