常微分方程的多区域配置方法

A Multi-Domain Collocation Method of Ordinary Differential Equations

摘要: 以等距节点为插值节点，构造常微分方程的Lagrange插值逼近算法格式，将常微分方程转化成矩阵方程求解。通过将区域分解的方式提高算法精度，数值实验证明本文所提算法的高精度，此方法可以广泛应用到其他常微分方程的求解中。

Abstract: Using equidistant nodes as interpolation nodes, construct a Lagrange interpolation approximation algorithm format for ordinary differential equations, and transform the ordinary differential equations into matrix equations for solution. By decomposing regions, the algorithm accuracy is improved, and numerical experiments have demonstrated the high accuracy of the proposed algorithm. This method can be widely applied to the solution of other ordinary differential equations.

1. 引言

常微分方程是现代数学的一个重要分支，在物理学，微分几何，计算数学，计算机图形学，图象处理以及大量的边缘科学诸如电磁流体力学、化学流体力学、动力气象学、海洋动力学、地下水动力学等学科中都有许多重要的应用。因此，对于常微分方程的求解就显得尤为重要。对于常微分方程的精确解求解办法有很多种方法 [1] 。对于常微分方程的数值求解方法 [2] 也有很多，有学者采用Milstein方法 [3] ，也有学者采用Birkhoff配点法 [4] ，还有学者采用龙格–库塔法，利用泰勒展开式求解二阶常微分方程 [5] 。最近，有人利用Lagrange插值 [6] 函数的优点，采用等距节点来求解偏微分方程 [7] 单区域时空方向上的数值解，也有学者对KdV方程采用Legendre-Hermite时空二元谱配置方法 [8] ，采用非等距节点求解偏微分方程，而文献 [9] 给出多区域非等距节点谱配法求解偏微分方程。而本文将考虑对常微分方程采用等距节点构造Lagrange插值算法，并给出多区域等距节点谱配法求解常微分方程的算法格式。

我们以下面的二阶常微分方程为例，记 $\frac{\partial }{\partial x}={\partial }_x$ ，考虑常微分方程的初边值问题：

