高阶线性常微分方程显式Euler法的数值分析

Numerical Analysis of High-Order Linear Ordinary Differential Equations via Explicit Euler Method

摘要: 针对高阶线性常微分方程，建立了显式Euler格式，分析了数值格式的局部截断误差，并在舍入误差条件下，获得了数值解的全局误差，同时，讨论了数值解的渐进稳定性。常系数与变系数二阶算例验证了理论分析结果的正确性与有效性。

Abstract: An explicit Euler scheme is established for high-order linear ordinary differential equations. The local truncation error of the numerical scheme is analyzed. Under the condition of rounding errors, the global error of the numerical solution is obtained. Also, the asymptotic stability of the numerical solution is discussed. The theoretical analysis results are verified to be correct and effective through constant-coefficient and variable-coefficient second-order ordinary differential equations.

1. 引言

在科学与工程实践中，许多物理问题的求解最终都可归结为常微分方程的求解。然而，绝大多数常微分方程很难获得解析解，部分方程甚至不存在解析表达式。因此，数值解法便成为求解常微分方程初边值问题的重要途径。

目前国内外对数值求解常微分方程已经有了系统性的研究和分析。文献[1]将高阶线性常微分方程转化为一阶线性常微分方程组，构建了含有一阶导数形式的最小二乘支持向量机(LS-SVM)回归模型。文献[2]提出了有效改进的中心支持向量机(Proximal Support Vector Machines, P-SVM)方法。文献[3]针对一阶同型微分方程，利用欧拉公式和改进的欧拉公式探讨了Lipschitz常数及不同的步长对整体误差的影响。文献[4]提出一种求解微分方程的神经网络方法，即通过物理约束耦合神经网络的方法。文献[5]给出三级Runge-Kutta (RK)方法的构造过程和相关细节。文献[6]在常微分方程组形式下使用PM算法，提出了分离原理将微分方程组对应的多目标优化问题转化为多个相对独立的单目标优化问题。文献[7]对线性常微分方程的初值问题的全局误差进行估计和优化求解。文献[8]从Newton-Cotes积分的角度构造了一个带参数的三级二阶显式Runge-Kutta格式和两个带参数的四级三阶显式Runge-Kutta格式。基于分数阶微积分基本定理和三次B样条理论，文献[9]构造了求解线性Caputo-Fabrizio型分数阶微分方程数值解的三次B样条方法针对一类带Caputo导数的变分数阶随机微分方程，[10]构造了Euler-Maruyama方法。文献[11]设计了一种基于文本的常微分方程求解框架，该框架结合了即时编译和耦合的CPU。文献[12]提出了新的加速梯度算法，并构建了Layapunov函数，证明了优化算法的收敛速率。

本文首先讨论了一阶常微分方程组解的存在唯一性条件；然后通过变量代换将高阶线性常微分方程转化为一阶线性常微分方程组，建立了一阶线性常微分方程的显式Euler格式，理论分析了其局部截断误差，并在舍入误差情况下分析了数值解的全局误差，以及显式Euler法的渐进稳定性；最后通过数值实验，验证了高阶线性常微分方程显式Euler法的稳定性与有效性。

2. 问题叙述

2.1. 一阶常微分方程组

一般地，一阶常微分方程组可表示为如下形式

$y^{\prime }\left(x\right)=f\left(x,y\right),y\left(x_0\right)=y_0,x\in \left[x_0,b\right]\subset ℝ$ , (2.1)

其中，未知向量函数，给定右端函数向量和初始条件分别为

$y\left(x\right)=\left(\begin{array}{c}{y}_1\left(x\right)\\ y_2\left(x\right)\\ ⋮\\ y_m\left(x\right)\end{array}\right),f\left(x,y\right)=\left(\begin{array}{c}{f}_1\left(x,y_1,y_2,\cdots ,y_m\right)\\ f_2\left(x,y_1,y_2,\cdots ,y_m\right)\\ ⋮\\ f_m\left(x,y_1,y_2,\cdots ,y_m\right)\end{array}\right),y_0=\left(\begin{array}{c}{y}_{1,0}\\ y_{2,0}\\ ⋮\\ y_{m,0}\end{array}\right)$ .

定义一阶常微分方程组雅可比矩阵为

$f_y\left(x,y\right)={\left[\frac{\partial f_i\left(x,y_1,\cdots ,y_m\right)}{\partial y_j}\right]}_{i,j=1}^m$ .

根据微分中值定理，对于任意分量函数$f_i\left(x,y\right)$ ，存在某个点${\xi }_i$ 介于y与z之间，使得

