用二阶ODE在Python中实现3/8 Runge-Kutta

Question

我想求解方程y''+5y'+6y=cos(t(。由于这是二阶，我首先创建了一个一阶常微分方程系统，其中dy/dt=z=f(t，y，z(和dz/dt=cos(t(-5z-6y=g(t，y，z(。我是python的新手，不确定如何用ODE系统准确地实现程序，但对于我的输入，我写了func=[f，g]和ic=[y0，dy0]。

其次，由于Runge-Kutta计算的变量表达式的值仅在运行时可用，所以我使用python的eval函数，尽管我收到了一个错误。以下是我的代码和错误消息：

def f(t,y,z): 
return z
def g(t,y,z): 
return np.cos(t)-5*z-6*y
t = np.linspace(0,20,100)
y = np.zeros((1,len(t)))
z = np.zeros((1,len(t)))
fun = [f(t,y,z),g(t,y,z)]
ic = [1,0] # y0,dy0
def ruku(fun,h):
t0=0
tf=20

t=t0
y=ic
fc=0
while t < tf:
if t+h > tf:
exit()
if h == tf-t:
exit()
k1=eval('fun',t,y)
k2=eval('fun',t+h/3,y+h*k1/3)
k3=eval('fun',t+2*h/3,y+h*(k2-k1/3))
k4=eval('fun',t+h,y+h*(k3-k2+k1))

y=y+h*(k1+3*(k2+k3)+k4)/8
fc=fc+4
t=t+h
return y

ruku(fun,0.01)
-->TypeError: globals must be a dict

提前感谢您的任何建议；我真的在努力让它发挥作用。

Answer 1

原来我低估了Python的简单性！将我的循环替换为：

while t +h <= tf:

k1=fun(t,y)
k2=fun(t+h/3,y+h*k1/3)
k3=fun(t+2*h/3,y+h*(k2-k1/3))
k4=fun(t+h,y+h*(k3-k2+k1))

y=y+h*(k1+3*(k2+k3)+k4)/8
fc=fc+4
t=t+h

并将我的输入转换为np.arrays