$\{\begin{array}{l}{\partial }_x^{2}u+2{\partial }_{x}u+u=e^{-x},x\in \left(-1,1\right),\\ u\left(-1\right)=2e,u\left(1\right)=0.\end{array}$ (1)

先通过构造单区域上的Lagrange插值逼近算法来逼近上述方程的精确解，再将区间(−1, 1)分解成两部分，构造多区域上的Lagrange插值逼近算法，将方程转化成矩阵方程求解。通过数值实验发现，该算法所求的近似解能够很好的逼近精确解，而且该算法实施较为容易，易于理解。

2. 基于等距节点的拉格朗日多项式的微分矩阵推导

记 $h=\left(b-a\right)/N$ ， $x_0=a$ ， $x_k=x_0+kh$ ， $k=0,1,\cdots ,N$ 。则以 $x_k$ 为节点的Lagrange插值基函数为：

${\phi }_n\left(x\right)=\frac{\left(x-x_0\right)\left(x-x_1\right)\cdots \left(x-x_{n-1}\right)\left(x-x_{n+1}\right)\cdots \left(x-x_N\right)}{\left(x_n-x_0\right)\left(x_n-x_1\right)\cdots \left(x_n-x_{n-1}\right)\left(x_n-x_{n+1}\right)\cdots \left(x_n-x_N\right)},n=0,1,\cdots ,N.$

令 $\omega \left(x\right)=\left(x-x_0\right)\left(x-x_1\right)\cdots \left(x-x_N\right)$ ，则有：

${\phi }_n\left(x\right)=\frac{\omega \left(x\right)}{\left(x-x_n\right){\partial }_x\omega \left(x_n\right)}.$

记 $P_N\left(a,b\right)$ 为次数 $\le N$ 的多项式集合，对于 $u\left(x\right)\in C\left(a,b\right)$ ，其插值多项式为：

$P_N\left(x\right)={\displaystyle \underset{n=0}{\overset{N}{\sum }}{u}_n{\phi }_n}\left(x\right),u_n=u\left(x_n\right),x\in \left[a,b\right].$

记 $P_N\left(a,b\right)$ 关于x求一阶导数，并令 $x=x_k,k=0,1,2,\cdots ,N$ ，得：

${\partial }_x{P}_N\left(x_k\right)={\displaystyle \underset{m=0}{\overset{N}{\sum }}{\partial }_x{\phi }_n\left(x_k\right)u_n}.$

设 $D=\left(d_{kn}\right)$ 是 $\left(N+1\right)\times \left(N+1\right)$ 矩阵，且 $d_{kn}={\partial }_x{\phi }_n\left(x_k\right)$ ，由文献 [7] 可知：

$d_{kn}=\{\begin{array}{l}\frac{{\left(-1\right)}^{n-k}k!\left(N-k\right)!}{h\left(k-n\right)n!\left(N-n\right)!},k\ne n,\\ \frac{1}{h}\left(\left(1+\frac{1}{2}+\cdots +\frac{1}{n}\right)-\left(1+\frac{1}{2}+\cdots +\frac{1}{N-n}\right)\right),0<k=n<N,\\ -\frac{1}{h}\left(1+\frac{1}{2}+\cdots +\frac{1}{N}\right),k=n=0,\\ \frac{1}{h}\left(1+\frac{1}{2}+\cdots +\frac{1}{N}\right),k=n=N.\end{array}$

对于 $d_{kn}$ 的具体推导在文献 [7] 中有详细推导，由文献 [10] 可知，m阶微分矩阵与一阶微分矩阵的关系是： $D^m=D^{\left(m\right)}$ 。

3. 常微分方程多区域Lagrange插值逼近算法格式

令 $\Omega =\left(-1,1\right)$ ，(1)式的Lagrange插值逼近方法就是求多项式 $u_M\left(x\right)\in P_M\left(\Omega \right)$ 满足：

$\{\begin{array}{l}{\partial }_x^2{u}_M+2{\partial }_x{u}_M+u_M=e^{-x},x\in \left(-1,1\right),\\ u_M\left(-1\right)=2e,u_M\left(1\right)=0.\end{array}$ (2)

其数值解可以写成：

$u_M\left(x\right)={\displaystyle \underset{n=0}{\overset{M}{\sum }}{\hat{u}}_n}{\phi }_n\left(x\right),{\hat{u}}_n=\hat{u}\left(x_n\right).$

3.1. 单区域

3.1.1. 单区域常微分方程Lagrange插值逼近算法构造

将数值解带入(2)中有：

$\{\begin{array}{l}{\displaystyle \underset{n=0}{\overset{M}{\sum }}{\hat{u}}_n}{\partial }_x^2{\phi }_n\left(x_k\right)+2{\displaystyle \underset{n=0}{\overset{M}{\sum }}{\hat{u}}_n}{\partial }_n{\phi }_n\left(x_k\right)+{\displaystyle \underset{n=0}{\overset{M}{\sum }}{\hat{u}}_n}{\phi }_n\left(x_k\right)=e^{-x_k},k=1,2,\cdots ,M-1,\\ {\displaystyle \underset{n=0}{\overset{M}{\sum }}{\hat{u}}_n}{\phi }_n\left(x_0\right)=2e,{\displaystyle \underset{n=0}{\overset{M}{\sum }}{\hat{u}}_n}{\phi }_n\left(x_M\right)=0.\end{array}$ (3)

由微分矩阵定义可以将(3)转化为(4)：

$\{\begin{array}{l}{\displaystyle \underset{n=0}{\overset{M}{\sum }}{\hat{u}}_n}{d}_{k,n}^{\left(2\right)}+2{\displaystyle \underset{n=0}{\overset{M}{\sum }}{\hat{u}}_n}{d}_{k,n}+{\hat{u}}_k=e^{-x_k},k=1,2,\cdots ,M-1,\\ {\hat{u}}_0=2e,{\hat{u}}_M=0.\end{array}$ (4)

则(4)为：

$\{\begin{array}{l}{\displaystyle \underset{n=1}{\overset{M-1}{\sum }}{\hat{u}}_n}{d}_{k,n}^{\left(2\right)}+{\hat{u}}_0{d}_{k,0}^{\left(2\right)}+{\hat{u}}_M{d}_{k,M}^{\left(2\right)}+2{\displaystyle \underset{n=1}{\overset{M-1}{\sum }}{\hat{u}}_n}{d}_{k,n}+{\hat{u}}_k=e^{-x_k}-2{\hat{u}}_0{d}_{k,0}2{\hat{u}}_M{d}_{k,M},k=1,2,\cdots ,M-1,\\ {\hat{u}}_0=2e,{\hat{u}}_M=0.\end{array}$ (5)

取 $k=1,2,\cdots ,M-1$ ，(5)式可以写成：

$\begin{array}{l}\left(\begin{array}{ccc}{d}_{1,1}^{\left(2\right)}& \cdots & d_{1,M-1}^{\left(2\right)}\\ ⋮& \ddots & ⋮\\ d_{M-1,1}^{\left(2\right)}& \cdots & d_{M-1,M-1}^{\left(2\right)}\end{array}\right)\left(\begin{array}{c}{\hat{u}}_1\\ ⋮\\ {\hat{u}}_{M-1}\end{array}\right)+2\left(\begin{array}{ccc}{d}_{1,1}& \cdots & d_{1,M-1}\\ ⋮& \ddots & ⋮\\ d_{M-1,1}& \cdots & d_{M-1,M-1}\end{array}\right)\left(\begin{array}{c}{\hat{u}}_1\\ ⋮\\ {\hat{u}}_{M-1}\end{array}\right)+\left(\begin{array}{c}{\hat{u}}_1\\ ⋮\\ {\hat{u}}_{M-1}\end{array}\right)\\ =\left(\begin{array}{c}{e}^{-x_1}\\ ⋮\\ e^{-x_{M-1}}\end{array}\right)-\left({\hat{u}}_0\left(\begin{array}{c}{d}_{1,0}^{\left(2\right)}\\ ⋮\\ d_{M-1,0}^{\left(2\right)}\end{array}\right)+{\hat{u}}_M\left(\begin{array}{c}{d}_{1,M}^{\left(2\right)}\\ ⋮\\ d_{M-1,M}^{\left(2\right)}\end{array}\right)\right)-2\left({\hat{u}}_0\left(\begin{array}{c}{d}_{1,0}\\ ⋮\\ d_{M-1,0}\end{array}\right)+{\hat{u}}_M\left(\begin{array}{c}{d}_{1,M}\\ ⋮\\ d_{M-1,M}\end{array}\right)\right).\end{array}$

我们不妨记：

$A=\left(\begin{array}{ccc}{d}_{1,1}^{\left(2\right)}& \cdots & d_{1,M-1}^{\left(2\right)}\\ ⋮& \ddots & ⋮\\ d_{M-1,1}^{\left(2\right)}& \cdots & d_{M-1,M-1}^{\left(2\right)}\end{array}\right),B=\left(\begin{array}{ccc}{d}_{1,1}& \cdots & d_{1,M-1}\\ ⋮& \ddots & ⋮\\ d_{M-1,1}& \cdots & d_{M-1,M-1}\end{array}\right),$

$X=\left(\begin{array}{c}{\hat{u}}_1\\ ⋮\\ {\hat{u}}_{M-1}\end{array}\right),D={\hat{u}}_0\left(\begin{array}{c}{d}_{1,0}^{\left(2\right)}\\ ⋮\\ d_{M-1,0}^{\left(2\right)}\end{array}\right)+{\hat{u}}_M\left(\begin{array}{c}{d}_{1,M}^{\left(2\right)}\\ ⋮\\ d_{M-1,M}^{\left(2\right)}\end{array}\right),$

$C=\left(\begin{array}{c}{e}^{-x_1}\\ ⋮\\ e^{-x_{M-1}}\end{array}\right),F=2\left({\hat{u}}_0\left(\begin{array}{c}{d}_{1,0}\\ ⋮\\ d_{M-1,0}\end{array}\right)+{\hat{u}}_M\left(\begin{array}{c}{d}_{1,M}\\ ⋮\\ d_{M-1,M}\end{array}\right)\right).$

我们可以得到如下矩阵方程：

$AX+2BX+X=C-D-F.$ (6)

3.1.2. 单区域数值结果

我们将(6)式转化为：

$\left(A+2B+E\right)X=C-D-F.$ (7)

其中E表示 $M-1$ 阶单位阵，则(7)式可以求得近似解为：

$X={\left(A+2B+E\right)}^{-1}\left(C-F-D\right).$

我们用 $L^{\infty }$ -来衡量精确解与数值解之间的误差：

$E_M=\underset{1\le k\le M-1}{\max }\left|u\left(x_k\right)-u_M\left(x_k\right)\right|.$

而方程(1)的精确解为： $u\left(x\right)=\frac{1}{2}{\left(x-1\right)}^2{e}^{-x}$ 。

可以得到误差图(图1)，对最大误差 $E_M$ 取对数可以得到随着插值节点数M的增大，最大误差 $E_M$ 呈下降趋势，且误差效果良好，可以说明上述所提算法格式良好，具有高精度。

3.2. 多区域

3.2.1. 区域常微分方程的Lagrange插值逼近算法构造

为了能够得到更高误差精度，我们将整个求解区域 $\left(-1,1\right)$ 分解成两个区域 ${\Omega }_1=\left(-1,c\right)$ ， ${\Omega }_2=\left(c,1\right)$ 。数值解展开为：

${u_M\left(x\right)|}_{{\Omega }_i}=u_M^i\left(x\right)={\displaystyle \underset{n=0}{\overset{M_i}{\sum }}{u}_M^{\left(i\right)}}{\phi }_n\left(x\right),x\in {\Omega }_i,i=1,2.$

Figure 1. Single area error plot

图1. 单区域误差图

由(1)的解 $u\in H^2\left(\Omega \right)$ ，在公共边界处有数值解的函数值和导数值相等，即：

函数值： ${\displaystyle \underset{n=0}{\overset{M_1}{\sum }}{\hat{u}}_n^{\left(1\right)}}{\phi }_n\left(x_{M_1}\right)={\displaystyle \underset{n=0}{\overset{M_2}{\sum }}{\hat{u}}_n^{\left(2\right)}}{\phi }_n\left({\tilde{x}}_0\right)$ 。

导数值： ${\displaystyle \underset{n=0}{\overset{M_1}{\sum }}{\hat{u}}_n^{\left(1\right)}}{\partial }_x{\phi }_n\left(x_{M_1}\right)={\displaystyle \underset{n=0}{\overset{M_2}{\sum }}{\hat{u}}_n^{\left(2\right)}{\partial }_x}{\phi }_n\left({\tilde{x}}_0\right)$ 。

由微分矩阵定义，将上述两式联立可以得到：

${\hat{u}}_{M_1}^{\left(1\right)}={\hat{u}}_0^{\left(2\right)}=\frac{{\displaystyle \underset{n=1}{\overset{M_2}{\sum }}{\hat{u}}_n^{\left(2\right)}{\tilde{d}}_{0,n}}-{\displaystyle \underset{n=0}{\overset{M_1-1}{\sum }}{\hat{u}}_n^{\left(1\right)}{d}_{M_1,n}}}{d_{M_1,M_1}-{\tilde{d}}_{0,0}}.$ (8)

其中 $d_{k,m}$ 表示区域 ${\Omega }_1$ 上的微分矩阵， ${\tilde{d}}_{k,m}$ 表示区域 ${\Omega }_2$ 上的微分矩阵， $x_k\in {\Omega }_1,k=1,2,\cdots ,M_1$ ，

${\tilde{x}}_k\in {\Omega }_2,k=1,2,\cdots ,M_2$ 。

在区域 ${\Omega }_1$ 用数值解逼近区域 ${\Omega }_1$ 的精确解得到：

$\{\begin{array}{l}{\displaystyle \underset{n=0}{\overset{M_1}{\sum }}{\hat{u}}_n^{\left(1\right)}}{\partial }_x^2{\phi }_n\left(x_k\right)+2{\displaystyle \underset{n=0}{\overset{M_1}{\sum }}{\hat{u}}_n^{\left(1\right)}}{\partial }_n{\phi }_n\left(x_k\right)+{\displaystyle \underset{n=0}{\overset{M_1}{\sum }}{\hat{u}}_n^{\left(1\right)}}{\phi }_n\left(x_k\right)=e^{-x_k},k=1,2,\cdots ,M_1-1,\\ {\displaystyle \underset{n=0}{\overset{M_1}{\sum }}{\hat{u}}_n^{\left(1\right)}}{\phi }_n\left(x_0\right)=2e.\end{array}$ (9)

对(9)利用微分矩阵定义并将端点值提出可以得到：

$\begin{array}{l}{\displaystyle \underset{n=1}{\overset{M_1-1}{\sum }}{\hat{u}}_n^{\left(1\right)}}{d}_{k,n}^{\left(2\right)}+{\hat{u}}_0^{\left(1\right)}{d}_{k,0}^{\left(2\right)}+{\hat{u}}_{M_1}^{\left(1\right)}{d}_{k,M_1}^{\left(2\right)}+2{\displaystyle \underset{n=1}{\overset{M_1-1}{\sum }}{\hat{u}}_n^{\left(1\right)}}{d}_{k,n}+{\hat{u}}_k^{\left(1\right)}\\ =e^{-x_k}-\left(2{\hat{u}}_0^{\left(1\right)}{d}_{k,0}+2{\hat{u}}_{M_1}^{\left(1\right)}{d}_{k,M_1}\right),k=1,2,\cdots ,M_1-1.\end{array}$ (10)

将(8)式带入到(10)式中有：

$\begin{array}{l}{\displaystyle \underset{n=1}{\overset{M_1-1}{\sum }}{\hat{u}}_n^{\left(1\right)}}{d}_{k,n}^{\left(2\right)}+2{\displaystyle \underset{n=1}{\overset{M_1-1}{\sum }}{\hat{u}}_n^{\left(1\right)}}{d}_{k,n}+{\hat{u}}_k^{\left(1\right)}+\frac{{\hat{u}}_{M_2}^{\left(2\right)}{\tilde{d}}_{0,M_2}-{\hat{u}}_0^{\left(1\right)}{d}_{M_1,0}}{d_{M_1,M_1}-{\tilde{d}}_{0,0}}\left(2d_{k,M_1}+d_{k,M_1}^{\left(2\right)}\right)\\ =e^{-x_k}-\frac{{\displaystyle \underset{n=1}{\overset{M_2-1}{\sum }}{\hat{u}}_n^{\left(2\right)}{\tilde{d}}_{0,n}}-{\displaystyle \underset{n=1}{\overset{M_1-1}{\sum }}{\hat{u}}_n^{\left(1\right)}{d}_{M_1,n}}}{d_{M_1,M_1}-{\tilde{d}}_{0,0}}\left(2d_{k,M_1}+d_{k,M_1}^{\left(2\right)}\right)-{\hat{u}}_0^{\left(1\right)}\left(d_{k,0}^{\left(2\right)}+2d_{k,0}\right),k=1,2,\cdots ,M_1-1.\end{array}$ (11)

将(11)式展开，并记：

$G=\left(\begin{array}{ccc}{d}_{1,1}^{\left(2\right)}& \cdots & d_{1,M_1-1}^{\left(2\right)}\\ ⋮& \ddots & ⋮\\ d_{M_1-1,1}^{\left(2\right)}& \cdots & d_{M_1-1,M_1-1}^{\left(2\right)}\end{array}\right),H=\left(\begin{array}{ccc}{d}_{1,1}& \cdots & d_{1,M_1-1}\\ ⋮& \ddots & ⋮\\ d_{M_1-1,1}& \cdots & d_{M_1-1,M_1-1}\end{array}\right),$

$X_1=\left(\begin{array}{c}{\hat{u}}_1^{\left(1\right)}\\ ⋮\\ {\hat{u}}_{M_1-1}^{\left(1\right)}\end{array}\right),X_2=\left(\begin{array}{c}{\hat{u}}_1^{\left(2\right)}\\ ⋮\\ {\hat{u}}_{M_2-1}^{\left(2\right)}\end{array}\right),F_1=\left(\begin{array}{c}{d}_{1,0}^{\left(2\right)}+2d_{1,0}\\ ⋮\\ d_{M_1-1,0}^{\left(2\right)}+2d_{M_1-1,0}\end{array}\right),$

$E_1=\left(\begin{array}{ccc}{d}_{M_1,1}& \cdots & d_{M_1,M_1-1}\end{array}\right),E_2=\left(\begin{array}{ccc}{\tilde{d}}_{0,1}& \cdots & {\tilde{d}}_{0,M_1-1}\end{array}\right),$

$C_1=\left(\begin{array}{c}{e}^{-x_1}\\ ⋮\\ e^{-x_{M_1-1}}\end{array}\right),D_1=\left(\begin{array}{c}{d}_{1,M_1}^{\left(2\right)}+2d_{1,M_1}\\ ⋮\\ d_{M_1-1,M_1}^{\left(2\right)}+2d_{M_1-1,M_1}\end{array}\right).$

由(1)式中边界值可知， ${\hat{u}}_0^{\left(1\right)}=2e$ ， ${\hat{u}}_{M_2}^{\left(2\right)}=0$ 。(11)式可以转化成下列矩阵方程：

$G{X}_1+2H{X}_1+X_1-\frac{2e{d}_{M_1,0}}{d_{M_1,M_1}-{\tilde{d}}_{0,0}}{D}_1=C_1-\frac{E_2{X}_2-E_1{X}_1}{d_{M_1,M_1}-{\tilde{d}}_{0,0}}{D}_1-2e{F}_1.$ (12)

同理，对区域 ${\Omega }_2$ 仿照区域 ${\Omega }_1$ 可以得到：

$\{\begin{array}{l}{\displaystyle \underset{n=0}{\overset{M_2}{\sum }}{\hat{u}}_n^{\left(2\right)}}{\partial }_x^2{\phi }_n\left({\tilde{x}}_k\right)+2{\displaystyle \underset{n=0}{\overset{M_2}{\sum }}{\hat{u}}_n^{\left(2\right)}}{\partial }_n{\phi }_n\left({\tilde{x}}_k\right)+{\displaystyle \underset{n=0}{\overset{M_2}{\sum }}{\hat{u}}_n^{\left(2\right)}}{\phi }_n\left({\tilde{x}}_k\right)=e^{-{\tilde{x}}_k},k=1,2,\cdots ,M_2-1,\\ {\displaystyle \underset{n=0}{\overset{M_2}{\sum }}{\hat{u}}_n^{\left(2\right)}}{\phi }_n\left({\tilde{x}}_{M_2}\right)=0.\end{array}$ (13)

有微分矩阵定义以及(8)式可以得到(14)式：

$\begin{array}{l}{\displaystyle \underset{n=1}{\overset{M_2-1}{\sum }}{\hat{u}}_n^{\left(2\right)}}{\tilde{d}}_{k,n}^{\left(2\right)}+2{\displaystyle \underset{n=1}{\overset{M_2-1}{\sum }}{\hat{u}}_n^{\left(2\right)}}{\tilde{d}}_{k,n}+{\hat{u}}_k^{\left(2\right)}+\frac{{\hat{u}}_{M_2}^{\left(2\right)}{\tilde{d}}_{0,M_2}-{\hat{u}}_0^{\left(1\right)}{d}_{M_1,0}}{d_{M_1,M_1}-{\tilde{d}}_{0,0}}\left(2{\tilde{d}}_{k,0}+{\tilde{d}}_{k,0}^{\left(2\right)}\right)\\ =e^{-{\tilde{x}}_k}-\frac{{\displaystyle \underset{n=1}{\overset{M_2-1}{\sum }}{\hat{u}}_n^{\left(2\right)}{\tilde{d}}_{0,n}}-{\displaystyle \underset{n=1}{\overset{M_1-1}{\sum }}{\hat{u}}_n^{\left(1\right)}{d}_{M_1,n}}}{d_{M_1,M_1}-{\tilde{d}}_{0,0}}\left(2{\tilde{d}}_{k,0}+{\tilde{d}}_{k,0}^{\left(2\right)}\right)-{\hat{u}}_{M_2}^{\left(2\right)}\left({\tilde{d}}_{k,M_2}^{\left(2\right)}+2{\tilde{d}}_{k,M_2}\right),k=1,2,\cdots ,M_2-1.\end{array}$ (14)

将(14)式展开，并记：

$P=\left(\begin{array}{ccc}{\tilde{d}}_{1,1}^{\left(2\right)}& \cdots & {\tilde{d}}_{1,M_2-1}^{\left(2\right)}\\ ⋮& \ddots & ⋮\\ {\tilde{d}}_{M_2-1,1}^{\left(2\right)}& \cdots & {\tilde{d}}_{M_2-1,M_2-1}^{\left(2\right)}\end{array}\right),Q=\left(\begin{array}{ccc}{\tilde{d}}_{1,1}& \cdots & {\tilde{d}}_{1,M_2-1}\\ ⋮& \ddots & ⋮\\ {\tilde{d}}_{M_2-1,1}& \cdots & {\tilde{d}}_{M_2-1,M_2-1}\end{array}\right),$

$C_2=\left(\begin{array}{c}{e}^{-{\tilde{x}}_1}\\ ⋮\\ e^{-{\tilde{x}}_{M_2-1}}\end{array}\right),D_2=\left(\begin{array}{c}{\tilde{d}}_{1,0}^{\left(2\right)}+2{\tilde{d}}_{1,0}\\ ⋮\\ {\tilde{d}}_{M_2-1,0}^{\left(2\right)}+2{\tilde{d}}_{M_2-1,0}\end{array}\right),F_2=\left(\begin{array}{c}{\tilde{d}}_{1,M_2}^{\left(2\right)}+2{\tilde{d}}_{1,M_2}\\ ⋮\\ {\tilde{d}}_{M_2-1,M_2}^{\left(2\right)}+2{\tilde{d}}_{M_2-1,M_2}\end{array}\right).$

将(14)式转化成如下矩阵方程：

$P{X}_2+2Q{X}_2+X_2-\frac{2e{d}_{M_1,0}}{d_{M_1,M_1}-{\tilde{d}}_{0,0}}{D}_2=C_2-\frac{E_2{X}_2-E_1{X}_1}{d_{M_1,M_1}-{\tilde{d}}_{0,0}}{D}_2.$ (15)

3.2.2. 多区域数值结果

将(12)和(15)式联立可以得到：

$\begin{array}{l}\left(\begin{array}{cc}G+2H& 0\\ 0& P+2Q\end{array}\right)\left(\begin{array}{c}{X}_1\\ X_2\end{array}\right)+\left(\begin{array}{c}{X}_1\\ X_2\end{array}\right)-\frac{2e{d}_{M_1,0}}{d_{M_1,M_1}-{\tilde{d}}_{0,0}}\left(\begin{array}{c}{D}_1\\ D_2\end{array}\right)\\ =\left(\begin{array}{c}{C}_1\\ C_2\end{array}\right)-2e\left(\begin{array}{c}{F}_1\\ 0\end{array}\right)-\frac{1}{d_{M_1,M_1}-{\tilde{d}}_{0,0}}\left(\begin{array}{c}{D}_1\\ D_2\end{array}\right)\left(\begin{array}{cc}-E_1& E_2\end{array}\right)\left(\begin{array}{c}{X}_1\\ X_2\end{array}\right).\end{array}$ (16)

(16)可以等价于：

$\begin{array}{l}\left(\left(\begin{array}{cc}G+2H& 0\\ 0& P+2Q\end{array}\right)+I+\frac{1}{d_{M_1,M_1}-{\tilde{d}}_{0,0}}\left(\begin{array}{c}{D}_1\\ D_2\end{array}\right)\left(\begin{array}{cc}-E_1& E_2\end{array}\right)\right)\left(\begin{array}{c}{X}_1\\ X_2\end{array}\right)\\ =\left(\begin{array}{c}{C}_1\\ C_2\end{array}\right)-2e\left(\begin{array}{c}{F}_1\\ 0\end{array}\right)+\frac{2e{d}_{M_1,0}}{d_{M_1,M_1}-{\tilde{d}}_{0,0}}\left(\begin{array}{c}{D}_1\\ D_2\end{array}\right)\end{array}$

不妨令：

$\begin{array}{l}R=\left(\left(\begin{array}{cc}G+2H& 0\\ 0& P+2Q\end{array}\right)+I+\frac{1}{d_{M_1,M_1}-{\tilde{d}}_{0,0}}\left(\begin{array}{c}{D}_1\\ D_2\end{array}\right)\left(\begin{array}{cc}-E_1& E_2\end{array}\right)\right),\\ S=\left(\begin{array}{c}{C}_1\\ C_2\end{array}\right)-2e\left(\begin{array}{c}{F}_1\\ 0\end{array}\right)+\frac{2e{d}_{M_1,0}}{d_{M_1,M_1}-{\tilde{d}}_{0,0}}\left(\begin{array}{c}{D}_1\\ D_2\end{array}\right),Y=\left(\begin{array}{c}{X}_1\\ X_2\end{array}\right).\end{array}$

则(16)可以写成下列方程：

$Y=R^{-1}S.$ (17)

其中I为 $M_1+M_2-2$ 阶单位矩阵。利用(17)式可以求出(1)式的数值解，我们令 $c=0.1$ ，仍然用用 $L^{\infty }$ -来衡量精确解与数值解之间的误差，其中 $E_{M_1}$ 表示区域 ${\Omega }_1$ 上的误差， $E_{M_2}$ 表示区域 ${\Omega }_2$ 上的误差， $E_M$ 表示整个区域上的误差, $M=M_1+M_2$ ：

$\begin{array}{l}{E}_{M_1}=\underset{1\le k\le M_1-1}{\max }\left|u_1\left(x_k\right)-u_{M_1}\left(x_k\right)\right|\\ E_{M_2}=\underset{1\le k\le M_2-1}{\max }\left|u_2\left(x_k\right)-u_{M_2}\left(x_k\right)\right|\\ E_M=\max \left(E_{M_1},E_{M_2}\right).\end{array}$

图2是多区域误差图，可以看出，误差量级与 $M$ 成线性关系，最大误差 $E_M$ 随着插值节点数 $M$ 的增大而减小，可以看到，多区域的误差精度能达到−12，明显比单区域高，由此可以说明多区域的Lagrange插值效果比单区域的Lagrange插值效果好,也体现了本文所提算法格式有效。

Figure 2. Multi region error plot

图2. 多区域误差图

4. 结论

本文通过构造Lagrange插值算法逼近格式求解常微分方程的数值解，将常微分方程转化成矩阵方程求解，通过数值实验结果可以得出，本文所提算法格式良好。为了进一步提高算法精度，通过将区域分解的方式，可以提高误差精度，数值结果也证明了这一结论。本文所提算法格式简单，易于操作，计算所用的节点个数少即可达到较高精度，对于求解其他常微分方程也适用。

基金项目

河南省大学生创新创业训练计划项目(No.202310464051)。

参考文献