$f\left(x,y\right)-f\left(x,z\right)={\displaystyle {\sum }_{j=1}^m\frac{\partial f_i\left(x,{\xi }_i\right)}{\partial y_j}\left(y_j-z_j\right)}$ ,

由此可得

$\left|f\left(x,y\right)-f\left(x,z\right)\right|\le {\displaystyle {\sum }_{j=1}^m\left|\frac{\partial f_i\left(x,{\xi }_i\right)}{\partial y_j}\right|\left|y_j-z_j\right|}$ . (2.2)

对向量$f\left(x,y\right)-f\left(x,z\right)$ 取无穷范数，有

${\Vert f\left(x,y\right)-f\left(x,z\right)\Vert }_{\infty }=\underset{1\le i\le m}{\max }\left|f_i\left(x,y\right)-f_i\left(x,z\right)\right|$ , (2.3)

结合式(2.2)与(2.3)，可进一步推得

${\Vert f\left(x,y\right)-f\left(x,z\right)\Vert }_{\infty }\le \underset{1\le i\le m}{\max }{\displaystyle {\sum }_{j=1}^m\left|\frac{\partial f_i\left(x,{\xi }_i\right)}{\partial y_j}\right|{\Vert y-z\Vert }_{\infty }}$ .

若函数$f\left(x,y\right)$ 关于$y$ 满足Lipschitz条件，即存在常数$k$ ，使得对任意$y,z\in {ℝ}^m$ 有

${\Vert f\left(x,y\right)-f\left(x,z\right)\Vert }_{\infty }\le k\cdot {\Vert y-z\Vert }_{\infty },\forall y,z\in {ℝ}^m$ .

则方程组(2.1)存在唯一解。Lipschitz常数$k$ 可通过雅可比矩阵的元素进行估计，其上界可取为

$k=\underset{t_0\le x\le b}{\max }\left(\underset{1\le i\le m}{\max }\left(\underset{y\in {ℝ}^m}{\sup }{\displaystyle {\sum }_{j=1}^m\left|\frac{\partial f_i\left(x,y\right)}{\partial y_j}\right|}\right)\right)$ .

特别地，一阶线性常微分方程组的形式为

$\frac{\text{d}{y}_i}{\text{d}x}={\displaystyle {\sum }_{j=1}^m{a}_{ij}\left(x\right)\cdot y_j}+f_i\left(x\right),i=1,2,\cdots ,m$ , (2.4)

其中，$A\left(x\right)={\left(a_{ij}\left(x\right)\right)}_{m\times m}$ 是系数矩阵，$f\left(x\right)={\left(f_1\left(x\right),\cdots ,f_m\left(x\right)\right)}^{\text{T}}$ 为给定的非齐次项向量函数。该方程组可表示为向量形式

$\frac{\text{d}y}{\text{d}x}=A\left(x\right)\cdot y+f\left(x\right)$ .

若${\Vert A\left(x\right)\cdot y-A\left(x\right)\cdot z\Vert }_{\infty }\le k\cdot {\Vert y-z\Vert }_{\infty }$ 成立，则

$k=\underset{x_0\le x\le b}{\max }\left(\underset{1\le i\le m}{\max }{\displaystyle \sum }_{j=1}^m\left|a_{ij}\left(x\right)\right|\right).$

2.2. 高阶线性常微分方程

一般地，高阶线性常微分方程表示形式如下

$y^{\left(m\right)}+a_1\left(x\right)\cdot y^{\left(m-1\right)}+\cdots +a_{m-1}\left(x\right)\cdot y^{\prime }+a_m\left(x\right)\cdot y=f\left(x\right),x\in \left[x_0,b\right].$ (2.5)

定义变量为

$y_1=y,y_2=y^{\prime },\cdots ,y_m=y^{\left(m-1\right)}.$

则原方程可以写成一阶方程组

$\frac{\text{d}y}{\text{d}x}=A\left(x\right)\cdot y+f\left(x\right)=F\left(x,A\left(x\right),f\left(x\right),y\left(x\right)\right),$ (2.6)

其中，系数矩阵$A\left(x\right)$ 的具体形式为

$A\left(x\right)={\left(\begin{array}{ccccc}0& 1& 0& \cdots & 0\\ 0& 0& 1& \cdots & 0\\ ⋮& ⋮& ⋮& & ⋮\\ 0& 0& 0& \cdots & 1\\ -a_m\left(x\right)& -a_{m-1}\left(x\right)& -a_{m-2}\left(x\right)& \cdots & -a_1\left(x\right)\end{array}\right)}_{m\times m},$

