1

我正在尝试实现 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()

我刚刚在评论之后更新了代码。所以，功能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

然而，结果如下：

在此处输入图像描述

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

在此处输入图像描述

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

4

1 回答 1

1

您正在做一个非常典型的错误，例如如何通过 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

以避免这个问题。

于 2016-01-28T20:37:06.970 回答