0

我试图了解我一起屠杀的代码有什么问题。下面的代码是我今天为求解 Lotka Volterra 微分方程（2 个系统）所做的众多实现之一，它是我带来的最接近预期结果的一个。

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



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

    USE:
        x = rk4(f, x0, t)

    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 of 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 
                h = t[i+1] - t[i] determines the step size h as suggested by the algorithm
                
                Example: t = np.linspace( 0, 500, 200 ), creates 200 time points between 0 and 500
                increasing the number of points in the intervall automatically decreases the 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 the 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 (with x + 0.5 * k1)
        # k3 is AGAIN the increment based on the slope at the midpoint (with x + 0.5 * k2)
        # 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 * k1, t[i] + 0.5 * h )
        k3 = f( x[i] + 0.5 * 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

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

    return x

 
# model
 
def model(state,t):
    """
    A function that creates an 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]
 


################################################################




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


# vector of times
t = np.linspace( 0, 500, 1000 )

result = rk4( 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()

上面的代码产生以下输出

在此处输入图像描述

但是，如果我将from pylab import *代码移到初始条件的正上方，我会得到正确的输出

在此处输入图像描述

为什么会发生这种情况，我该如何解决？

4

2 回答 2

1

pylab定义了它自己的实现rk4，它取自matplotlib：

In [1]: import pylab

In [2]: pylab.rk4
Out[2]: <function matplotlib.mlab.rk4>

当您执行通配符导入from pylab import *时，您将覆盖任何具有相同名称的本地函数。特别是，您在这里重新定义了自己的rk4实现（即，您编写的代码从未使用过）。

这就是为什么你不应该做这样的通配符导入。 pylab尤其成问题，因为它定义了几个函数（例如any和all），它们的输出与某些输入的 python 内置函数完全不同。

无论如何，您的问题的根本原因似乎是您的 RK4 实现不正确。

您需要在计算中使用步长k_n。例如，这是我自己的 RK4 实现的一小段（我承认，它是针对速度而不是可读性进行调整的）：

    while not target(xs):
        ...

        # Do RK4
        self.f(xs, self.k1)
        self.f(xs + halfh*self.k1, self.k2)
        self.f(xs + halfh*self.k2, self.k3)
        self.f(xs + self.h*self.k3, self.k4)

        xs += sixthh*(self.k1 + self.k2 + self.k2 + self.k3 + self.k3 \
                + self.k4)

您会注意到整个状态向量乘以h，而不仅仅是时间分量。尝试在您自己的代码中修复它并查看结果是否相同。

（在我看来，wiki 等将时间视为特殊情况的习惯是导致许多此类问题的原因。您的时间向量 ，ts，只是一个特殊的导数，其中t' = 1。）

因此，对于您自己的代码，我相信但尚未测试过，这样的东西应该可以工作：

    k1 = f( x[i], t[i] )
    k2 = f( x[i] + 0.5 * h * k1, t[i] + 0.5 * h )  ## changed to use h
    k3 = f( x[i] + 0.5 * h * k2, t[i] + 0.5 * h )  ## changed to use h
    k4 = f( x[i] + h * k3, t[i] + h )
于 2015-01-06T00:59:51.473 回答
0

尝试

import pylab
help(pylab.rk4)

您会发现 pylab.rk4 命令的详细说明。

这就是为什么使用from x import *. 这样做要好得多import pylab as py ，然后这将不是问题。

请注意，即使移动您的导入命令，您以后可能必须进行的任何调用都rk4将失败。

于 2015-01-06T01:00:29.057 回答