$f\left(x\right)$ 为向量形式的非齐次项，仅最后一行为$f\left(x\right)$ ，其余分量为零。

为了保证原方程(2.5)所对应一阶方程组(2.6)解的存在性与唯一性，需要满足Lipschitz条件，即对任意$y,z\in {ℝ}^m$ ，有

${\Vert F\left(x,A\left(x\right),f\left(x\right),y\left(x\right)\right)-F\left(x,A\left(x\right),f\left(x\right),z\left(x\right)\right)\Vert }_{\infty }\le k\cdot {\Vert y-z\Vert }_{\infty },$

$k=\max \left\{1,{\displaystyle {\sum }_{i=1}^m\left|a_i\left(x\right)\right|}\right\}$ .

3. 显式Euler格式数值分析

3.1. 显式Euler格式局部截断误差

显式欧拉格式用于求解一阶方程组(1.13)，其格式为

$y_{n+1}=y_n+\text{Δ}x\left(A\left(x_n\right)\cdot y_n+f\left(x_n\right)\right),$ (3.1)

其中，$\Delta x=x_{n+1}-x_n=h$ 为步长，$A\left(x_n\right)$ 和$f\left(x_n\right)$ 分别表示$A\left(x\right)$ 和$f\left(x\right)$ 在$x_n$ 取值。

定理1 若方程(2.6)右端项$F\left(x,A\left(x\right),f\left(x\right),y\left(x\right)\right)$ 关于$y$ 满足Lipschitz条件，假设第$n$ 步精确成立，即$y\left(x_n\right)=y_n$ ，则显式Euler法的局部截断误差

$y\left(x_{n+1}\right)-y_{n+1}=O\left(h^2\right)$ .

证明：将$y\left(x_{n+1}\right)$ 在$x_n$ 处泰勒展开为：

