高阶非线性常微分方程边值问题的Lagrange插值逼近

Lagrange Interpolation Approximation for Edge Value Problems of Higher-Order Nonlinear Ordinary Differential Equations

摘要: 本文通过选取等距节点为插值节点，利用Lagrange插值的微分矩阵，对高阶非线性常微分方程进行Lagrange插值逼近，最终转化为求解非线性方程组问题，利用不动点迭代法进行求解，并计算与解析解的误差的数量级来说明此方法的精确性，此方法对求解高阶非线性常微分方程具有重要意义。

Abstract: In this paper, equidistant nodes are selected as interpolation nodes, and the differential matrix of Lagrange interpolation is used to approximate the high-order nonlinear ordinary differential equa-tion, which is finally transformed into solving the problem of nonlinear equations. The fixed point iteration method is used to solve the problem, and the order of error between the calculation and analytical solution is calculated to illustrate the accuracy of this method. This method is of great significance to solving high-order nonlinear ordinary differential equations.

1. 引言

常微分方程边值问题(Boundary Value Problem, BVP)是数学、物理、工程等领域中的重要问题，涉及到许多实际应用，如机械、航空、天文等。总之，BVP方程是一个十分重要的研究领域，在实际中有着广泛的应用。因此求解非线性常微分方程边值问题的解对于学术研究和计算工程都有着重要的意义。对于常微分方程边值问题的求解，方法多种多样，有学者采用无波动定理，对常微分方程进行数值求解 [1] ，有学者采用非线性泛函分析法 [2] [3] ，可以求解几类高阶非线性常微分方程边值问题的正解，也有学者利用锥上的不动点指数理论，证明了一般的三阶常微分方程的正2π-周期解的存在性 [4] ，还有学者利用函数迭代法，推导出三类新的高阶非线性常微分方程的求解定理 [5] 。但在实际运用中，复杂的微分方程，本文通常不便于求出精确解，于是本文设法采用数值方法来逼近其精确解。目前，有学者利用Legendre谱配置法，来解决二阶非线性常微分方程初值问题 [6] ，也有学者研究了拉格朗日插值方法在求解热传导方程的应用 [7] 。对于较为一般的方程，通常的求解方法有有限差分法 [8] ，文献 [9] 也给出了求解二阶线性常微分方程的办法，但是未解决非线性的问题。但是以上文献的局限性在于，或只能适用于特定方程，或求解过程较为复杂，或精度不足。对于较为一般的高阶非线性常微分方程，目前求解的方法并不多。当前，有许多学者正在不断探索和创新，以提出更加高效和精确的数值方法和算法。

本文利用Lagrange插值的方法，通过构造微分矩阵，来求解高阶非线性的常微分方程，其优点有：

精度高：通过调整节点数量或节点间距，可以获得更加精确的数值解。当节点数量足够多时，这种方法的逼近解可以达到高精度的水平。

灵活性好：Lagrange插值方法，可以求解一般的非线性常微分方程，实用性强，可以适用于多种问题的求解。

计算简便：与其他数值方法相比，该方法计算量相对较小，推导过程较为简单，有利于减少求解问题所用的时间。

结果可视化：该方法可以求出问题在各个节点上的逼近解，因此可以通过直接绘制图像来表示数值结果，并可以方便地与精确解相比较。

适用范围广：本文以四阶非线性微分方程为例，对此方法进行运用，但此方法可以利用于更高阶的其他非线性形式的方程。

2. 高阶非线性常微分方程的拉格朗日插值逼近求解

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

取等距节点 $x_m=x_0+mh\left(m=0,1,\cdots ,M\right)$ ，其中h为步长，则可以得到M − 1次的Lagrange插值基函数为：

$l_m\left(x\right)=\frac{\left(x-x_0\right)\left(x-x_1\right)\cdots \left(x-x_{m-1}\right)\left(x-x_{m+1}\right)\cdots \left(x-x_M\right)}{\left(x_m-x_0\right)\left(x_m-x_1\right)\cdots \left(x_m-x_{m-1}\right)\left(x_m-x_{m+1}\right)\cdots \left(x_m-x_M\right)}$

令： $\omega \left(x\right)=\left(x-x_0\right)\left(x-x_1\right)\cdots \left(x-x_M\right)$

