我试图做一个谐振子的简单例子,它将由Runge-Kutta四阶方法解决。要求解的二阶常微分方程(ODE)和初始条件为:
y''+ y = 0
y(0)= 0且y'(0)= 1 / pi
范围介于0和1之间,共有100个步骤。我使用u作为辅助变量将我的二阶ODE分成两个一阶ODE:
你'=你
你'= -y
解析解是正弦y(x)=(1 / pi)^ 2 sin(pi * x)。
我的Python代码如下:
from math import pi
from numpy import arange
from matplotlib.pyplot import plot, show
# y' = u
# u' = -y
def F(y, u, x):
return -y
a = 0
b = 1.0
N =100
h = (b-a)/N
xpoints = arange(a,b,h)
ypoints = []
upoints = []
y = 0.0
u = 1./pi
for x in xpoints:
ypoints.append(y)
upoints.append(u)
m1 = h*u
k1 = h*F(y, u, x) #(x, v, t)
m2 = h*(u + 0.5*k1)
k2 = h*F(y+0.5*m1, u+0.5*k1, x+0.5*h)
m3 = h*(u + 0.5*k2)
k3 = h*F(y+0.5*m2, u+0.5*k2, x+0.5*h)
m4 = h*(u + k3)
k4 = h*F(y+m3, u+k3, x+h)
y += (m1 + 2*m2 + 2*m3 + m4)/6
u += (k1 + 2*k2 + 2*k3 + k4)/6
plot(xpoints, ypoints)
show()
所有代码都按照LutzL的建议进行了修正。见下面的评论。
代码正在运行,但我的数值解决方案与分析解决方案不匹配。我制作了一个显示以下两种解决方案的图表。我在互联网上将我的脚本与其他代码(https://math.stackexchange.com/questions/721076/help-with-using-the-runge-kutta-4th-order-method-on-a-system-of-2-first-order-od)进行了比较,但我看不出错误。在链接中,有两个代码,一个是Matlab,另一个是Fortran。即便如此,我也找不到自己的错误。谁能帮我?