求解二阶 ODES 的 Runge-Kutta 四阶方法 [英] Runge-Kutta 4th order method to solve second-order ODES
问题描述
我正在尝试做一个简单的谐振子示例,它将通过 Runge-Kutta 4 阶方法解决.要求解的二阶常微分方程 (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从 numpy 导入范围从 matplotlib.pyplot 导入图,显示#y' = u#u' = -ydef F(y, u, x):返回 -y一 = 0b = 1.0N = 100h = (b-a)/Nxpoints = arange(a,b,h)ypoints = []upoints = []y = 0.0u = 1./pi对于 x 中的 x 点:ypoints.append(y)upoints.append(u)m1 = h*uk1 = 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)/6u += (k1 + 2*k2 + 2*k3 + k4)/6绘图(xpoints,ypoints)表演()
所有代码都按照 LutzL 的建议进行了更正.请参阅下面的评论.
代码正在运行,但我的数值解与解析解不匹配.我做了一个图表,显示了下面的两个解决方案.我将我的脚本与其他一些代码进行了比较(
我的代码是正确的.解析解是错误的.正确的分析答案是
sin(x)/pi
正如 LutzL 指出的那样.下面,可以看到解析解和数值解.限制是从 a=0 到 b=6.5.
I am trying to do a simple example of the harmonic oscillator, which will be solved by Runge-Kutta 4th order method. The second-order ordinary differential equation (ODE) to be solved and the initial conditions are:
y'' + y = 0
y(0) = 0 and y'(0) = 1/pi
The range is between 0 and 1 and there are 100 steps. I separated my 2nd order ODE in two first-order ODEs, using u as auxiliary variable:
y' = u
u' = -y
The analytical solution is sinusoidal y(x) = (1/pi)^2 sin(pi*x).
My Python code is below:
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()
All code was corrected as suggested by LutzL. See comments below.
The code is running but my numerical solution does not match with the analytical solution. I made a graph showing the two solutions below. I compared my script with some other's codes (https://math.stackexchange.com/questions/721076/help-with-using-the-runge-kutta-4th-order-method-on-a-system-of-2-first-order-od) on internet and I cannot see the error. In the link, there are two codes, a Matlab one and Fortran one. Even then, I cannot find my mistake. Can anyone help me?
My code is correct. The analytical solution was wrong. The correct analytical answer is
sin(x)/pi
as LutzL pointed. Below, one can see the analytical and numerical solutions. The limits are from a=0 to b=6.5.
这篇关于求解二阶 ODES 的 Runge-Kutta 四阶方法的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!