点击上方「蓝字」关注我们
在自然科学的许多领域特别是科学与工程计算中，经常遇到常微分方程的求解问题。然而，只有非常少数且十分简单的微分方程可以用初等方法求得它们的解，多数情形只能利用近似方法求解，如幂级数解法、皮卡（Picard）逐步逼近法等。这些方法可以给出解的近似表达式，通常称为近似解析方法，主要是通过手算和符号计算软件如 Maple、Mathematica 实现。还有一类近似方法称为数值方法，这些方法可以给出解在一些离散点上的近似值，主要是利用计算机编写程序来处理。我们所要研究的正是这类方法。
前文中我们已经介绍了经典的欧拉方法，
在这里继续研究梯形方法。
问题描述
本文主要研究的对象为一阶常微分方程初值问题
其中， 为 的已知函数， 为给定的初始值，将上述问题的精确解记为 。
梯形方法
如果从另外的角度看欧拉方法，就可以引申出其他一些新的方法。事实上，除了用泰勒公式将导数用差商来近似外，我们可以直接对式 (1) 中的方程两边从 到 积分，得
上式右端的积分项中被积函数既含有 ，又含有未知函数 ，无法求出积分值，只能采用数值的方法进行计算。于是对上式右端的积分项用数值积分中的左矩形公式逼近，即
并用 作为 的近似值，这样式 (2) 近似为差分格式
即欧拉方法公式，故欧拉方法也称为矩形法。
从上面的分析也可以看出本质上欧拉方法效果不是很好的主要原因是它采用了精度较差的矩形法计算右端积分。此外，如果用右矩形公式计算式（2）右端的积分，则类似可得差分格式
称为后退欧拉方法。
可以想象，如果改用梯形公式计算式（2）右端的积分，可期望得到较高的精度。这时
将这个结果代入式（2）并将其中的 用 近似代替，则得新的差分格式：
这个方法称为梯形方法。
上述 3 个计算格式，即式 (3)、式 (4) 及式（5）都是由 计算 ，这种只用前一步 即可算出 的方法称为单步法。其中，式（3）可由 逐次求出 的值，称为显式方法，而式 (4) 及式 (5) 左右两端都含有 ，无法直接从 求出 ，而需要解含有 的方程，求解很不容易，常用迭代法求解，这类方法称为隐式方法。
迭代公式
注意到梯形方法是隐式方法，在实际求解时需要用迭代法。而迭代的初值则由欧拉方法提供，即有以下梯形方法的迭代公式:
使用迭代法时，需要设立迭代误差限 以控制迭代步数，也就是说，先用欧拉方法由 得出 的初始近似值 ，然后用式 (6) 中的第二个公式进行迭代，反复迭代直到 为止，并把 取作 的近似值 。
显然，应用梯形方法，如果序列 收敛，它的极限 便满足方程
比较式 (5) ，可知序列的极限可取作 。
例题
继续前文的例题，用梯度方法求解常微分方程初值问题：
该方程为伯努利方程，存在精确解 。利用式 (6) 编程，这里取步长为 0.1 ，迭代误差限为 。
Fortran 实现
利用 Fortran 语言，我们可以高效地实现该问题的数值解法。为了便于比较数值算法的效果，我们同时给出原方程的精确解。具体方法参见以下的代码和注释：
program main
implicit none
integer :: N, i, k
real(kind=8) :: a, b, h, ye, err
real(kind=8) ytemp1, ytemp2, epsilon
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
epsilon = 1e-4 !迭代误差阈值
!打开存储文件
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
ytemp1 = y(i)+h*f(x(i),y(i)) !梯形法
do k = 1, 100 !迭代求解
ytemp2 = y(i)+h*(f(x(i),y(i))+f(x(i+1),ytemp1))/2.0
if (abs(ytemp1-ytemp2)<epsilon) exit
ytemp1=ytemp2
enddo
y(i+1)=ytemp2
ye = yexact(x(i+1)) !计算精确解
err = abs(y(i+1) - ye) !计算节点处的误差
!保存结果
write(10,"(4e12.5,i3)") x(i+1), y(i+1), ye, err, k
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.002141, 对比前文中欧拉方法的 0.052720，梯形方法的结果明显好于后者。
Gnuplot 绘图
为了更好地展示计算结果，我们使用 Gnuplot 来读取 Fortran 程序的输出文件并绘图。在 Fortran 程序同一目录下，创建并运行以下 plt 脚本文件，就可以轻松得到 png 格式的图片文件。运行结果和对应的脚本代码如下：
set terminal png
set output "res.png"
set title "Trapezoidal 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
结语
数值求解常微分方程是一项充满挑战的任务，它不仅要求我们对数学理论有深刻理解，还需要掌握适当的编程技巧。无论是科学研究还是工业设计，数值方法都是探索未知世界的宝贵工具。希望本文能为你打开探索常微分方程数值求解世界的大门，激发你对该领域的兴趣与热情。
往期推荐
推荐阅读
FEtch 系统是笔者团队开发的新一代有限元软件开发平台。只需按照有限元语言格式填写脚本文件，即可在线自动生成基于现代 Fortran 的有限元计算程序，从而大幅提高 CAE 软件的开发效率。欢迎私信交流。
有任何疑问或建议，欢迎加Q群 "FEtch有限元开发系统(519166061)" 留言讨论。我们长期开展 FEtch 系统的试用活动，感兴趣的朋友入群后可直接联系管理员，免费获取许可证文件。
