私はPythonの初心者/中級者です。4 次のルンゲクッタ法 (RK4)を Python にコーディングしました。基本的には振り子を解いていますが、それはここでは重要ではありません。
RK4 メソッドを次のように改善したいと考えています。関数 f を RK4 関数に直接渡すことができるようにしたいです。つまり、RK4(y_0, n, h) は RK4(f,y_0,n,h) になる必要があります。これには、この 1 つの振り子だけでなく、他のシステムを記述する他の f 関数に RK4 を使用できるという大きな利点があります。
単純な関数を RK4 に渡すだけで遊んでみましたが、何か間違ったことをしています。Pythonでこれを行うにはどうすればよいですか?
import numpy as np
def RK4(y_0, n, h):
#4th order Runge-Kutta solver, takes as input
#initial value y_0, the number of steps n and stepsize h
#returns solution vector y and time vector t
#right now function f is defined below
t = np.linspace(0,n*h,n,endpoint = False) #create time vector t
y = np.zeros((n,len(y_0))) #create solution vector y
y[0] = y_0 #assign initial value to first position in y
for i in range(0,n-1):
#compute Runge-Kutta weights k_1 till k_4
k_1 = f(t[i],y[i])
k_2 = f(t[i] + 0.5*h, y[i] + 0.5*h*k_1)
k_3 = f(t[i] + 0.5*h, y[i] + 0.5*h*k_2)
k_4 = f(t[i] + 0.5*h, y[i] + h*k_3)
#compute next y
y[i+1] = y[i] + h / 6. * (k_1 + 2.*k_2 + 2.*k_3 + k_4)
return t,y
def f(t,vec):
theta=vec[0]
omega = vec[1]
omegaDot = -np.sin(theta) - omega + np.cos(t)
result = np.array([omega,omegaDot])
return result
test = np.array([0,0.5])
t,y = RK4(test,10,0.1)