常微分方程组初值问题的Lagrange插值逼近方法

Lagrange Interpolation Approximation Method for Initial Value Problems of Ordinary Differential Equations

摘要: 常微分方程的快速发展和很多学科有着紧密的关系，随着众多的数学家对其研究的不断加深，极大地推动了现代数学的发展。因此，本文研究拉格朗日插值(Lagrange)逼近方法在常微分方程组初值问题上的应用，推导出Lagrange插值法的逼近算法，求解常微分方程组初值问题的数值解，也就是通过具体的例子，利用Lagrange插值逼近方法构造相对应的数值格式，寻找误差和多项式次数的关系，将结果可视化，最后对相应的误差结果进行分析。数值结果表明Lagrange插值逼近方法具有较高的精度。

Abstract: The rapid development of ordinary differential equations is closely related to many disciplines, and with the continuous deepening of the research of many mathematicians, the development of modern mathematics has been greatly promoted. Therefore, this paper studies the application of the Lagrange interpolation approximation method in the initial value problem of ordinary differential equations, deduces the approximation algorithm of the Lagrange interpolation approximation method, and solves the numerical solution of the initial value problem of ordinary differential equations. The numerical results show that the Lagrange interpolation approximation method has high accuracy.

1. 引言

常系数微分方程是微分方程在社会上应用的较为基础的内容，即本文所说的常微分方程。常微分方程(ODE)是指微分方程的自变量只有一个的方程 [1] ，常可被用于刻画自然界中很多事物的运动规律，在自然科学和社会科学等各个方面应用广泛。因此，对常微分方程初值问题求解方法的研究显得尤为重要。

近年来，对常微分方程已经有一些理论研究：文献 [2] 利用卷积逼近和Bihari不等式等工具，在函数满足一定条件下，提出并证明了常微分方程初值问题的解存在且唯一；文献 [3] 用单调迭代方法给出了二阶常微分方程处置问题的解的存在性唯一性。在常微分方程中初值问题一直都是人们学习研究的方向，常用于求解常微分方程中初值问题的数值解的解法有Runge-Kudta法、Euler法和Adas法等 [4] [5] [6] ，学者们对数值解法做了许多研究：文献 [7] 研究了具有无穷多个导数的常微分方程的初值问题，给出了该问题的求解方法；文献 [8] 这里讨论了简化后的带奇异系数的二阶常微分方程初值问题的有限元解和它的收敛性；文献 [9] 研究了Nystromn方法在一类奇异初值问题上的应用；文献 [10] 研究了运用Euler方法求解二阶非线性微分方程。在上述研究的基础上，研究利用Lagrange插值逼近方法来求解常微分方程组初值问题。

考虑下面的常微分方程的初值问题：

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=2x+3y\left(t\right)+f_1\left(t\right),t\in \left[0,1\right),f_1\left(t\right)=5t,\\ \frac{\text{d}y}{\text{d}t}=3x+2y+f_2\left(t\right),t\in \left(0,1\right]f_2\left(t\right)=8e^t,\\ x\left(1\right)=2-\frac{13}{5}-3e,y\left(0\right)=\frac{17}{5}.\end{array}$ (1)

针对(1)式，基于相关的理论知识，利用Lagrange插值法构造相对应的数值格式，寻找误差和多项式次数的关系，将结果可视化，并对相应的误差结果进行分析。

2. 预备知识

2.1. Lagrange插值方法

法国一位世界闻名的数学家拉格朗日，发明了一种多项式和函数插值计算方法，即为人们所熟知的拉格朗日插值法。利用拉格朗日插值法就可得到某个多项式，能够正好地从每个观察的节点取到所观察到的最大值，通过坐标构建插值函数，与原函数相结合进而求得近似已知曲线的函数。

设在节点 $a\le x_0<x_1<\cdots <x_n\le b$ 处，函数值 $y_i=f\left(x_i\right),i=0,1,\cdots ,n$ ，寻求函数 $p\left(x\right)$ 使其满足 $p\left(x_i\right)=y_i,i=0,1,2,\cdots ,n$ 。为了求解函数 $p\left(x\right)$ ，因此，拉格朗日插值法被提出。

