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

Posted 2023-03-24

【中文标题】Runge-Kutta 代码不与内置方法收敛

【英文标题】：Runge-Kutta code not converging with builtin method

【发布时间】：2016-05-06 09:55:12

【问题描述】：

我正在尝试实施 runge-kutta 方法来解决 Lotka-Volterra 系统，但代码（波纹管）无法正常工作。我遵循了我在 *** 的其他主题中找到的建议，但结果与内置的 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 方法不同。请问有人可以帮我吗？

【问题讨论】：

通过该修改，变体 1，您需要将 x[i+1] 中的因子修改为 (h/6)。

你是对的，我的错。尽管如此，我只是实现了这一点，它与 Pylab 内置方法仍然存在差异。

你确定常量和初始值是一样的吗？

我刚刚在代码上将result = meurk4( model, [x0,y0], t ) 更改为result = rk4( model, [x0,y0], t )。

只能猜测 pylab/matplotlib rk4 的内部结构。我假设它具有某种步长控制。步长 0.5 太大。 0.01 或 0.001 是合适的。但似乎没有帮助。

【参考方案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

为了避免这个问题。

【讨论】：

感谢您纠正我在 numpy 中使用 x.* 的错误，正如您建议的那样，我正在考虑 Matlab。无论如何，你的回答是对的，我的大脑因一整天的工作而煎熬。

我刚刚删除了 k 之前的 h 因子，即变体 2，程序运行。即使收敛为零。 （也许选择参数使平衡点位于第一象限） -- 请在问题中添加一个部分，说明您更正了什么以及新错误的描述。

感谢 LutzL。我刚刚用一个新部分更新了这个问题，并根据您的建议修改了代码。