$y\left(x_{n+1}\right)=y\left(x_n\right)+y\text{'}\left(x_n\right)\cdot \Delta x+\frac{1}{2}\Delta x^2{y}^{″}\left({\theta }^n\right)$ ,

其中，$y^{\prime }\left(x_n\right)={\left({y^{\prime }}_1\left(x_n\right),{y^{\prime }}_2\left(x_n\right),\cdots ,{y^{\prime }}_m\left(x_n\right)\right)}^{\text{T}}$ ，$y^{″}\left({\theta }_n\right)={\left({y^{″}}_1\left({\theta }_1^n\right),{y^{″}}_2\left({\theta }_2^n\right),\cdots ,{y^{″}}_m\left({\theta }_m^n\right)\right)}^{\text{T}}$ ， ${\theta }^n={\left({\theta }_1^n,{\theta }_2^n,\cdots ,{\theta }_m^n\right)}^{\text{T}}\in {ℝ}^m$ ，${\theta }_i^n\in \left(x_n,x_{n+1}\right)$ 。又由方程(2.1)，则局部截断误差为

$\begin{array}{l}y\left(x_{n+1}\right)-y_{n+1}\\ =y\left(x_n\right)+y^{\prime }\left(x_n\right)\cdot \Delta x+\frac{1}{2}\Delta x^2{y}^{″}\left({\theta }^n\right)-\left(y_n+\Delta x\left(A\left(x_n\right)\cdot y_n+f\left(x_n\right)\right)\right)\\ =\frac{1}{2}\Delta x^2{y}^{″}\left({\theta }^n\right)=O\left(h^2\right).\end{array}$

因此，显式Euler格式的局部截断误差是$O\left(h^2\right)$ 。

3.2. 显式Euler格式全局误差

用计算机进行计算时，由于计算机的字长有限，无法完整显示所有数据，并且计算过程有可能会产生新的误差，这些误差称为舍入误差。本文在舍入误差都可控的情况下，分析显示Euler格式的全局误差。即存在正常数$\gamma$ ，使得每步的舍入误差向量${\eta }^n\in {ℝ}^m$ 满足${\text{sup}}_n{\Vert {\eta }^n\Vert }_{\infty }\le \gamma$ ，则有下述定理。

定理2 若方程(2.6)右端项$F\left(x,A\left(x\right),f\left(x\right),y\left(x\right)\right)$ 关于$y$ 满足Lipschitz条件，并且满足${\text{max}}_{a\le x\le b}{\Vert y^{\prime }\left(x\right)\Vert }_{\infty }\le M$ ，记$x_{n+1}=x_0+\left(n+1\right)\Delta x$ ，定义全局误差为

$E_{n+1}=y\left(x_{n+1}\right)-{\tilde{y}}_{n+1},$

其中$y\left(x_{n+1}\right)$ 是精确解，${\tilde{y}}_{n+1}$ 是考虑了舍入误差的数值解，则有

${\Vert E_{n+1}\Vert }_{\infty }\le {\text{e}}^{k\left(x_{n+1}-x_0\right)}{\Vert E_0\Vert }_{\infty }+\left({\text{e}}^{k\left(x_{n+1}-x_0\right)}-1\right)\left(\frac{\Delta xM}{2k}+\frac{\gamma }{\Delta xk}\right).$

证明：根据泰勒公式，函数$y\left(x\right)$ 在$x_{n+1}=x_n+\Delta x$ 处展开为

$y\left(x_{n+1}\right)=y\left(x_n\right)+\Delta x\cdot y^{\prime }\left(x_n\right)+\frac{{\left(\Delta x\right)}^2}{2}{y}^{″}\left({\xi }^n\right),$

其中${\xi }^n={\left({\xi }_1^n,{\xi }_2^n,\cdots ,{\xi }_m^n\right)}^{\text{T}}\in {ℝ}^m,{\xi }_i^n\in \left(x_n,x_{n+1}\right)$ ，又因为 $y^{\prime }\left(x\right)=F\left(x,A\left(x\right),f\left(x\right),y\right)$ ，所以有

$y\left(x_{n+1}\right)=y\left(x_n\right)+\Delta x\cdot F\left(x_n,A\left(x_n\right),f\left(x_n\right),y\left(x_n\right)\right)+\frac{{\left(\Delta x\right)}^2}{2}y″\left({\xi }^n\right).$

显式Euler格式的数值解${\tilde{y}}_{n+1}$ 为

${\tilde{y}}_{n+1}={\tilde{y}}_n+\Delta x\cdot F\left(x_n,A\left(x_n\right),f\left(x_n\right),{\tilde{y}}_n\right)+{\eta }^n,$

定义全局误差$E_{n+1}=y\left(x_{n+1}\right)-{\tilde{y}}_{n+1}$ ，$E_n=y\left(x_n\right)-{\tilde{y}}_n$ ，代入可得

$\begin{array}{c}{E}_{n+1}=y\left(x_n\right)+\Delta x\cdot F\left(x_n,A\left(x_n\right),f\left(x_n\right),y\left(x_n\right)\right)+\frac{{\left(\Delta x\right)}^2}{2}{y}^{″}\left({\xi }^n\right)+\frac{{\left(\Delta x\right)}^2}{2}{y}^{″}\left({\xi }^n\right)\\ -\left({\tilde{y}}_n+\Delta x\cdot F\left(x_n,A\left(x_n\right),f\left(x_n\right),{\tilde{y}}_n\right)+{\eta }^n\right)\\ =E_n+\Delta x\left(F\left(x_n,A\left(x_n\right),f\left(x_n\right),y\left(x_n\right)\right)-F\left(x_n,A\left(x_n\right),f\left(x_n\right),{\tilde{y}}_n\right)\right)+\frac{{\left(\Delta x\right)}^2}{2}{y}^{″}\left({\xi }^n\right)-{\eta }^n.\end{array}$

因为$F\left(x,A\left(x\right),f\left(x\right),y\right)$ 关于$y$ 满足Lipschitz条件，即

${\Vert F\left(x_n,A\left(x_n\right),f\left(x_n\right),y\left(x_n\right)\right)-F\left(x_n,A\left(x_n\right),f\left(x_n\right),{\widehat{y}}_n\right)\Vert }_{\infty }\le k{\Vert y\left(x_n\right)-{\widehat{y}}_n\Vert }_{\infty }=k{\Vert E_n\Vert }_{\infty },$

对$E_{n+1}$ 取无穷范数，根据范数性质可得

$\begin{array}{c}{\Vert E_{n+1}\Vert }_{\infty }\le {\Vert E_n\Vert }_{\infty }+\Delta x{\Vert F\left(x_n,A\left(x_n\right),f\left(x_n\right),y\left(x_n\right)\right)-F\left(x_n,A\left(x_n\right),f\left(x_n\right),{\widehat{y}}_n\right)\Vert }_{\infty }\\ +\frac{{\left(\Delta x\right)}^2}{2}{\Vert y^{″}\left({\xi }^n\right)\Vert }_{\infty }+{\Vert {\eta }^n\Vert }_{\infty }\\ \le {\Vert E_n\Vert }_{\infty }+\Delta x\cdot k{\Vert E_n\Vert }_{\infty }+\frac{{\left(\Delta x\right)}^2}{2}{\Vert y^{″}\left({\xi }^n\right)\Vert }_{\infty }+{\Vert {\eta }^n\Vert }_{\infty },\end{array}$

已知${\text{max}}_{a\le x\le b}{\Vert y^{\prime }\left(t\right)\Vert }_{\infty }\le M$ 且${\sup }_n{\Vert {\eta }^n\Vert }_{\infty }\le \gamma$ ，则

${\Vert E_{n+1}\Vert }_{\infty }\le \left(1+\Delta x\cdot k\right){\Vert E_n\Vert }_{\infty }+\frac{{\left(\Delta x\right)}^2}{2}M+\gamma ,$

令$\alpha =1+\Delta x\cdot k$ ，$\beta =\frac{{\left(\Delta x\right)}^2}{2}M+\gamma$ ，则有${\Vert E_{n+1}\Vert }_{\infty }\le \alpha {\Vert E_n\Vert }_{\infty }+\beta$ ，通过多次递推可得

${\Vert E_{n+1}\Vert }_{\infty }\le {\alpha }^{n+1}{\Vert E_0\Vert }_{\infty }+\left({\displaystyle \sum }_{i=0}^n{\alpha }^i\right)\beta ,$

由等比数列求和公式${\displaystyle {\sum }_{i=0}^n{\alpha }^i}=\frac{{\alpha }^{n+1}-1}{\alpha -1}$ ($\alpha \ne 1$ )，且$\alpha -1=\Delta x\cdot k$ ，所以

${\Vert E_{n+1}\Vert }_{\infty }\le {\alpha }^{n+1}{\Vert E_0\Vert }_{\infty }+\frac{{\alpha }^{n+1}-1}{\Delta x\cdot k}\beta .$

因为$\alpha =1+\Delta x\cdot k$ ，根据重要极限${\text{lim}}_{x\to 0}{\left(1+x\right)}^{\frac{1}{x}}=\text{e}$ ，当$\Delta x$ 较小时，有${\left(1+\Delta x\cdot k\right)}^{\frac{1}{\Delta x\cdot k}}\approx \text{e}$ 。

又因为$\left(n+1\right)\Delta x=x_{n+1}-x_0$ ，所以

${\alpha }^{n+1}={\left(1+\Delta x\cdot k\right)}^{n+1}={\left(1+\Delta x\cdot k\right)}^{\frac{\left(n+1\right)\Delta xk}{\Delta t}}\approx {\text{e}}^{k\left(x_{n+1}-x_0\right)}$ ，将其代入上式可得

${\Vert E_{n+1}\Vert }_{\infty }\le {\text{e}}^{k\left(x_{n+1}-x_0\right)}{\Vert E_0\Vert }_{\infty }+\frac{{\text{e}}^{k\left(x_{n+1}-x_0\right)}-1}{\Delta x\cdot k}\beta ,$

把$\beta =\frac{{\left(\Delta x\right)}^2}{2}M+\gamma$ 代入即得到定理结果：

${\Vert E_{n+1}\Vert }_{\infty }\le {\text{e}}^{k\left(x_{n+1}-x_0\right)}{\Vert E_0\Vert }_{\infty }+\left({\text{e}}^{k\left(x_{n+1}-x_0\right)}-1\right)\left(\frac{\Delta xM}{2k}+\frac{\gamma }{\Delta xk}\right).$

定理3 若定理2的条件成立，若不考虑初始误差和数据的舍入误差，即$E_0=0$ ，$\gamma =0$ 时，则显式Euler格式的全局误差为

${\Vert E_{n+1}\Vert }_{\infty }\le \left({\text{e}}^{k\left(x_{n+1}-x_0\right)}-1\right)\frac{\Delta xM}{2k}.$

证明：将$E_0=0$ ，$\gamma =0$ 代入定理2的结论

${\Vert E_{n+1}\Vert }_{\infty }\le {\text{e}}^{k\left(x_{n+1}-x_0\right)}\times 0+\left({\text{e}}^{k\left(x_{n+1}-x_0\right)}-1\right)\left(\frac{\Delta xM}{2k}+\frac{0}{\Delta x\cdot k}\right),$

化简后可得

${\Vert E_{n+1}\Vert }_{\infty }\le \left({\text{e}}^{k\left(x_{n+1}-x_0\right)}-1\right)\frac{\Delta xM}{2k}.$

3.3. 显式Euler格式的渐进稳定性

定义1 称Euler方法是渐进稳定的，如果存在常数$C$ 和$h_0$ ，使得$\forall h=\Delta x$ 满足$0<h<h_0$ 以及区间$\left[a,b\right]$ 内的离散点$x_n=a+nh$ ，Euler方法的解$y_n$ 和$z_n$ 满足以下不等式：

${\Vert y_n-z_n\Vert }_{\infty }\le C\cdot {\Vert y_0-z_0\Vert }_{\infty }.$

定理4 在定理2的假设条件下，显式Euler法是渐进稳定的。

证明：考虑Euler方法的迭代公式：

$y_{n+1}=y_n+\text{Δ}x\cdot F\left(x_n,A\left(x_n\right),f\left(x_n\right),y_n\right),$

和

$z_{n+1}=z_n+\text{Δ}x\cdot F\left(x_n,A\left(x_n\right),f\left(x_n\right),z_n\right),$

两式相减，利用Lipschitz条件，可得

$\begin{array}{c}{\Vert y_{n+1}-z_{n+1}\Vert }_{\infty }\le \left(1+\Delta x\cdot k\right)\cdot {\Vert y_n-z_n\Vert }_{\infty }\\ \le {\left(1+\Delta x\cdot k\right)}^{n+1}\cdot {\Vert y_0-z_0\Vert }_{\infty }.\end{array}$

利用指数函数的性质${\left(1+h\cdot k\right)}^{n+1}\le {\text{e}}^{\left(n+1\right)\cdot h\cdot k}$ ，进一步推导

$\begin{array}{c}{\Vert y_{n+1}-z_{n+1}\Vert }_{\infty }\le {\text{e}}^{\left(n+1\right)\cdot h\cdot k}\cdot {\Vert y_0-z_0\Vert }_{\infty }\\ \le {\text{e}}^{\left(b-a\right)\cdot k}\cdot {\Vert y_0-z_0\Vert }_{\infty }.\end{array}$

从而，当$a\le x_n=a+nh\le b$ 时，令$C={\text{e}}^{\left(b-a\right)\cdot K}$ 有

${\Vert y_n-z_n\Vert }_{\infty }\le {\text{e}}^{\left(b-a\right)\cdot k}\cdot {\Vert y_0-z_0\Vert }_{\infty }=C{\Vert y_0-z_0\Vert }_{\infty }.$

故对$\forall h\in \left(0,h_0\right)$ ，显式Euler法的解满足渐进稳定性条件，且连续依赖于初值，因此，显式Euler法是渐进稳定的。

4. 数值实验

4.1. 常系数二阶线性常微分方程

$\{\begin{array}{l}{y}^{″}-2y^{\prime }+2y={\text{e}}^{2x}\sin x,0\le x\le 1\\ y\left(0\right)=-0.4,y^{\prime }\left(0\right)=-0.6\end{array}$ .

设$y_1=y,y_2=y^{\prime }$ ，则原方程可转化为一阶方程组

$\{\begin{array}{l}{y^{\prime }}_1=y_2\\ {y^{\prime }}_2=2y_2-2y_1+{\text{e}}^{2x}\sin x\end{array}$

记

$\overline{y}=\left(\begin{array}{c}{y}_1\\ y_2\end{array}\right),A=\left(\begin{array}{cc}0& 1\\ -2& 2\end{array}\right),\overline{f}\left(x\right)=\left(\begin{array}{c}0\\ {\text{e}}^{2x}\sin x\end{array}\right)$.

则上述方程可写为一阶线性常微分方程组

${\overline{y}}^{\prime }=A\overline{y}+\overline{f}\left(x\right),{\overline{y}}_0=\left(\begin{array}{c}{y}_1\left(0\right)\\ y_2\left(0\right)\end{array}\right)=\left(\begin{array}{c}-0.4\\ -0.6\end{array}\right).$

将区间$\left[0,1\right]$ 等分为N个子区间，设步长为$\Delta x=\frac{1}{N}$ ，记$x_n=nh$ 。显式Euler格式为

${\overline{y}}_{n+1}={\overline{y}}_n+\Delta x\left(A{\overline{y}}_n+\overline{f}\left(x_n\right)\right),n=0,1,\cdots ,N-1.$

其中

${\overline{f}}_n=\overline{f}\left(x_n\right)=\left(\begin{array}{c}0\\ {\text{e}}^{2x_n}\sin x_n\end{array}\right),{\overline{y}}_n\approx \overline{y}\left(x_n\right)=\left(\begin{array}{c}{y}_n^1\\ y_n^2\end{array}\right),$

递推过程为

$\begin{array}{l}{\overline{y}}_1={\overline{y}}_0+\Delta x\left(A{\overline{y}}_0+{\overline{f}}_0\right)\\ {\overline{y}}_2={\overline{y}}_1+\Delta x\left(A{\overline{y}}_1+{\overline{f}}_1\right)\\ ⋮\\ {\overline{y}}_{n+1}={\overline{y}}_n+h\cdot \left(A{\overline{y}}_n+{\overline{f}}_n\right)\end{array}$

因此，${\overline{y}}_n$ 的第一个分量$y_n^1$ 即为原始常微分方程在$x_n$ 点的数值解。

图1展示了显式Euler法在不同步数(N = 10, 100, 1000, 10,000)下的数值解与精确解的对比情况。可以明显观察到，随着步数的增加(即步长$\Delta x$ 减小，$\Delta x=\frac{1}{N}$ )，显式Euler法得到的数值解逐渐逼近精确解。当N = 1000时，数值解几乎与精确解重合，表明在足够小的步长条件下，显式Euler法能够较为准确地逼近精确解。

Figure 1. Comparison of numerical solution and exact solution of explicit Euler method under non-synchronization number

图1. 显式Euler法数值解与精确解在不同步数下的对比

当N = 1000时，精确解与数值解的结果：

Table 1. Comparison of numerical solution and exact solution at N = 1000

表1. 数值解与精确解在N = 1000下的对比结果

x	精确解	数值解
0.0	−0.40000	−0.40000
0.1	−0.46173	−0.46172
0.2	−0.52556	−0.52556
0.3	−0.58860	−0.58866
0.4	−0.64661	−0.64678
0.5	−0.69356	−0.69394
0.6	−0.72115	−0.72185
0.7	−0.71815	−0.71934
0.8	−0.66971	−0.67160
0.9	−0.55644	−0.55933
1.0	−0.35339	−0.35764

为验证显式Euler法的有效性，本文选取步数N = 1000，将数值解与精确解关键节点处进行对比。由表1可见，数值解与精确解在各节点处误差极小，显式Euler法能够很好地逼近精确解。

Table 2. Infinite norm of error vector under asynchronous numbers

表2. 不同步数下误差向量的无穷范数

N	${\Vert �差\Vert }_{\infty }=\underset{1\le n\le N}{\max }\left|y_1\left(x_n\right)-y_n^1\right|$
50	0.08153
100	0.04165
200	0.02105
400	0.01058
800	0.00531

为了考察显式Euler法收敛性，本文选取不同的步数(N = 50, 100, 200, 400, 800)进行计算，并记录相应的误差向量的无穷范数(即数值解与精确解在各离散节点上的最大偏差)。结果如表2所示，随着步数N的增加，步长随之减小，误差呈现出单调递减的趋势。值得注意的是，误差几乎呈线性下降的规律，即当步数加倍(步长减半)时，误差无穷范数近似减小为原来的一半。例如，从N = 100增加至N = 200时，误差由0.04165下降至0.02105，几乎减少一半。因此，显式Euler格式保持一阶收敛，这与定理3的理论分析结果一致。

4.2. 变系数二阶线性常微分方程

$\{\begin{array}{l}{y}^{″}-\frac{1}{x}{y}^{\prime }=0,1\le x\le 2\\ y\left(1\right)=1,y^{\prime }\left(1\right)=2\end{array}$

设$y_1=y,y_2=y^{\prime }={y^{\prime }}_1$ ，则原方程可转化为一阶方程组

$\{\begin{array}{l}{y^{\prime }}_1=y_2\\ {y^{\prime }}_2=\frac{1}{x}{y}_2\end{array}$

记

$\overline{y}=\left(\begin{array}{c}{y}_1\\ y_2\end{array}\right),{\overline{y}}^{\prime }=\left(\begin{array}{c}{y^{\prime }}_1\\ {y^{\prime }}_2\end{array}\right),$

构造变系数矩阵

$A\left(x\right)=\left(\begin{array}{cc}0& 1\\ 0& \frac{1}{x}\end{array}\right).$

则上述方程可写为一阶线性常微分方程组

${\overline{y}}^{\prime }=A\left(x\right)\overline{y},\overline{y}\left(x_0\right)={\overline{y}}_0=\left(\begin{array}{c}y\left(1\right)\\ y^{\prime }\left(1\right)\end{array}\right)=\left(\begin{array}{c}1\\ 2\end{array}\right)$.

将区间$\left[1,2\right]$ 等分为N个子区间，设步长为$\Delta x=\frac{1}{N}$ ，记$x_n=nh$ 。显式Euler格式为

${\overline{y}}_{n+1}={\overline{y}}_n+\Delta xA\left(x_n\right){\overline{y}}_n,$

其中${\overline{y}}_n\approx \overline{y}\left(x_n\right)=\left(\begin{array}{c}{y}_n^1\\ y_n^2\end{array}\right)$ 。

递推过程为

$\begin{array}{l}{\overline{y}}_1={\overline{y}}_0+\Delta xA\left(x_0\right){\overline{y}}_0\\ {\overline{y}}_2={\overline{y}}_1+\Delta xA\left(x_1\right){\overline{y}}_1\\ ⋮\\ {\overline{y}}_{n+1}={\overline{y}}_n+\Delta xA\left(x_n\right){\overline{y}}_n\end{array}$

因此，${\overline{y}}_n$ 的第一个分量$y_n^1$ 即为原始变系数二阶微分方程在$x_n$ 点的数值解。

图2展示了当不同步数(N = 10, 20, 50, 100, 200, 400, 10,000)时所得数值解与精确解之间的对比情况。从图2中可以看出，随着步数的逐渐增加(即步长减小)，显式Euler法得到的数值解逐渐逼近精确解。当N ≥ 200时，数值解与精确解几乎重合，说明该方法在步长足够小时能够有效逼近真实解。

Figure 2. Comparison of numerical solution and exact solution of explicit Euler method under different steps N

图2. 不同步数N下显式Euler法数值解与精确解的对比

当N = 200时，精确解与数值解的结果：

Table 3. The comparison between the numerical solution and the exact solution of the Euler method when N = 200

表3. N = 200时显示Euler法数值解与精确解的对比

x	数值解	精确解
1.0	1.0000	1.00000
1.1	1.20980	1.21000
1.2	1.43481	1.44000
1.3	1.68940	1.69000
1.4	1.95361	1.96000
1.5	2.24900	2.25000
1.6	2.55880	2.56000
1.7	2.88860	2.89000
1.8	3.23840	3.24000
1.9	3.60061	3.61000
2.0	3.99800	4.00000

为验证显式Euler法的有效性，本文选取步数N = 200，即步长$h=\frac{1}{200}$ ，并与该初值问题的精确解进行比较。由表3可见，数值解与精确解在各节点处误差极小，最大误差不超过0.002，说明显式Euler法具有良好的数值精度。

当N取不同值时，误差向量的无穷范数如下表所示：

Table 4. Infinite norm of numerical error under different steps N

表4. 不同步数N下数值误差无穷范数

N	${\Vert �差\Vert }_{\infty }=\underset{1\le n\le N}{\max }\left|y_1\left(x_n\right)-y_n^1\right|$
50	0.02000
100	0.01000
200	0.00500
400	0.00250
800	0.00125

以不同的步数(N = 50, 100, 200, 400, 800)进行数值计算，记录数值解与精确解的最大误差。结果如表4所示，误差向量的无穷范数随着步数的增大而减小，呈现线性收敛趋势，即当步数加倍，误差缩小为原来的一半。例如，当N从50增长至100时，误差由0.02000减小至0.01000。因此，显式Euler格式是线性收敛的，即一阶收敛，这与理论分析结果一致(见定理3)。

5. 结论

本文围绕高阶线性常微分方程，分析了显式Euler法在舍入误差与非舍入误差下的误差估计，并讨论了Euler格式的渐进稳定性。在数值实验部分，分别选取常系数与变系数高阶线性常微分方程，结合不同步长对显式Euler格式的数值误差与渐进稳定性进行了验证。结果表明，显式Euler法在步长足够小的情况下具有良好的渐进稳定性，数值误差随步长减小呈线性收敛，符合理论分析结果，数值实验验证了显式Euler法在数值求解高阶线性常微分方程的稳定性和有效性。

显式Euler法具有结构简单、实现便捷的优势，但仅有一阶数值精度。在后续研究中，将进一步研究梯形格式、Runge-Kutta格式等，提高高阶常微分方程数值解的收敛速度。

基金项目

重庆科技大学2024年度市级大学生科技创新训练项目(S2024115001029)；重庆市教育委员会科学技术研究项目(KJQN202301521)。

NOTES

^*通讯作者。