则有：

$l_m\left(x\right)=\frac{\omega \left(x\right)}{\left(x-x_m\right)\left({{\omega }^{\prime }\left(x\right)|}_{x=x_m}\right)}$ (1)

对于M次多项式 $P_N\left(x\right)$ ， $u\left(x\right)\in C\left(a,b\right)$ ，那么他的插值多项式可表示为：

$P_N\left(x\right)={\displaystyle \underset{m=0}{\overset{M}{\sum }}{u}_m}{l}_m(\; x\; )$

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

${P^{\prime }}_M\left(x\right)={\displaystyle \underset{m=0}{\overset{M}{\sum }}{l^{\prime }}_m}\left(x_k\right)u_m$

当 $k\ne m$ ，即 $x_k-x_m\ne 0$ 时，将(1)对x求导，得到：

${l^{\prime }}_m\left(x\right)=\frac{{\omega }^{\prime }\left(x\right)\left(x-x_m\right)-\omega \left(x\right)}{\left[{{\omega }^{\prime }\left(x\right)|}_{x=x_m}\right]{\left(x-x_m\right)}^2}$ (2)

由于：

${{\omega }^{\prime }\left(x\right)|}_{x=x_n}=\left(x_n-x_0\right)\left(x_n-x_1\right)\cdots \left(x_n-x_M\right)$

因为节点为等距节点，故有：

$x_a-x_b=\left(a-b\right)h$ (3)

令 $x=x_k$ ，并将(3)应用于(2)可得：

${l^{\prime }}_m\left(x_k\right)=\frac{{\left(-1\right)}^{m-k}k!\left(M-k\right)!}{m!\left(M-m\right)!\left(k-m\right)h}$

当 $k=m$ ，即 $x_k-x_m=0$ 时，同理可得：

${l^{\prime }}_m\left(x\right)=\frac{{\omega }^{″}\left(x_m\right)}{2{\omega }^{\prime }(\; x\; )}$

设 $d_{k,m}$ 为 $\left(M+1\right)\times \left(M+1\right)$ 矩阵，令 $d_{k,m}={l^{\prime }}_m(\; x\; k\; )$

即：

