点击上方「蓝字」关注我们
在自然科学的许多领域特别是科学与工程计算中，经常遇到常微分方程的求解问题。然而，只有非常少数且十分简单的微分方程可以用初等方法求得它们的解，多数情形只能利用近似方法求解，如幂级数解法、皮卡（Picard）逐步逼近法等。这些方法可以给出解的近似表达式，通常称为近似解析方法，主要是通过手算和符号计算软件如 Maple、Mathematica 实现。还有一类近似方法称为数值方法，这些方法可以给出解在一些离散点上的近似值，主要是利用计算机编写程序来处理。我们所要研究的正是这类方法。
前文中我们已经介绍了经典的欧拉方法、梯形方法和改进的欧拉方法，
在这里继续研究龙格-库塔法（Runge-Kutta methods）。这类方法能够提供比欧拉方法更高的精度，并且相对简单易用。
问题描述
本文主要研究的对象为一阶常微分方程初值问题
其中， 为 的已知函数， 为给定的初始值，将上述问题的精确解记为 。
龙格-库塔法
龙格-库塔法（Runge-Kutta methods）具有较高的精度和稳定性，适用于大多数非刚性常微分方程问题的求解。
它的基本思想是通过计算不同点上的函数值，并对这些函数值作线性组合来构造近似公式，再把近似公式与泰勒级数方法进行比较，使前面的若干项相同，从而使近似公式达到一定的阶数。
常用的最经典的是四阶龙格-库塔法。这种方法的局部截断误差为 ，整体误差为 。
对于每个网格节点 ，计算 的四阶龙格-库塔公式为：
其中
给出了从 到 的初始斜率。 和 都是在 处计算的，但使用了不同的 值，这使得方法能够更好地捕捉到 的变化趋势。 使用了 处的信息，以及 提供的 值。 最终的 是基于 ，， 和 的加权平均，这种加权方式确保了该方法的高阶精度。
以上公式的推导过程较为复杂，这里不再给出，感兴趣的读者可以下载并阅读以下 pdf 文档：
https://pan.baidu.com/s/1og3wQTTDfo_fuQUuPEhMLw?pwd=d6js
例题
用四阶龙格-库塔法求解常微分方程初值问题：
取步长为 0.2。该方程为伯努利方程，存在精确解 。
Fortran 实现
利用 Fortran 语言，我们可以高效地实现该问题的数值解法。为了便于比较数值算法的效果，我们同时给出原方程的精确解。具体方法参见以下的代码和注释：
program main
implicit none
integer :: N, i
real(kind=8) :: a, b, h, y_predict, ye, err
real(kind=8) :: k1, k2, k3, k4
real(kind=8), allocatable, dimension(:) :: x, y
a = 0.0D0
b = 1.0D0
N = 5 !总的剖分数
h = (b - a) / dble(N) !步长
!动态分配长度为(N+1)的数组
!x存放节点坐标
!y存放对应节点的数值解
allocate(x(N+1), y(N+1))
!节点坐标
x = [(a + i * h, i = 0, N)]
y(1) = 1.0D0 !初值
err = 0.0D0
!打开存储文件
open(10, file="res.txt")
write(10,*) "# x y yexact err"
!保存初值
write(10,"(4e12.5)") x(1), y(1), yexact(x(1)), err
do i = 1, N !4阶龙格-库塔法
x(i+1) = x(i) + h
k1 = h*f(x(i), y(i))
k2 = h*f(x(i) + 0.5D0*h, y(i) + 0.5D0*k1)
k3 = h*f(x(i) + 0.5D0*h, y(i) + 0.5D0*k2)
k4 = h*f(x(i+1), y(i) + k3)
y(i+1) = y(i) + (k1 + 2D0 * k2 + 2D0 * k3 + k4) / 6.0D0
ye = yexact(x(i+1)) !计算精确解
err = abs(y(i+1) - ye) !计算节点处的误差
!保存结果
write(10,"(4e12.5)") x(i+1), y(i+1), ye, err
end do
close(10)
contains
!右端项函数
function f(x, y) result (res)
real(kind=8), intent(in) :: x, y
real(kind=8) :: res
res = y - 2.0D0 * x / y
end function f
!精确解y
function yexact(x) result (res)
real(kind=8), intent(in) :: x
real(kind=8) :: res
res = sqrt(1.0D0 + 2.0D0 * x)
end function yexact
end program main
编译并运行，将计算结果总结如下。
Gnuplot 绘图
为了更好地展示计算结果，我们使用 Gnuplot 来读取 Fortran 程序的输出文件并绘图。在 Fortran 程序同一目录下，创建并运行以下 plt 脚本文件，就可以轻松得到 png 格式的图片文件。运行结果和对应的脚本代码如下：
set terminal png
set output "res.png"
set title "Runge-Kutta Method"
set xlabel "x"
set ylabel "y"
set key center left spacing 2
plot "res.txt" using 1:2 title "Numerical", \
"res.txt" using 1:3 title "Analytical" \
with lines linewidth 2
小结
本文介绍了四阶龙格-库塔法的基本原理，并通过一个具体的例子展示了如何使用 Fortran 语言来实现这一方法。这种方法不仅适用于一阶常微分方程，也可以通过适当的转换应用到更高阶的微分方程上。通过这种方式，我们可以解决许多现实世界中（例如物理、工程和经济学等领域）的数学建模问题。
往期推荐
推荐阅读
FEtch 系统是笔者团队开发的新一代有限元软件开发平台。只需按照有限元语言格式填写脚本文件，即可在线自动生成基于现代 Fortran 的有限元计算程序，从而大幅提高 CAE 软件的开发效率。欢迎私信交流。
有任何疑问或建议，欢迎加Q群 "FEtch有限元开发系统(519166061)" 留言讨论。我们长期开展 FEtch 系统的试用活动，感兴趣的朋友入群后可直接联系管理员，免费获取许可证文件。
