点击上方「蓝字」关注我们
在自然科学的许多领域特别是科学与工程计算中，经常遇到常微分方程的求解问题。然而，只有非常少数且十分简单的微分方程可以用初等方法求得它们的解，多数情形只能利用近似方法求解，如幂级数解法、皮卡（Picard）逐步逼近法等。这些方法可以给出解的近似表达式，通常称为近似解析方法，主要是通过手算和符号计算软件如 Maple、Mathematica 实现。还有一类近似方法称为数值方法，这些方法可以给出解在一些离散点上的近似值，主要是利用计算机编写程序来处理。我们所要研究的正是这类方法。
问题描述
本文主要研究的对象为一阶常微分方程初值问题
其中， 为 的已知函数， 为给定的初始值，将上述问题的精确解记为 。
正如上面所提到的，即使形如式（1）这样形式简单的方程也未必能找到精确解。事实上，C. Henry Edwards 和 David E. Penney 在其所著的《常微分方程基础》一书中这样描述：
对于具有以上一般形式的微分方程，如果只用初等的方法就可以精确地求解出来，那么这种情况也不过是一个例外，并非普遍规则。
It is the exception rather than a rule when a differential equation of the above general form can be solved exactly and explicitly by elementary methods.
作为我们要研究的对象，原则上只有当初值问题式（1）的解存在且唯一时，使用数值解法去求解才有意义，这一前提条件由下面的存在唯一性定理保证。
存在唯一性定理
设函数 在区域 ， 上连续，且在区域 内关于 满足李普希兹（Lipschitz）条件，即存在正数 ，使得对于 内任意两点 与 ，恒有
则初值问题式（1）的解 存在且唯一。
（此定理的证明可在普通的《常微分方程》教材上找到，此处从略。）
当然，有时候我们也会遇到一些实际问题，在理论上难以证明其解的存在唯一性，或者根本不知道原问题解，这些问题属于微分方程基本理论研究的范畴，这时通常可以采用数值方法先进行试探性的数值模拟，然后根据数值方法的结果，来提供一些理论上的证明方向与思路。
数值求解思路
用数值方法求解问题式（1）的基本思想是：
在求解区域 上，任取 个节点，在这些节点上采用离散化方法（通常用数值积分、微分、泰勒展开式等），将上述初值问题化成关于离散变量的相应问题，并且将所得离散问题的解 作为精确解 的近似值，这样求得的 就是上述初值问题在节点 上的数值解。
一般说来，不同的离散过程将得到不同的数值方法，从而得到不同的解。
数值方法的具体操作过程如下：
网格剖分（在求解区域上作剖分）。
在求解区域 上，取 个节点
这里 称为 到 的步长， 。这些 可以不相等，但有时为了计算方便，特别是编写程序方便，一般取成相等，此时就有等步长 ，这样的区域剖分称为等距剖分或者一致网格剖分， 称为网格节点。整个剖分过程相当于对一维的求解区域进行离散化。对连续方程作节点离散，使连续方程仅在离散节点上成立，从而将原方程弱化，得到节点离散方程。
在节点上采用离散化方法（忽略高阶项），并用 作为精确解 的近似值，可得差分格式或差分方程。
将差分方程改写为迭代格式，或者线性方程组的形式，并求解，得出精确解 的近似值 。
从理论上研究数值格式的局部截断误差（即相容性）、稳定性以及数值解的收敛性与整体误差。
分析数值结果与理论分析是否一致，考察有无局限性及可改进的方法。
注意，有时可以将第 4 步和第 5 步交换次序，先在理论上分析第 3 步得到的数值格式的数值效果。如果效果不好的话就没必要进行第 4 步的实际操作了。当然，以上操作过程比较抽象，我们需要通过具体的方法来实现这些过程，以加深大家对上述操作步骤的理解 。
欧拉（Euler）方法
首先，对问题式(1)的求解区域即区间 作等距剖分，得到等距节点 ， ，步长为 。
其次，将原连续方程 即 弱化，使之在散节点处成立，即 。
接着，采用差商近似微商的离散化方法，在离散节点 处应有节点离散方程：
忽略高阶项，并用数值解 作为精确解 的近似值，可得差分方程：
这样，原来的连续方程初值问题式（1）可用以下差分格式作数值逼近
这就是著名的欧拉方法，也叫欧拉折线法。
求解初值问题式（1）其实就是在 平面上求一条通过点 的曲线（也称为积分曲线） ，并使曲线上任意一点 处的切线斜率为 。
事实上，欧拉方法有明显的几何意义，就是从点 出发作一条斜率为 的直线，此直线交直线 于点 ， 点的纵坐标 就是 的近似值；再从点 作一条斜率为 的直线，交直线 于点 ， 点的纵坐标 就是 的近似值；按同样的过程继续下去，可得到一条折线 。该折线就是数值解的图形，它是精确解 的近似折线，这也是欧拉折线法名称的由来。
例题
用欧拉方法求解常微分方程初值问题：
该方程为伯努利方程，存在精确解 。
这里取步长为 ，相当于对 等距剖分 10 份，得到节点坐标
利用欧拉方法，得到差分方程
Fortran 实现
利用 Fortran 语言，我们可以高效地实现该问题的数值解法。为了便于比较数值算法的效果，我们同时给出原方程的精确解。具体方法参见以下的代码和注释：
program main
implicit none
integer :: N, i
real(kind=8) :: a, b, h, ye, err
real(kind=8), allocatable, dimension(:) :: x, y
a = 0.0D0
b = 1.0D0
N = 10 !总的剖分数
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.0
!打开存储文件
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
y(i+1) = y(i) + h * f(x(i), y(i)) !欧拉方法
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
编译并运行，将计算结果总结如下。
从表中可见，用欧拉方法得到的数值结果最大误差为 0.052720, 效果并不是很好。为此，我们可以通过减小步长来进一步提高精度，大家可自行尝试。
Gnuplot 绘图
为了更好地展示计算结果，我们使用 Gnuplot 来读取 Fortran 程序的输出文件并绘图。在 Fortran 程序同一目录下，创建并运行以下plt脚本文件，就可以轻松得到 png 格式的图片文件。运行结果和对应的脚本代码如下：
set terminal png
set output "res.png"
set title "Euler Method"
set xlabel "x"
set ylabel "y"
set key center left spacing 2
plot "res.txt" using 1:2 title "Numerical" \
with lines linewidth 2, \
"res.txt" using 1:3 title "Analytical"
结语
数值求解常微分方程是一项充满挑战的任务，它不仅要求我们对数学理论有深刻理解，还需要掌握适当的编程技巧。无论是科学研究还是工业设计，数值方法都是探索未知世界的宝贵工具。希望本文能为你打开探索常微分方程数值求解世界的大门，激发你对该领域的兴趣与热情。
往期推荐
推荐阅读
FEtch 系统是笔者团队开发的新一代有限元软件开发平台。只需按照有限元语言格式填写脚本文件，即可在线自动生成基于现代 Fortran 的有限元计算程序，从而大幅提高 CAE 软件的开发效率。欢迎私信交流。
有任何疑问或建议，欢迎加Q群 "FEtch有限元开发系统(519166061)" 留言讨论。我们长期开展 FEtch 系统的试用活动，感兴趣的朋友入群后可直接联系管理员，免费获取许可证文件。