拉格朗日插值多项式的定义：

将不同插值节点处的函数值与相应插值基函数的乘积求和，构成一个多项式 $L_n\left(x\right)$ ，我们称 $L_n\left(x\right)$ 为拉格朗日插值多项式。

n次插值拉格朗日多项式可表示为：

$L_n\left(x\right)={\displaystyle \underset{k=0}{\overset{n}{\sum }}{y}_k{l}_k\left(x\right)}$ (2)

其中： $l_k\left(x\right)$ 称为n阶插值基函数，表达式如下

$l_k(x)=\frac{\left(x-x_0\right)\cdots \left(x-x_{k-1}\right)\left(x-x_{k+1}\right)\cdots \left(x-x_n\right)}{\left(x_k-x_0\right)\cdots \left(x_k-x_{k-1}\right)\left(x_k-x_{k+1}\right)\cdots \left(x_k-x_n\right)}\begin{array}{cc},& k=0,1,\cdots ,n\end{array}$

并且满足：

$l_k\left(x_i\right)=\{\begin{array}{c}1\begin{array}{cc},& k=i,\end{array}\\ 0\begin{array}{cc},& k\ne i.\end{array}\end{array}$

显然， $l_k\left(x\right)$ 是一个n次多项式。

引入 ${\omega }_{n+1}\left(x\right)=\left(x-x_0\right)\left(x-x_1\right)\cdots \left(x-x_n\right)$ ，则n次插值多项式为

$L_n\left(x\right)={\displaystyle \underset{k=0}{\overset{n}{\sum }}{y}_k\frac{{\omega }_{n+1}\left(x\right)}{\left(x-x_k\right){\omega }_{x+1}\left(x_k\right)}}.$

形如(2)式的插值多项式 $L_n\left(x\right)$ 称为拉格朗日插值多项式。

2.2. 基于等距节点的Lagrange插值多项式的微分矩阵

在定区间 $\left[c,d\right]$ ，记 $p=\left(d-c\right)/N$ ， $x_0=c$ ， $x_j=x_0+kp$ ， $k=0,1,\cdots ,N$ 。对于任意次数以 $x_j$ ， $j=0,1,\cdots ,N$ 为节点的N次插值基函数为：

$l_j\left(x\right)=\frac{\omega \left(x\right)}{\left(x-x_j\right){\partial }_{x}w\left(x_j\right)},j=0,1,\cdots ,N$

其中 $\omega \left(x\right)=\left(x-x_0\right)\left(x-x_1\right)\cdots \left(x-x_N\right)$ 。

其相应的插值多项式为：

$L_n\left(x\right)={\displaystyle \underset{j=0}{\overset{N}{\sum }}{y}_j{l}_j\left(x\right)},x\in \left[c,d\right]$

对 $L_n\left(x\right)$ 关于x求一阶导数，并令 $x=x_k,k=0,1,\cdots ,N$ 得

${\partial }_x{L}_n\left(x_k\right)={\displaystyle \underset{j=0}{\overset{N}{\sum }}{\partial }_x{l}_j\left(x_k\right)y_j}$

由文献 [11] 可知：

设 $D=d_{kj}$ 是 $\left(N+1\right)\times \left(N+1\right)$ 矩阵，且 $d_{kj}={\partial }_x{l}_j\left(x_k\right)$ ，则有：

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

那么 $L_n\left(x_k\right)$ 的一阶导数可表示为

${\partial }_x{L}_n\left(x_k\right)={\displaystyle \underset{j=0}{\overset{N}{\sum }}{d}_{kj}{y}_j}$

再对 $L_n\left(x_k\right)$ 求二阶导数，可得：

${\partial }_x^2{L}_n\left(x_k\right)={\displaystyle \underset{j=0}{\overset{N}{\sum }}{d}_{kj}^2{y}_j},d_{kj}^{(2)}={\partial }_x^2{l}_j\left(x_k\right).$

这里 $\hat{D}=\left({\hat{D}}_{kj}\right)=D^2$ 。

3. 常微分方程组初值问题的数值方法

在这一节中，针对常微分方程组初值问题，利用Lagrange插值逼近方法构造相应的数值格式，求得方程组近似解和相应精确解之间的误差情况，并且将通过数值实验所得的误差结果进行分析。

3.1. Lagrange插值逼近方法的数值算法

本节给出问题(1)的Lagrange插值逼近方法的数值算法。对区间 $\left(0,1\right)$ 做等距划分， $t_j=jh$ ， $h=\frac{1}{N}$ ， $j=0,1,\cdots ,N$ 。记 ${\varphi }_j\left(t\right)$ 为 $t\in \left(0,1\right)$ 上的插值基函数，问题(1)的相应的数值解分别为：

$\begin{array}{l}{p}_N\left(t\right)={\displaystyle \underset{j=0}{\overset{N}{\sum }}{x}_j}\left(t_j\right){\varphi }_j\left(t\right)={\displaystyle \underset{j=0}{\overset{N}{\sum }}{x}_j}{\varphi }_j\left(t\right),x_j=x\left(t_j\right),\\ q_N\left(t\right)={\displaystyle \underset{j=0}{\overset{N}{\sum }}y}\left(t_j\right){\varphi }_j^{}\left(t\right)={\displaystyle \underset{i=0}{\overset{N}{\sum }}{y}_j}{\varphi }_j\left(t\right),y_j=y\left(t_j\right).\end{array}$ (3)

用数值解(3)式来逼近问题(1)的函数 $x\left(t\right)$ 和 $y\left(t\right)$ 的精确解，将(3)式代入到(1)式中可得：

$\{\begin{array}{l}\frac{\text{d}{p}_N\left(t\right)}{\text{d}t}=2p_N\left(t\right)+3q_N\left(t\right)+5t,t=t_k,t\in \left[0,1\right),k=0,\cdots N-1,\\ \frac{\text{d}{q}_N\left(t\right)}{\text{d}t}=2p_N\left(t\right)+3q_N\left(t\right)+8e^t,t=t_k,t\in \left(0,1\right],k=1,\cdots ,N,\\ p_N\left(t_N\right)=p_N\left(1\right)=2-\frac{13}{5}-3e,q_N\left(t_0\right)=q_N\left(0\right)=\frac{17}{5}.\end{array}$ (4)

将 $p_N\left(t\right)$ 和 $q_N\left(t\right)$ 的表达式代入到(4)式中可得：

$\{\begin{array}{l}{\displaystyle \underset{j=0}{\overset{N}{\sum }}{x}_j}\frac{\text{d}{\varphi }_j\left(t_k\right)}{\text{d}t}=2{\displaystyle \underset{j=0}{\overset{N}{\sum }}{x}_j}{\varphi }_j\left(t_k\right)+3{\displaystyle \underset{j=0}{\overset{N}{\sum }}{y}_j}{\varphi }_j\left(t_k\right)+5t,k=0,1\cdots N-1,\\ {\displaystyle \underset{j=0}{\overset{N}{\sum }}{y}_j}\frac{\text{d}{\varphi }_j\left(t_k\right)}{\text{d}t}=3{\displaystyle \underset{j=0}{\overset{N}{\sum }}{x}_j}{\varphi }_j\left(t_k\right)+2{\displaystyle \underset{j=0}{\overset{N}{\sum }}{y}_j}{\varphi }_j\left(t_k\right)+8e^t,k=1,2\cdots ,N\\ {\displaystyle \underset{j=1}{\overset{N}{\sum }}{x}_j}{\varphi }_j\left(t_N\right)=-\frac{3}{5}-3e,{\displaystyle \underset{j=0}{\overset{N}{\sum }}{y}_j}{\varphi }_j\left(t_0\right)=\frac{17}{5},\end{array}$ (5)

整理上式可以得

$\{\begin{array}{l}{\displaystyle \underset{j=0}{\overset{N-1}{\sum }}{d}_{kj}}{x}_j=2x_k+3y_k+5t_k-x_N{d}_{k,N},k=0,1,\cdots ,N-1,\\ {\displaystyle \underset{j=1}{\overset{N}{\sum }}{d}_{kj}}{y}_j=3x_k+2y_k+8e^{t_k}-y_0{d}_{k,0},k=1,2,\cdots ,N,\\ x_N=2-\frac{13}{5}-3e,y_0=\frac{17}{5}.\end{array}$ (6)

其中： $d_{kj}=\frac{\text{d}{\varphi }_j\left(t_k\right)}{\text{d}t},\left(k=0,1,2,\cdots ,N;j=0,1,2,\cdots ,N\right)$ 。将初值条件代入(6)式中的前面两个方程中，(6)式可等价为

$\{\begin{array}{l}{\displaystyle \underset{j=1}{\overset{N-1}{\sum }}\left(\left(d_{00}-2\right)d_{kj}-d_{k0}{d}_{0j}\right)x_j}=2\left(d_{00}-2\right)x_k+3\left(d_{00}-2\right)y_k\\ +\left(5t_k-d_{kN}{x}_N\right)\left(d_{00}-2\right)-\left(3y_0+5t_0-d_{0N}{x}_N\right)d_{k0},k=1,\cdots ,N-1,\\ {\displaystyle \underset{j=1}{\overset{N-1}{\sum }}\left(\left(d_{NN}-2\right)d_{kj}-d_{kN}{d}_{Nj}\right)y_j}=3\left(d_{NN}-2\right)x_k+2\left(d_{NN}-2\right)y_k\\ +\left(8e^{t_k}-d_{k0}{y}_0\right)\left(d_{NN}-2\right)-\left(3x_N+8e^{t_N}-d_{N0}{y}_0\right)d_{kN},k=1,\cdots N-1.\end{array}$ (7)

整理后可等价为

$\{\begin{array}{l}\left({\bar{D}}_1-2\left(d_{00}-2\right)E_1\right)X-\left(3\left(d_{00}-2\right)E_2\right)Y=F_1,\\ \left(-3\left(d_{NN}-2\right)E_1\right)X+\left({\bar{D}}_2-2\left(d_{NN}-2\right)E_2\right)Y=F_2.\end{array}$ (8)

其中： $E_1$ 和 $E_2$ 是 $N-1$ 阶单位矩阵，

$\begin{array}{l}{x}_0=\frac{1}{d_{00}-2}\left(3y_0+5t_0-d_{0N}{x}_N-{\displaystyle \underset{j=1}{\overset{N-1}{\sum }}{d}_{0j}}{x}_j\right),\\ y_N=\frac{1}{d_{NN}-2}\left(3x_N+8e^{t_N}-d_{N0}{y}_0-{\displaystyle \underset{j=1}{\overset{N-1}{\sum }}{d}_{Nj}}{y}_j\right),\end{array}$

$\begin{array}{l}X=\left(\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_{N-1}\end{array}\right),Y=\left(\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_{N-1}\end{array}\right),\\ \overline{D_1}=\left(d_{00}-2\right)\left(\begin{array}{cccc}{d}_{11}& d_{12}& \cdots & d_{1,N-1}\\ d_{21}& d_{22}& \cdots & d_{2,N-1}\\ \cdots & \cdots & \cdots & \cdots \\ d_{N-1,1}& d_{N-1,2}& \cdots & d_{N-1,N-1}\end{array}\right)-\left(\begin{array}{c}{d}_{10}\\ d_{20}\\ ⋮\\ d_{N-10}\end{array}\right)\left(d_{01}\begin{array}{cc}{d}_{02}\cdots & d_{0N-1}\end{array}\right),\end{array}$

$\begin{array}{l}\overline{D_2}=\left(d_{NN}-2\right)\left(\begin{array}{cccc}{d}_{11}& d_{12}& \cdots & d_{1,N-1}\\ d_{21}& d_{22}& \cdots & d_{2,N-1}\\ \cdots & \cdots & \cdots & \cdots \\ d_{N-1,1}& d_{N-1,2}& \cdots & d_{N-1,N-1}\end{array}\right)-\left(\begin{array}{c}{d}_{1,N}\\ d_{2,N}\\ ⋮\\ d_{N-1,\text{N}}\end{array}\right)\left(d_{N,1}\begin{array}{cc}{d}_{N,2}\cdots & d_{N,N-1}\end{array}\right),\\ F_1=\left(d_{00}-2\right)\left(\begin{array}{c}5t_1-d_{1,N}{x}_N\\ 5t_2-d_{2,N}{x}_N\\ ⋮\\ 5t_{N-1}-d_{N-1,N}{x}_N\end{array}\right)-\left(3y_0+5t_0-d_{0,N}{x}_N\right)\left(\begin{array}{c}{d}_{10}\\ d_{20}\\ ⋮\\ d_{N-1,0}\end{array}\right),\\ F_2=\left(d_{NN}-2\right)\left(\begin{array}{c}8e^{t_1}-d_{1,0}{y}_0\\ 8e^{t_2}-d_{2,0}{y}_0\\ ⋮\\ 8e^{t_{N-1}}-d_{N-1,0}{y}_0\end{array}\right)-\left(3x_N+8e^{t_N}-d_{N,0}{y}_0\right)\left(\begin{array}{c}{d}_{1,N}\\ d_{2,N}\\ ⋮\\ d_{N-1,\text{N}}\end{array}\right).\end{array}$

整理化简(8)式可得相应的矩阵方程：

$\left(\begin{array}{c}X\\ Y\end{array}\right)={\left(\begin{array}{cc}{\bar{D}}_1-2\left(d_{00}-2\right)E_1& -3\left(d_{00}-2\right)E_2\\ -3\left(d_{NN}-2\right)E_1& {\bar{D}}_2-2\left(d_{NN}-2\right)E_2\end{array}\right)}^{-1}\left(\begin{array}{c}{F}_1\\ F_2\end{array}\right)$ (9)

3.2. 基于等距节点的Lagrange插值法的数值结果

上一节我们已经将(1)式转化为矩阵方程(8)式。在本节中，我们将用不动点迭代法求解(8)的近似解，进而求解(1)式，给出利用Lagrange插值逼近方法所求的数值解及误差分析。

问题(1)的精确解为 $x=2t-\frac{13}{5}-3e^t,y=-3t+\frac{12}{5}+e^t$ 。

用 $L-\infty$ 范数来度量数值误差，表示为：

$E_N\left(u\right)=\underset{0\le j\le N}{\max }\left(\left|u\left(x_j\right)\right|\right).$

记

$E_{1,N}=E_N\left(x\left(t\right)-p_N\left(t\right)\right),$

$E_{2,N}=E_N\left(y\left(t\right)-q_N\left(t\right)\right),$

其中 $E_{1,N}$ 和 $E_{2,N}$ 分别表示在 $L-\infty$ 范数意义下精确解 $u\left(x\right)$ ， $v\left(x\right)$ 与数值解 $p_N\left(x\right)$ ， $q_N\left(x\right)$ 之间的误差。

利用MATLAB软件进行处理，可以得到以下结果：

图1是误差对数 ${\log }_{10}{E}_{1,N}$ 随多项式次数N的变化情况。由图可知：误差对数随着多项式次数的增加，误差值巨幅变小，证明在多项式次数在十一附近拟合程度较好。另外，误差下降快速，也说明了此方法逼近效果良好。

图2是误差对数 ${\log }_{10}{E}_{2,N}$ 随多项式次数N的变化情况。由图可知：误差对数随着多项式次数N的增加呈现线性递减的趋势，误差结果达到 ${10}^{-12}$ 量级，证明在多项式次数在N = 11时附近拟合程度较好，误差下降快速，说明逼近效果良好。

4. 总结

本文研究常微分方程组初值问题的拉格朗日(Lagrange)插值逼近方法，应用到具体的例子，利用Lagrange插值逼近方法构造Lagrange插值函数和相对应的数值格式，寻找误差和多项式次数的关系，

Figure 1. The change in the logarithm of error ${\log }_{10}{E}_{1,N}$ with the degree of the polynomial N

图1. 误差对数 ${\log }_{10}{E}_{1,N}$ 随多项式次数N的变化情况

Figure 2. The change in the logarithm of error ${\log }_{10}{E}_{2,N}$ with the degree of the polynomial N

图2. 误差对数 ${\log }_{10}{E}_{2,N}$ 随多项式次数N的变化情况

利用MATLAB软件节省了大量计算过程，最后对相应的误差结果进行分析。数值实验结果表明了Lagrange插值逼近方法具有较高的精度。