1.引言
在数值分析中，龙格-库塔法（Runge-Kutta methods）是用于非线性常微分方程的解的重要的一类隐式或显式迭代法。这些技术由数学家卡尔·龙格和马丁·威尔海姆·库塔于1900年左右发明。
龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法，其中包括著名的欧拉法，用于数值求解微分方程。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。
在各种龙格-库塔法当中有一个方法十分常用，以至于经常被称为“RK4”或者就是“龙格－库塔法”。该方法主要是在已知方程导数和初值信息，利用计算机仿真时应用，省去求解微分方程的复杂过程。
对一般微分方程有：
{
y
′
=
f
(
x
,
y
)
y
(
x
0
)
=
y
0
\begin{cases} y'=f(x,y)\\ y(x_0)=y_0 \end{cases}
{y′=f(x,y)y(x0)=y0
在x的取值范围内将其离散为
n
n
n段，定义步长，令第
n
n
n步对应的函数值为
y
n
y_n
yn。于是通过一系列的推导可以得到下一步的
y
n
+
1
y_{n+1}
yn+1值为
y
n
+
1
=
y
n
+
h
6
(
K
1
+
2
K
2
+
2
K
3
+
K
4
)
y_{n+1}=y_n+\frac{h}{6} (K_1+2K_2+2K_3+K_4)
yn+1=yn+6h(K1+2K2+2K3+K4)
其中
{
K
1
=
f
(
x
n
,
y
n
)
K
2
=
f
(
x
n
+
h
2
,
y
n
+
h
2
K
1
)
K
3
=
f
(
x
n
+
h
2
,
y
n
+
h
2
K
2
)
K
4
=
f
(
x
n
+
h
,
y
n
+
h
K
3
)
\begin{cases} K_1=f(x_n, y_n) \\ K_2=f(x_n+\frac{h}{2}, y_n+\frac{h}{2}K_1) \\ K_3=f(x_n+\frac{h}{2},y_n+\frac{h}{2}K_2) \\ K_4=f(x_n+h,y_n+hK_3) \end{cases}
⎩
⎨
⎧K1=f(xn,yn)K2=f(xn+2h,yn+2hK1)K3=f(xn+2h,yn+2hK2)K4=f(xn+h,yn+hK3)
2.Python代码实现
以如下微分方程为例：
{
y
′
=
y
−
2
x
y
y
(
0
)
=
1
\begin{cases} y'=y - \frac{2x}{y} \\ y(0)=1 \end{cases}
{y′=y−y2xy(0)=1
其中，
0
≤
x
≤
2
0\le x \le 2
0≤x≤2。
已知，上述微分方程的函数解析表达式为：
y
=
2
x
+
1
y=\sqrt{2x+1}
y=2x+1
##################################
import matplotlib.pyplot as plt
import math
def f_function(x,y): #完成对微分方程的表达
y_pie=0.0
#y_pie=math.cos(x) #对于不同的微分方程，这里采用不同的表达式
y_pie=y-2.0*x/y
return y_pie
def f_accurate(x):
y_accu=math.sqrt(2.0*x+1)#计算精确解
return y_accu
#存储各步进状态时的值
xlist=[]
ylist=[]
#存储解析解的值
y_accurate_list=[]
#赋初值
x=0.0
y=1.0
h=0.1
#对x范围内步进
while x <=2.0:
#将步进结果放进列表中
xlist.append(x)
ylist.append(y)
#四步龙格-库塔核心步骤
k1=f_function(x,y)
k2=f_function(x+h/2,y+h/2*k1)
k3=f_function(x+h/2,y+h/2*k2)
k4=f_function(x+h,y+h*k3)
y= y+h/6*(k1+2*k2+2*k3+k4) #得到yn+1的值
x=x+h #x步进
#输出各步进状态的计算结果
for i in range(len(xlist)):
print(xlist[i],ylist[i])
y_accurate_list.append(f_accurate(xlist[i]))#求精确解
#将结果绘图
plt.plot(xlist,ylist,'x-',color='red',label='R-K Results')
plt.plot(xlist,y_accurate_list,'o',color='green',label='Accurate Results')
plt.legend()
#plt.title('Computation Comare')
plt.xlabel("x")
plt.ylabel("y")
plt.show()
print('over')
3.计算过程如下：
其中第一列为x，第二列为龙格-库塔法计算结果，第三列为精确解。
0.0 1.00000 1.00000
0.1 1.09545 1.09545
0.2 1.18322 1.18322
0.3 1.26491 1.26491
0.4 1.34164 1.34164
0.5 1.41422 1.41421
0.6 1.48324 1.48324
0.7 1.54920 1.54919
0.8 1.61246 1.61245
0.9 1.67332 1.67332
1.0 1.73206 1.73205
1.1 1.78886 1.78885
1.2 1.84392 1.84391
1.3 1.89738 1.89737
1.4 1.94937 1.94936
1.5 2.00001 2.00000
1.6 2.04941 2.04939
1.7 2.09764 2.09762
1.8 2.14478 2.14476
1.9 2.19092 2.19089
计算结束
结合下图可以看出，四阶龙格-库塔法的计算精度是比较高的。
参考文献
https://www.bilibili.com/read/cv15448483/
