点击上方「蓝字」关注我们
在自然科学的许多领域特别是科学与工程计算中，经常遇到常微分方程的求解问题。然而，只有非常少数且十分简单的微分方程可以用初等方法求得它们的解，多数情形只能利用近似方法求解，如幂级数解法、皮卡（Picard）逐步逼近法等。这些方法可以给出解的近似表达式，通常称为近似解析方法，主要是通过手算和符号计算软件如 Maple、Mathematica 实现。还有一类近似方法称为数值方法，这些方法可以给出解在一些离散点上的近似值，主要是利用计算机编写程序来处理。我们所要研究的正是这类方法。
前文中我们已经介绍了经典的欧拉方法和梯形方法，
在这里继续研究改进的欧拉方法。
问题描述
本文主要研究的对象为一阶常微分方程初值问题
其中， 为 的已知函数， 为给定的初始值，将上述问题的精确解记为 。
改进的欧拉方法
从理论上讲，梯形方法虽然提高了计算精度，但计算成本增加了，因为每迭代一次都要重新计算函数 的值，而且迭代又要反复进行若干次，计算量较大。
为了降低计算成本，在实际计算时，通常是采用迭代一次的梯形方法，即先用欧拉公式求得一个初步的近似值 ，称为预估值。预估值可能精度很差，于是就用梯形公式对它校正一次得 ，这个结果称为校正值，而这样建立的预估一校正系统称为改进的欧拉方法。具体地，改进的欧拉方法计算公式为：
上式称为改进的欧拉公式，也叫作预估校正公式，其中第一式叫预估公式，第二式叫校正公式。
或者
例题
用改进的欧拉方法求解常微分方程初值问题：
取步长为 0.1。该方程为伯努利方程，存在精确解 。
Fortran 实现
利用 Fortran 语言，我们可以高效地实现该问题的数值解法。为了便于比较数值算法的效果，我们同时给出原方程的精确解。具体方法参见以下的代码和注释：
program main
implicit none
integer :: N, i
real(kind=8) :: a, b, h, y_predict, 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_predict = y(i)+h*f(x(i),y(i)) !改进欧拉方法
y(i+1) = y(i)+h*(f(x(i),y(i))+f(x(i+1),y_predict))/2.0
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.005817，对比前文中梯形方法的 0.002141，误差是相同数量级的，计算效果相差不多。由于改进的欧拉方法不需要迭代，计算量就明显会少一些，所以，在实际应用中，改进的欧拉方法更加实用些。
Gnuplot 绘图
为了更好地展示计算结果，我们使用 Gnuplot 来读取 Fortran 程序的输出文件并绘图。在 Fortran 程序同一目录下，创建并运行以下 plt 脚本文件，就可以轻松得到 png 格式的图片文件。运行结果和对应的脚本代码如下：
set terminal png
set output "res.png"
set title "Modified Euler 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
网格加密
对于改进的欧拉方法，如果将上例中的节点数加倍，即取步长为 ，相当于对原网格 进行网格加密，重新用上述程序计算，可以得到 20 个节点 处的近似值，为了看出加密后的效果，我们对比两次网格计算结果并列表如下。
从表中可以看出，网格加密以后误差更小。这与理论分析是一致的，充分地体现了该方法的收敛性。
结语
数值求解常微分方程是一项充满挑战的任务，它不仅要求我们对数学理论有深刻理解，还需要掌握适当的编程技巧。无论是科学研究还是工业设计，数值方法都是探索未知世界的宝贵工具。希望本文能为你打开探索常微分方程数值求解世界的大门，激发你对该领域的兴趣与热情。
往期推荐
推荐阅读
FEtch 系统是笔者团队开发的新一代有限元软件开发平台。只需按照有限元语言格式填写脚本文件，即可在线自动生成基于现代 Fortran 的有限元计算程序，从而大幅提高 CAE 软件的开发效率。欢迎私信交流。
有任何疑问或建议，欢迎加Q群 "FEtch有限元开发系统(519166061)" 留言讨论。我们长期开展 FEtch 系统的试用活动，感兴趣的朋友入群后可直接联系管理员，免费获取许可证文件。
