基于欧拉方法求解微分方程初值问题

Solving Initial Value Problems of Differential Equations Based on Euler Method

摘要: 微分方程在数学实际应用以及物理、化学、生物等众多领域中被应用得十分广泛。欧拉法是求解微分方程最基本、最常用的方法，而改进欧拉法则是在欧拉法的基础上进行相应的改进而形成的方法。本文对比关于欧拉法的三种方法近似解与精确解的关系，并根据其具体计算结果进行对比，得到改进欧拉法的计算效率更高、精度更高。

Abstract: Differential equations are widely used in many fields such as mathematics, physics, chemistry, bi-ology and so on. Euler method is the most basic and commonly used method to solve differential equation, and the improved Euler method is a method based on the corresponding improvement of Euler method. In this paper, the relationship between the approximate solution and the exact solution of the three Euler methods is compared, and according to the specific calculation results, the improved Euler method is more efficient and more accurate.

1. 引言

目前，在科学工程技术及自然科学领域中，很多问题都可以用常微分方程的初值问题来解决。微分方程逐渐随着微积分学一起发展起来，不仅在解决与导数有关的问题中发挥重要的作用，而且在物理学中许多涉及变力的运动学 [1]、动力学 [2] 问题，如空气的阻力为速度函数的落体运动等问题，以及航空航天、工程学、天文学等领域建立模型 [3] [4] [5] [6] 中都有十分广泛的应用。但是，只有少数简单的微分方程可以求得精确解。在无法求得精确解时，可以通过求出近似解的方式代替这些精确解，以此达到解决这些实际问题的目的，而这种求出近似解的方法，便是微分方程的初值解法，有一定的有效性且有一定的准确度 [7]。

近年来，微分方程初值问题的基本数值方法被应用得越来越广泛，包括欧拉法，后退欧拉法，梯形方法以及改进的欧拉方法等。其中欧拉法是求解常微分方程中最基本、最常用的方法，而改进欧拉法则是在欧拉法的基础上进行相应的改进而形成的方法。本文针对具体的数值算例对基于欧拉法三种基本的数值方法所求出的近似解与精确解的关系进行分析。得到改进的欧拉法相对于其他两种数值方法，不仅计算精度高，而且更加稳定，可以提高计算效率。

2. 数值方法

本文考虑具有如下初值问题的微分方程：