$d_{k,m}=\{\begin{array}{l}\frac{{\left(-1\right)}^{m-k}k!\left(M-k\right)!}{m!\left(M-m\right)!\left(k-m\right)h},k\ne m\hfill \\ \frac{1}{h}\left[\left(1+\frac{1}{2}+\cdots +\frac{1}{m}\right)-\left(1+\frac{1}{2}+\cdots +\frac{1}{M-m}\right)\right],0<k=m<M\hfill \\ \frac{1}{h}{\left(-1\right)}^M\left(1+\frac{1}{2}+\cdots +\frac{1}{M}\right),k=m=0\hfill \\ \frac{1}{h}\left(1+\frac{1}{2}+\cdots +\frac{1}{M}\right),k=m=M\hfill \end{array}$ (4)

这里，m阶微分矩阵 $d^{\left(m\right)}$ 与一阶微分矩阵 $d^{\left(1\right)}$ 存在以下关系 [10] ：

$d_{k,m}^{\left(m\right)}={\left(d_{k,m}^{\left(1\right)}\right)}^m$

2.2. 高阶非线性常微分方程的Lagrange插值逼近求解算法

本文以以下四阶非线性边值问题为例进行求解：

$\{\begin{array}{l}\frac{{\text{d}}^{4}u}{\text{d}{x}^4}=6{\text{e}}^{-4u}-\frac{12}{{\left(1+x\right)}^4},0<x<1\hfill \\ u\left(0\right)=0,u\left(1\right)=\ln 2\hfill \\ u^{\prime }\left(0\right)=1,u^{\prime }\left(1\right)=\frac{1}{2}\hfill \end{array}$ (5)

其精确解为 $u\left(x\right)=\ln \left(1+x\right)$ 。

$l_m\left(x\right)$ 为插值基函数，则 $u_m\left(x\right)$ 可以表示为：

$u_m\left(x\right)={\displaystyle \underset{m=0}{\overset{M}{\sum }}{\hat{u}}_m}{l}_m\left(x\right)$ (6)

将(6)带入到(5)中，并令 $x=x_k$ ，可得以下方程：

$\{\begin{array}{l}{\displaystyle \underset{m=0}{\overset{M}{\sum }}{\hat{u}}_m}{l}_m^{\left(4\right)}\left(x_k\right)=6{\text{e}}^{-4{\displaystyle \underset{m=0}{\overset{M}{\sum }}{\hat{u}}_m}{l}_m\left(x_k\right)}-\frac{12}{{\left(1+x_k\right)}^4},0<x<1\hfill \\ {\displaystyle \underset{m=0}{\overset{M}{\sum }}{\hat{u}}_m}{l}_m\left(x_0\right)=0\hfill \\ {\displaystyle \underset{m=0}{\overset{M}{\sum }}{\hat{u}}_m}{l}_m\left(x_M\right)=\ln 2\hfill \\ {\displaystyle \underset{m=0}{\overset{M}{\sum }}{\hat{u}}_m}{l^{\prime }}_m\left(x_0\right)=1\hfill \\ {\displaystyle \underset{m=0}{\overset{M}{\sum }}{\hat{u}}_m}{l^{\prime }}_m\left(x_M\right)=\frac{1}{2}\hfill \end{array}$ (7)

其中 $k,m=0,1,\cdots ,M$ 。

根据上文一阶微分矩阵(4)的描述，本文可以将方程(7)写成矩阵的形式：

$\{\begin{array}{l}{\displaystyle \underset{m=0}{\overset{M}{\sum }}{\hat{u}}_m}{d}_{k,m}^{\left(1\right)}=6{\text{e}}^{-4{\hat{u}}_m}-\frac{12}{{\left(1+x_k\right)}^4},0<x<1\hfill \\ {\hat{u}}_0=0,{\hat{u}}_M=\ln 2\hfill \\ {\hat{u}}_0{d}_{k,0}^{\left(1\right)}=1,{\hat{u}}_0{d}_{k,M}^{\left(1\right)}=\frac{1}{2}\hfill \end{array}$ (8)

其中 $k,m=0,1,\cdots ,M$ 。

此时，本文只要解矩阵方程(8)，即求出 $u_m$ ，就可以获得插值节点上的值。

2.3. 不动点迭代法求解

对于上述矩阵方程(8)，本文采用不动点迭代的方式求解非线性方程组的近似解，构造迭代格式(9)：

$d_{k,m}^{\left(4\right)}{u}_{n+1}=6{\text{e}}^{-4u_n}-12/{\left(1+X\right)}^4$ (9)

其中 $u_n={\left[u_n\left(x_0\right),u_n\left(x_1\right),\cdots ,u_n\left(x_M\right)\right]}^{\text{T}}$ ， $X={\left[x_0,x_1,\cdots ,x_M\right]}^{\text{T}}$ ，并且上式中“/”，“ ${\left(1+X\right)}^4$ ”为矩阵按照

元素进行除法、乘方运算。

代入边界条件，利用MATLAB可以计算出方程的近似解，如图1：

Figure 1. The interpolation effect when M = 10

图1. 当M = 10时的插值效果

由图1可以看出，以M = 10为例，即选取11个插值节点，插值效果已经比较好，下面做出误差数量级随插值节点个数变化的图像，如图2：

Figure 2. The magnitude of error changes with the number of interpolated nodes

图2. 误差数量级随插值节点个数变化走势图

2.4. 结论

对于高阶非线性常微分方程的数值求解，本文采用等距节点的Lagrange插值方法进行逼近求解，并通过具体的数值结果表明，本方法获得的数值解可以有效地逼近精确解。本方法优点在于普适性好，可以用于其他阶数以及各种非线性形式的方程的求解。此外，该方法算法结构简单，易于编程计算，求解时间短，本方法对求解高阶非线性常微分方程具有重要意义。

本方法的不足之处在于，对于不同的方程，由于每个方程都具有独特的算法格式，不便使用于求解方程数目多，种类复杂的问题，本文后续也会尝试推导更加适用于大规模其求解问题的方法。

基金项目

河南省自然科学基金(No. 202300410156)；

河南科技大学研究训练计划(SRTP)项目(批准号：2022230)。