これらは私が書いた関数です:
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 (4 次ルンゲ クッタ法) を使用して得られた 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 として使用できるようにするにはどうすればよいですか?