$\{\begin{array}{l}\frac{\text{d}y}{\text{d}x}=f\left(x,y\right),x\in \left[x_0,b\right),\\ y\left(x_0\right)=y_0.\end{array}$ (1)

首先对(1)式进行离散化，定义离散节点为 $x_n=x_0+nh$， $n=0,1,2,\cdots$，其中 $h=x_{n+1}-x_n$ 为等距步长。通过数值方法求得微分方程的精确解，记为 $y\left(x_i\right)$， $i=1,2,\cdots ,n,n+1,\cdots$。在离散节点 $x_1,x_2,\cdots ,x_n,x_{n+1},\cdots$ 上的近似值记为 $y_i$， $i=1,2,\cdots ,n,n+1,\cdots$。对于常微分方程初值求解问题来说，无论是常规欧拉法还是改进后的欧拉法，都是通过数值积分方法来实现离散化处理的。

2.1. 欧拉法

对常微分方程(1)式中的导数用均差近似，即

$\frac{y\left(x_{n+1}\right)-y\left(x_n\right)}{h}\approx y^{\prime }\left(x_n\right)=f\left(x_n,y\left(x_n\right)\right),n=0,1,2,\cdots .$ (2)

进而得到形如下式的欧拉法公式，即

$\{\begin{array}{l}{y}_n{}_{+1}=y_n+hf\left(x_n,y_n\right),n=0,1,2,\cdots ,\\ y\left(x_0\right)=y_0.\end{array}$ (3)

若对微分方程(1)式从 $x_n$ 到 $x_{n+1}$ 积分，则有

$y\left(x_{n+1}\right)=y\left(x_n\right)+{\displaystyle {\int }_{x_n}^{x_{n+1}}f\left(t,y\left(t\right)\right)\text{d}t},$ (4)

其中右端积分以左矩形公式 $hf\left(x_n,y\left(x_n\right)\right)$ 近似，左端第一项以 $y_n$ 代替，并以 $y_{n+1}$ 代替 $y\left(x_{n+1}\right)$ 也可得到欧拉法公式，即(3)式。

由于其初值 $y_0$ 已知，则由(3)式可以逐次循环计算得到

$\begin{array}{l}{y}_1=y_0+hf\left(x_0,y_0\right),\\ y_2=y_1+hf\left(x_1,y_1\right),\\ ⋮\\ y_n=y_{n-1}+hf\left(x_{n-1},y_{n-1}\right),\\ ⋮\end{array}$ (5)

下面利用泰勒展开将 $y\left(x_{n+1}\right)$ 在 $x_n$ 处展开，则有

$y\left(x_{n+1}\right)=y\left(x_n+h\right)=y\left(x_n\right)+y^{\prime }\left(x_n\right)h+\frac{h^2}{2}{y}^{″}\left({\zeta }_n\right),\zeta \in \left(x_n,x_{n+1}\right)$ (6)

由于在 $y_n=y\left(x_n\right)$ 的前提下，有 $f\left(x_n,y_n\right)=f\left(x_n,y\left(x_n\right)\right)=y^{\prime }\left(x_n\right)$。进而可以得到欧拉法的局部截断误差为

$y\left(x_{n+1}\right)-y_{n+1}=\frac{h^2}{2}{y}^{″}\left({\zeta }_n\right)\approx \frac{h^2}{2}{y}^{″}\left(x_n\right).$ (7)

欧拉法的局部截断误差也可以表示为 $O\left(h^{\wedge }2\right)$。

2.2. 后退欧拉法

下面将(4)式中的右端积分以右矩形公式 $hf\left(x_{n+1},y\left(x_{n+1}\right)\right)$ 近似，则可得到

$y_{n+1}=y_n+hf\left(x_{n+1},y_{n+1}\right).$ (8)

上式即为后退的欧拉法；由于欧拉法是关于 $y_{n+1}$ 的一个直接的计算公式，被称为是显式的；然而(8)式右端中含有未知项 $y_{n+1}$，是关于 $y_{n+1}$ 的一个函数方程，这类公式被称为是隐式的。因此，后退的欧拉法也被称为隐式欧拉法。其精度和欧拉法的精度基本相同，通常采用迭代法求解，而迭代过程的实质是逐步显式化 [8]。

2.3. 改进欧拉法

下面对对欧拉法进行进一步相应的改进，以提高欧拉法的求解精度。首先，利用欧拉法初步求出微分方程初值问题的一个近似值 ${\bar{y}}_{n+1}$，记该值为预测值。结合梯形公式 $y_{n+1}=y_n+\frac{h}{2}\left[f\left(x_n,y_n\right)+f\left(x_{n+1},y_{n+1}\right)\right]$ 将其校正，进而建立如下的预测–校正系统：

$\{\begin{array}{l}预测{\bar{y}}_{n+1}=y_n+hf\left(x_n,y_n\right),\\ 校正y_{n+1}=y_n+\frac{h}{2}\left[f\left(x_n,y_n\right)+f\left(x_{n+1},{\bar{y}}_{n+1}\right)\right].\end{array}$ (9)

上式中的预测–校正系统可表示为下列平均化形式，即为

$\{\begin{array}{l}{y}_p=y_n+hf\left(x_n,y_n\right),\underset{{}_{}}{}\\ y_c=y_n+hf\left(x_{n+1},y_p\right),\\ y_{n+1}=\frac{1}{2}\left(y_p+y_c\right).\end{array}$ (10)

(10)式被称为改进的欧拉公式，其局部截断误差也可以表示为 $O\left(h^{\wedge }3\right)$。

3. 数值算例

设其计算区间为 $\left(0,1\right)$，计算步长为 $h=0.1$，应用本文总结的欧拉法，后退欧拉法，改进欧拉法分别求解具有如下初值问题的常微分方程：

$\{\begin{array}{l}\frac{\text{d}y}{\text{d}x}=y-\frac{2x}{y},x\in \left(0,1\right),\\ y\left(0\right)=1.\end{array}$

得到的近似解与精确解的数值结果关系即如图1所示。

Figure 1. The relationship between approximate solution and exact solution obtained by different methods (a) Euler method; (b) Backward Euler method; (c) Improved Euler method

图1. 不同方法求出的近似解与精确解的关系 (a) 欧拉法；(b) 后向欧拉法；(c) 改进欧拉法

并且，将通过欧拉法与改进欧拉法计算所得到近似解的计算结果与微分方程的精确解进行对比，得到的结果如表1所示。

Table 1. Comparison between Euler method and improved Euler method

表1. 欧拉法与改进欧拉法计算结果对比

由上述结果可以看出，利用改进欧拉法求解微分方程初值问题所得到的近似解与精确解之间的误差较小，能够更加接近精确解。并且精度高，计算效率更高，使得微分方程的初值求解问题更加高效，减小了数值计算的工作强度，因此改进欧拉法更加适合计算机的计算。

4. 总结

本文针对具体的数值算例对欧拉法、后退欧拉法、改进欧拉法三种基本数值方法所求出的近似解与精确解进行对比分析。得到改进欧拉法相对于其他两种数值方法，不仅计算精度最高，而且可以提高计算效率，更加适合微分方程初值问题的求解，能够被广泛应用于各个自然科学领域。