这些是我写的函数:
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?
最佳答案
您的第一个 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])
关于python - python中的Runge Kutta方法,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/28509925/