Runge-Kutta 四阶法求解二阶 ODES

Question

我想做一个简谐振子的简单例子，用龙格-库塔四阶法求解。待解的二阶常微分方程（ODE）及初始条件为：

y'' + y = 0

y(0) = 0 和 y'(0) = 1/pi

范围在0到1之间，有100步。我将二阶 ODE 分成两个一阶 ODE，使用 u 作为辅助变量：

y' = u

u' = -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) 进行了比较，但我看不到错误。在link中有两个代码，一个是Matlab的，一个是Fortran的。即使那样，我也找不到我的错误。谁能帮帮我？

Answer 1

我的代码是正确的。解析解错了。正确的解析答案是

sin(x)/pi

正如 LutzL 指出的那样。下面，可以看到解析解和数值解。限制是从 a=0 到 b=6.5.

Runge-Kutta 四阶法求解二阶 ODES

Runge-Kutta 4th order method to solve second-order ODES

python

ode

runge-kutta