Я пытаюсь сделать простой пример гармонического осциллятора, который будет решен методом Рунге-Кутты 4-го порядка. Обыкновенное дифференциальное уравнение второго порядка (ОДУ), подлежащее решению, и начальные условия:
y '' + y = 0
y (0) = 0 и y '(0) = 1 / pi
Диапазон составляет от 0 до 1 и 100 шагов. Я разделил свой ODE 2-го порядка на два 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) в Интернете, и я не вижу ошибки. В ссылке есть два кода, один Matlab и Fortran. Даже тогда я не могу найти свою ошибку. Может кто-нибудь мне помочь?
