python中的Runge-Kutta方法

Question

以下是我编写的函数：

def rk4(f, t0, y0, h, N):
    t = t0 + arange(N+1)*h
    y = zeros((N+1, size(y0)))
    y[0] = y0
    for n in range(N):
        xi1 = y[n]
        f1 = f(t[n], xi1)
        xi2 = y[n] + (h/2.)*f1
        f2 = f(t[n+1], xi2)
        xi3 = y[n] + (h/2.)*f2
        f3 = f(t[n+1], xi3)
        xi4 = y[n] + h*f3
        f4 = f(t[n+1], xi4)
        y[n+1] = y[n] + (h/6.)*(f1 + 2*f2 + 2*f3 + f4)
    return y

def lorenzPlot():
    y = rk4(fLorenz, 0, array([0,2,20]), .01, 10000)
    fig = figure()
    ax = Axes3D(fig)
    ax.plot(*y.T)

def fLorenz(t,Y):   
    def Y(x,y,z):
        return array([10*(y-x), x*(28-z)-y, x*y-(8./3.)*z])

通过对lorenzPlot的实现，将用rk4（四阶Runge-Kutta方法）得到的fLorenz（Lorenz方程组）的数值解用图形表示出来。我对弗洛伦兹的功能有问题。当我把它定义为上面的并调用lorenzPlot时，我得到一个错误，它说

File "C:/...", line 38, in rk4
    xi2 = y[n] + (h/2.)*f1
TypeError: unsupported operand type(s) for *: 'float' and 'NoneType'

我猜这与不能正确地乘法数组有关。
但是，当我把弗洛伦兹换成

def fLorenz(t,x,y,z):   
    return array([10*(y-x), x*(28-z)-y, x*y-(8./3.)*z])

调用lorenzPlot会给我一个错误，指出fLorenz需要4个参数，但只给出了2个。

此外，rk4和lorenzPlot对于由奇异方程组成的函数都是正确的。

我应该如何改变fLorenz，使它可以作为f在rk4和lorenzPlot中使用？

Answer 1

第一个fLorenz函数定义了一个子函数，但实际上没有返回任何内容。你的第二个需要，除了它需要四个参数（t, x, y, z），你只需要给它t, Y。

据我所知，Y是一个三元组；您只需在使用值之前将其解包即可：

def fLorenz(t, Y):
    x, y, z = Y
    return array([10*(y-x), x*(28-z)-y, x*y-(8./3.)*z])

Answer 2

您的fLorenz不返回任何内容，也就是说，由于有很好的文档记录的Python约定，它返回None。

你的错误信息告诉你

unsupported operand type(s) for *: 'float' and 'NoneType'

错误发生在

    xi2 = y[n] + (h/2.)*f1

在前一行中，您分配给f1调用fLorenz返回的内容，即None。

如果更改函数fLorenz使其返回数值（标量或数组），则可以继续计算。

Answer 3

除了实现错误之外，您对RK4方法的理解是不完整的。中点斜率的计算必须在所有分量（包括时间分量）的中点进行。因此，在评估f2和f3时，应该是t[n]+h/2，而不是t[n+1]。