Runge-Kutta 代码不与内置方法收敛

Question

我正在尝试实施 runge-kutta 方法来解决 Lotka-Volterra 系统，但代码（波纹管）无法正常工作。我遵循了我在 StackOverflow 的其他主题中找到的建议，但结果与内置的 Runge-Kutta 方法不收敛，例如 Pylab 中可用的 rk4 方法。有人可以帮助我吗？

import matplotlib.pyplot as plt
import numpy as np
from pylab import *

def meurk4( f, x0, t ):
    n = len( t )
    x = np.array( [ x0 ] * n )    

    for i in range( n - 1 ):

        h =  t[i+1] - t[i]

        k1 = h * f( x[i], t[i] )
        k2 = h * f( x[i] + 0.5 * h * k1, t[i] + 0.5 * h )
        k3 = h * f( x[i] + 0.5 * h * k2, t[i] + 0.5 * h )
        k4 = h * f( x[i] + h * k3, t[i] + h)

        x[i+1] = x[i] + ( k1 + 2 * ( k2 + k3 ) + k4 ) * 6**-1 

    return x

def model(state,t):

    x,y = state     

    a = 0.8
    b = 0.02
    c = 0.2
    d = 0.004
    k = 600

    return np.array([ x*(a*(1-x*k**-1)-b*y) , -y*(c - d*x) ]) # corresponds to [dx/dt, dy/dt]

# initial conditions for the system
x0 = 500
y0 = 200

# vector of time
t = np.linspace( 0, 50, 100 )

result = meurk4( model, [x0,y0], t )
print result

plt.plot(t,result)

plt.xlabel('Time')
plt.ylabel('Population Size')
plt.legend(('x (prey)','y (predator)'))
plt.title('Lotka-Volterra Model')
plt.show()

---

我刚刚更新了 cmets 之后的代码。所以，函数meurk4：

def meurk4( f, x0, t ):
        n = len( t )
        x = np.array( [ x0 ] * n )    

        for i in range( n - 1 ):

            h =  t[i+1] - t[i]

            k1 = h * f( x[i], t[i] )
            k2 = h * f( x[i] + 0.5 * h * k1, t[i] + 0.5 * h )
            k3 = h * f( x[i] + 0.5 * h * k2, t[i] + 0.5 * h )
            k4 = h * f( x[i] + h * k3, t[i] + h)

            x[i+1] = x[i] + ( k1 + 2 * ( k2 + k3 ) + k4 ) * 6**-1 

        return x

现在变成（更正）：

def meurk4( f, x0, t ):
    n = len( t )
    x = np.array( [ x0 ] * n )    

    for i in range( n - 1 ):

        h =  t[i+1] - t[i]

        k1 = f( x[i], t[i] )
        k2 = f( x[i] + 0.5 * h * k1, t[i] + 0.5 * h )
        k3 = f( x[i] + 0.5 * h * k2, t[i] + 0.5 * h )
        k4 = f( x[i] + h * k3, t[i] + h)

        x[i+1] = x[i] + ( k1 + 2 * ( k2 + k3 ) + k4 ) * (h/6)

    return x

不过，结果如下：

enter image description here

而 buitin 方法 rk4（来自 Pylab）的结果如下：

enter image description here

所以，当然我的代码仍然不正确，因为它的结果与内置的 rk4 方法不同。请问有人可以帮我吗？

Answer 1

您正在犯一个非常典型的错误，例如见How to pass a hard coded differential equation through Runge-Kutta 4 或这里Error in RK4 algorithm in Python

要么是

k2 = f( x+0.5*h*k1, t+0.5*h )
...
x[i+1]=x[i]+(k1+2*(k2+k3)+k4)*(h/6)

或

k2 = h*f( x+0.5*k1, t+0.5*h )

以此类推，x[i+1] 原样，但不能同时使用两种变体。

---

更新：一个更隐蔽的错误是初始值的推断类型以及x 向量数组的结果。根据原始定义，两者都是整数，因此

x = np.array( [ x0 ] * n )    

创建一个整数向量列表。因此更新步骤

    x[i+1] = x[i] + ( k1 + 2 * ( k2 + k3 ) + k4 ) * (h/6)

将始终舍入为整数。由于存在两个值都低于1 的阶段，因此积分稳定在零。从而修改为

# initial conditions for the system
x0 = 500.0
y0 = 200.0

为了避免这个问题。