一维Sobolev方程的重心插值配点法

Barycentric Interpolation Collocation Method for One-Dimensional Sobolev Equation

摘要: 本文使用重心Lagrange插值配点法求解一维Sobolev方程的数值解，分别采用等距节点和Chebyshev节点进行数值计算。实验结果表明：在使用重心Lagrange插值求解一维Sobolev方程的数值解时，采用第二类Chebyshev节点可取得更高精度的数值解。

Abstract: This paper uses the center of gravity Lagrange interpolation method to solve the numerical solution of the one-dimensional Sobolev equation. Equidistant nodes and Chebyshev nodes are used for numerical calculation. The experimental results show that: when using the center of gravity Lagrange interpolation to solve the numerical solution of the one-dimensional Sobolev equation, using the second type of Chebyshev node can obtain a higher precision numerical solution.

1. 引言

本文讨论如下Sobolev方程的数值解：

$\{\begin{array}{l}{u}_t-u_{xx}=f\left(x,t\right),\left(x,t\right)\in \left(0,1\right)\times \left(0,1\right]\\ u\left(x,0\right)=0,x\in \left(0,1\right)\\ u\left(1,t\right)=u\left(0,t\right)=0,t\in \left[0,1\right]\end{array}$ (1)

其中 $\epsilon$ 和 $\gamma$ 为两个正常数， $f\left(x,t\right)$， $\phi \left(x,t\right)$ 和 $u_0\left(x\right)$ 为三个已知函数。

Sobolev方程在许多实际工程领域中都有着广泛的应用 [1] [2] [3]。近几十年来，Sobolev方程一直受到密切的关注，目前已经有许多关于Sobolev方程的数值研究 [4] [5] [6] [7]。 1990年，Yanping Lin研究了二维齐次Dirichlet边界条件下非线性Sobolev方程解的Galerkin逼近，得到了连续Crank-Nicolson和外推Crank-Nicolson逼近的误差估计 [4]。郭会和芮洪兴于2006年提出将最小二乘Galerkin有限元法应用于Sobolev方程，并通过误差估计表明该方法在L2范数意义下具有最优收敛阶 [5]。2013年，李宏和周文文从时间二阶精度的Crank-Nicolson时间半离散格式出发，构造出了Crank-Nicolson全离散化的有限元格式，并导出了误差估计 [6]。2017年，经鑫在文献 [7] 中建立了Sobolev方程的一个两层隐格式，证明了Sobolev方程在该格式下解的存在唯一性，并通过数值实验进行验证。以上介绍的这些方法均为有限差分法和有限元素法。但在实际求解过程中，这些方法不易实现，数值模拟也比较困难，尤其在使用有限元求解偏微分方程时，很难构造出满足特定条件的有限元格式。而且当方程维数为二维甚至三维的时候，数值模拟实验也比较复杂。本文我们所要求解的Sobolev方程是与时间t有关的一类方程，使用差分法和有限元法更会增加求解的难度。

W. J. Taylor于1945年发现了计算插值多项式的重心公式 [8]。2004年，J. P. Berrut和L. N. Trefethen等人将Lagrange公式改为重心公式的形式，得到了重心Lagrange插值公式 [9]。2007年，Michael S. Floater和Kai Hormann提出了重心有理插值 [10]，该方法操作简单，已经应用于许多偏微分方程的求解，都得到了精度较高的数值解 [11] [12] [13] [14]。重心插值配点法在求解二维或者三维偏微分方程的数值解时，数值模拟难度也不是很大。因此，本文将分别采用重心Lagrange插值配点法和重心有理插值配点法求解一维Sobolev方程，并对这两种插值法所得数值结果进行比较与分析。

2. 重心Lagrange插值与重心有理插值

给定 $n+1$ 个不同插值节点 $x_i\left(i=0,\cdots ,n\right)$ 及其相应的一组实数 $f_i$，那么Lagrange插值多项式公式可写为

$p\left(x\right)={\displaystyle \underset{i=0}{\overset{n}{\sum }}{l}_i\left(x\right)f_i}$ (2)

其中

$l_i\left(x\right)=\frac{{\displaystyle \underset{k=0,k\ne i}{\overset{n}{\prod }}\left(x-x_k\right)}}{{\displaystyle \underset{k=0,k\ne i}{\overset{n}{\prod }}\left(x_i-x_k\right)}}$ (3)

定义重心插值权：

${\omega }_i=\frac{1}{{\displaystyle \underset{i\ne k}{\prod }\left(x_i-x_k\right)}},i=0,1,\cdots ,n$ (4)

则可将(2)式重写为

$p\left(x\right)={\displaystyle \underset{j=0}{\overset{n}{\sum }}{\xi }_j\left(x\right)f_j},j=0,1,2,\cdots ,n$ (5)

其中

${\xi }_j\left(x\right)=\frac{\frac{{\omega }_j}{x-x_j}}{{\displaystyle \underset{j=0}{\overset{n}{\sum }}\frac{{\omega }_j}{x-x_j}}}$ (6)

公式(5)即为重心Lagrange插值公式。

重心有理插值是用有理函数进行插值。我们需要考虑，对于给定的节点分布及其相应的函数值，是否存在满足插值条件的有理函数。将式(5)重新定义为：

${\omega }_k={\displaystyle \underset{i\in J_k}{\sum }{\left(-1\right)}^i}{\displaystyle \underset{j=i,j\ne k}{\overset{i+d}{\prod }}\frac{1}{x_k-x_j}},k=0,1,2,\cdots ,n$ (7)

此处d为一给定的正整数，且 $0\le d\le n$。 $I=\left\{0,1,2，\cdots ,n\right\}$， $J_k=\left\{i\in I:k-d\le i\le k\right\}$。

从而得到重心有理插值公式：

$r\left(x\right)=\frac{{\displaystyle \underset{j=0}{\overset{n}{\sum }}\frac{{\omega }_j}{x-x_j}{f}_j}}{{\displaystyle \underset{j=0}{\overset{n}{\sum }}\frac{{\omega }_j}{x-x_j}}}$ (8)

在使用重心有理插值法求解偏微分方程时，参数d的选取非常重要，它会直接影响计算精度。

3. 一维Sobolev方程的计算

对方程(1)，采用第二类Chebyshev节点，将空间方向和时间方向分别离散为m，n个插值节点 $a=x_1<\cdots <x_m=b$， $0=t_1<\cdots <t_n=T$，共有 $m\times n$ 个张量型插值节点 $u\left(x_i,t_j\right)$，其中 $i=1,2,\cdots ,m$， $j=1,2,\cdots ,n$。则

$u\left(x,t\right)={\displaystyle \underset{i=1}{\overset{m}{\sum }}{\displaystyle \underset{j=1}{\overset{n}{\sum }}{L}_i\left(x\right)T_j\left(t\right)u_{ij}}}$ (9)

其中 $u\left(x_i,t_j\right)=u_{ij}$。 $L_i\left(x\right),T_j\left(t\right)$ 分别为x方向和t方向的插值基函数。

将上式代入Sobolev方程(1)得到

${\displaystyle \underset{i=1}{\overset{m}{\sum }}{\displaystyle \underset{j=1}{\overset{n}{\sum }}\left(L_i\left(x\right){T^{\prime }}_j\left(t\right)-\epsilon {L^{″}}_i\left(x\right){T^{\prime }}_j\left(t\right)-\gamma {L^{″}}_i\left(x\right)T_j\left(t\right)\right)u_{ij}}}=f\left(x,t\right)$ (10)

记 $f\left(x_i,t_j\right)=f_{ij}$， $L_j\left(x_i\right)={\delta }_{ij},T_k\left(t_j\right)={\delta }_{jk}$，用 $\otimes$ 表示矩阵的Kronecker积，则可将(10)式写为

$\left(I_m\otimes T^{\left(1\right)}-\epsilon \left(L^{\left(2\right)}\otimes T^{\left(1\right)}\right)-\gamma \left(I_m\otimes T^{\left(1\right)}\right)\right)U=F$ (11)

其中 $L^{\left(r\right)}$ 为插值节点 $x_1,x_2,\cdots ,x_m$ 的r阶微分矩阵， $T^{\left(r\right)}$ 为插值节点 $t_1,t_2,\cdots ,t_l$ 的r阶微分矩阵。且式(11)中的U和F分别为如下形式

$U=\left[u_{11},u_{12},\cdots ,u_{1n},u_{21},u{}_{22},\cdots ,u_{2n},\cdots ,u_{m1},u_{m2},\cdots ,u_{mn}\right]$

$F=\left[f_{11},f_{12},\cdots ,f_{1n},f_{21},f_{22},\cdots ,f_{2n},\cdots ,f_{m1},f_{m2},\cdots ,f_{mn}\right]$

令 $L=I_m\otimes T^{\left(1\right)}-\epsilon \left(L^{\left(2\right)}\otimes T^{\left(1\right)}\right)-\gamma \left(I_m\otimes T^{\left(1\right)}\right)$，则方程(21)可以写为

$LU=F$ (12)

4. 数值算例

算例4.1 方程(1)的精确解为

$u\left(x,t\right)=\left(0.5x^2-0.5x\right)t$

下面分别用这两种插值法求解方程(1)的数值解，节点类型选取第二类Chebyshev节点和等距节点，将 $x,t$ 离散化， $0=x_1<\cdots <x_m=1$， $0=t_1<\cdots <t_n=1$，其中m，n为剖分节点数。通过数值模拟，所得绝对误差如表1和表2所示。

Table 1. Absolute error of Lagrange interpolation

表1. 重心Lagrange插值的绝对误差

由表1的结果可以看出，在计算一维Sobolev方程的数值解时，当所选节点数目较小时，采用等距节点逼近效果较好，但随着节点数目增加，会出现“龙格现象”。整体来看，采用第二类Chebyshev节点所得数值解逼近效果更好。因此，在使用重心Lagrange插值配点法求解偏微分方程数值解时，为了取得较高精度的数值解，通常会选取第二类Chebyshev节点进行数值计算。

由表1结果可知，选取等距节点产生的结果并不理想。并且通过参考大量文献，可以发现选取等距节点所得实验结果往往都是病态的，所以在使用重心有理插值求解Sobolev方程时，我们将只选取第二类Chebyshev节点，其所得结果见表2。由于参数d的选择对结果的精度影响很大，所以需谨慎选择参数d。此处我们令d = 5。

Table 2. Absolute error of rational interpolation

表2. 重心有理插值的绝对误差

表2结果表明，在求解Sobolev方程时，重心有理插值也有很好的逼近效果。将表2结果与表1的第二列比较，可以发现，这两种插值方法在选取第二类Chebyshev节点时都可产生较高精度的结果，且重心有理插值产生的数值解精度要略高一些。由表1和表2结果还可以发现，不仅节点类型对数值结果有影响，节点数目的大小对计算结果也有较大的影响。随着节点数目的增大，数值精度反而在逐渐下降，最终趋于稳定。由表1第三列可以看出，节点数目的大小对结果的影响非常大，从整体来看，选取等距节点产生的结果是很不稳定的。因此，具体的方程数值解求解过程中，一般不选取等距节点。

5. 结论

采用这两种插值方法求解一维Sobolev方程时，为了获得逼近效果好的数值解，通常建议选取第二类Chebyshev节点。本文中有理插值所产生的数值结果精度略高，但并不能因此说明重心有理插值一定能够获得比重心Lagrange插值精度高的数值解。可将这两种插值法推广到二维和三维Sobolev方程，甚至其他偏微分和常微分方程的数值求解，并分析比较这两种方法各自的优势。