3

这些是我编写的函数：

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（四阶龙格库塔法）获得的 fLorenz（洛伦兹方程组）的数值解。我在使用 fLorenz 函数时遇到问题。当我如上所述定义它并调用 lorenzPlot 时，我收到一条错误消息

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

我猜这与无法正确乘以数组有关。
但是，当我将 fLorenz 更改为

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 以便它可以用作 rk4 和 lorenzPlot 中的 f？

4

3 回答 3

2

您的第一个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])
于 2015-02-13T23:19:34.350 回答
1

你fLorenz没有返回任何东西，即，由于有据可查的 Python 约定，它返回None.

你的错误信息告诉你

* 不支持的操作数类型：“float”和“NoneType”

并且错误发生在语句中

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

在上一行中，您分配f1了调用返回的内容fLorenz，即None.

如果您更改函数fLorenz以使其返回数值，无论是标量还是数组，也许您可​​以继续进行计算。

于 2015-02-13T23:10:30.617 回答
1

除了实现错误，你对RK4方法的理解是不完整的。中点斜率的评估必须发生在所有分量的中点，包括时间分量。因此，它应该t[n]+h/2代替 和的t[n+1]评估。f2f3

于 2015-02-26T21:37:19.727 回答