需要帮助修复我的 RK4 实施

Question

如果在实施方面更有经验的人能帮助我发现我当前代码中的逻辑缺陷，我将不胜感激。在过去的几个小时里，我一直在执行和测试以下 RK4 函数的各种步长来解决 Lotka-Volterra Differential equation.

我尽最大努力确保代码的可读性并注释掉关键步骤，因此下面的代码应该清晰。

import matplotlib.pyplot as plt
import numpy as np

def model(state,t):
    """
    A function that creates an 1x2-array containing the Lotka Volterra Differential equation

    Parameter assignement/convention:
    a   natural growth rate of the preys
    b   chance of being eaten by a predator
    c   dying rate of the predators per week
    d   chance of catching a prey 
    """

    x,y = state     # will corresponding to initial conditions  
                    # consider it as a vector too 

    a = 0.08
    b = 0.002
    c = 0.2
    d = 0.0004

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


def rk4( f, x0, t):
    """ 
    4th order Runge-Kutta method implementation to solve x' = f(x,t) with x(t[0]) = x0.

    INPUT:
        f     - function of x and t equal to dx/dt. 

        x0    - the initial condition(s).  
                Specifies the value of x @ t = t[0] (initial).  
                Can be a scalar or a vector (NumPy Array)

                Example: [x0, y0] = [500, 20]

        t     - a time vector (array) at which the values of the solution are computed at.
                t[0] is considered as the initial time point 
                the step size h is dependent on the time vector, choosing more points will 
                result in a smaller step size. 


    OUTPUT:
        x     - An array containing the solution evaluated at each point in the t array.

    """

    n = len( t )
    x = np.array( [ x0 ] * n )      # creating an array of length n 

    for i in xrange( n - 1 ):
        h = t[i+1]- t[i]            # step size, dependent on time vector

        # starting below - the implementation of the RK4 algorithm:
        # for further informations visit http://en.wikipedia.org/wiki/Runge-Kutta_methods

        # k1 is the increment based on the slope at the beginning of the interval (same as Euler)
        # k2 is the increment based on the slope at the midpoint of the interval 
        # k3 is AGAIN the increment based on the slope at the midpoint 
        # k4 is the increment based on the slope at the end of the interval

        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  )

        # finally computing the weighted average and storing it in the x-array
        t[i+1] = t[i] + h
        x[i+1] = x[i] + h * ( ( k1 + 2.0 * ( k2 + k3 ) + k4 ) / 6.0 )

    return x


################################################################
# just the graphical output

# initial conditions for the system

x0 = 500
y0 = 20


# vector of times
t = np.linspace( 0, 200, 150 )

result = rk4( model,[x0,y0], t )


plt.plot(t,result)

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

当前输出看起来 'okay-ish' 一小段时间然后变为 'berserk'。奇怪的是，当我选择较大的步长而不是较小的步长时，代码似乎表现得更好，这表明我的实现一定是错误的，或者我的模型可能已关闭。我自己无法发现错误。

输出（错误）：

这是所需的输出，可以通过使用 Scipys 集成模块之一轻松获得。请注意，在时间间隔 [0,50] 上，模拟似乎是正确的，然后每一步都变得更糟。

Answer 1

不幸的是，您陷入了我偶尔也会陷入的同一个陷阱：您的初始 x0 数组包含整数，因此，所有结果 x[i] 值将在计算后转换为整数.

这是为什么？因为int是你初始条件的类型：

x0 = 500
y0 = 20

解决方案当然是明确地使它们 floats:

x0 = 500.
y0 = 20.

那么为什么 scipy 在您输入整数起始值时会正确执行呢？在开始实际计算之前，它可能会将它们转换为 float。例如，您可以这样做：

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

然后您仍然可以安全地使用整数初始条件而不会出现问题。至少这样一劳永逸地在函数内部完成了转换，半年以后你再用一次（或者，别人用过），就不会再落入那个圈套了。

需要帮助修复我的 RK4 实施

Need help fixing my implementation of RK4

python

python-2.7

debugging

numerical-methods

runge-kutta