半直线上Kortewego-de Vries方程的Laguerre谱配置方法

Generalized Laguerre Spectral-Collocation Method for KdV Equations on the Half Line

摘要: 以Laguerre-Gauss-Radau节点为配置点，用带松弛因子的Lagrange插值函数逼近半直线上的Kortewego-de Vries方程初边值问题的理论解，给出算法格式和相应的数值结果，表明所提算法格式的有效性和高精度。对理论解中参数的不同取值，通过适当地选择插值函数中的松弛因子，数值解可以很好地匹配理论解，而且所给算法对长时间的计算仍然有效。

Abstract: Interpolation function approximations with relaxation factor by using Laguerre-Gauss-Radau nodes as collocation points to the Korteweg-de Vries equation on semi-infinite intervals are considered. The validity and high accuracy of the proposed algorithm are demonstrated. By choosing the relaxation factor of the interpolation function properly, the numerical solution can match the theoretical solution well, and the algorithm is still valid for a long time.

1. 引言

记 $R^{+}=\left\{x|0<x<\infty \right\}$ 为简单，记 $\frac{\partial u}{\partial z}={\partial }_{z}u$。下式为半直线上Kortewego-de Vries (KdV)方程初边值

问题 [1] [2] [3]：

$\{\begin{array}{l}{\partial }_{t}u+\gamma u{\partial }_{x}u+\nu {\partial }_x^{3}u=0,x\in R^{+},t>0,\\ u\left(x,t\right),{\partial }_{x}u\left(x,t\right)\to 0,x\to +\infty ,t\ge 0,\\ u\left(x,0\right)=g_0\left(x\right),x\in R^{+},\\ u\left(0,t\right)=g_1\left(t\right),t\ge 0.\end{array}$ (1)

其中 $\gamma$ 和 $\nu$ 是常数。文献 [4] [5] [6] 用谱方法或有限元/b-样条有限元方法研究有界区域上KdV方程的数值解，针对方程(1)中的参数 $\gamma ,\nu$ 的不同取值，文献 [7] [8] 利用空间Hermite谱配置方法、时间有限差分方法或Chebyshev-Hermite多项式时空谱配置方法求解其Cauchy问题的数值解；文献 [9] 用Laguerre拟谱方法研究了半直线上非线性热传导方程的数值解。针对问题(1)的孤波解的性态，用含有因子 ${\text{e}}^{-\beta x/2}\left(\beta >0\right)$ 的插值函数可以更好地逼近问题的理论解，同时通过适当选取伸缩因子 $\beta$ 可以改进数值误差精度。另外，用Gauss型节点得到的Lagrang插值多项式，与其相关的高阶微分矩阵是一阶微分矩阵的乘积，这为实际计算带来极大的方便。

2. 基于Gauss型节点的Lagrange插值函数及其微分矩阵

次数为l的广义拉盖尔多项式定义为 [9]：

$L_l^{\left(\alpha ,\beta \right)}\left(x\right)=\frac{1}{l!}{x}^{-\alpha }{\text{e}}^{\beta x}{\partial }_x^l\left(x^{l+\alpha }{\text{e}}^{-\beta x}\right),l=0,1,\cdots .$ (2)

广义拉盖尔函数定义为：

${\tilde{L}}_l^{\left(\alpha ,\beta \right)}\left(x\right)={\text{e}}^{-\frac{1}{2}\beta x}{L}_l^{\left(\alpha ,\beta \right)}\left(x\right),l=0,1,\cdots .$ (3)

特别， $L_l^{\left(\beta \right)}\left(x\right)=L_l^{\left(0,\beta \right)}\left(x\right),{\tilde{L}}_l^{\left(\beta \right)}\left(x\right)={\text{e}}^{-\frac{1}{2}\beta x}{L}_l^{\left(\beta \right)}\left(x\right)$。

令 $x_j\left(0\le j\le N\right)$ 是 $x{\partial }_x{L}_{N+1}^{\left(\beta \right)}\left(x\right)=0$ 的根。以 $x_j$ 为节点的通常的Lagrange插值基函数为 [9] [10] [11]

${\phi }_j^{\left(\beta \right)}\left(x\right)=\frac{{\text{e}}^{-\beta x/2}}{\beta \left(N+1\right){\text{e}}^{-\beta x_j/2}}\frac{x{\partial }_x{L}_{N+1}^{\left(\beta \right)}\left(x\right)}{\left(x_j-x\right)L_{N+1}^{\left(\beta \right)}\left(x_j\right)},j=0,1,\cdots ,N$

对任意 $u\left(x\right)\in C\left(R^{+}\right)$，其Lagrange插值多项式为：

$p_N^{\left(\beta \right)}\left(x\right)={\displaystyle \underset{j=0}{\overset{N}{\sum }}u\left(x_j\right){\phi }_j^{\left(\beta \right)}\left(x\right)},x\in R^{+}.$

对 $p_N\left(x\right)$ 求m阶导数得，

${\partial }_x^m{p}_N^{\left(\beta \right)}\left(x\right)={\displaystyle \underset{j=0}{\overset{N}{\sum }}u\left(x_j\right){\partial }_x^m{\phi }_j^{\left(\beta \right)}}\left(x\right).$

令

$D^{\left(m\right)}={\left({\partial }_x^m{\phi }_j^{\left(\beta \right)}\left(x_k\right)\right)}_{0\le k,j\le N},D:=D^{\left(1\right)}=\left(d_{k,j}^{\left(\beta \right)}={\partial }_x{\phi }_j^{\left(\beta \right)}\left(x_k\right)\right),$

则由文献 [9] [10] [11]

$d_{k,j}^{\left(\beta \right)}=\{\begin{array}{l}\frac{{\tilde{L}}_{N+1}^{\left(\beta \right)}\left(x_k\right)}{\left(x_k-x_j\right){\tilde{L}}_{N+1}^{\left(\beta \right)}\left(x_j\right)},k\ne j,\\ -\frac{\left(N+1\right)\beta }{2},k=j=0,\\ 0,k=j\ne 0.\end{array}$ (4)

$D^{\left(m\right)}=D^m.$ (5)

特别，

$D^{\left(3\right)}={\left({\tilde{d}}_{k,j}^{\left(\beta \right)}={\partial }_x^3{\phi }_j^{\left(\beta \right)}\left(x_k\right)\right)}_{0\le k,j\le N}.$

3. KdV方程的谱配置方法

3.1. KdV方程的谱配置格式

在式(1)中去 $\gamma =\nu =1$，考虑KdV方程的初边值问题：

$\{\begin{array}{l}{\partial }_{t}u+u{\partial }_{x}u+{\partial }_x^{3}u=0,x\in R^{+},t>0,\\ u\left(x,t\right),{\partial }_{x}u\left(x,t\right)\to 0,x\to +\infty ,t\ge 0,\\ u\left(x,0\right)=g_0\left(x\right),x\in R^{+},\\ u\left(0,t\right)=g_1\left(t\right),t\ge 0.\end{array}$ (6)

用 $w_N\left(x,t\right)={\displaystyle \underset{j=0}{\overset{N}{\sum }}u\left(x_j,t\right){\phi }_j^{\left(\beta \right)}\left(x\right),x\in R^{+}}$，逼近式(6)的解 $u\left(x,t\right)$，将其代入式(6)，得

$\{\begin{array}{l}{\partial }_t{w}_N\left(x,t\right)+w_N\left(x,t\right){\partial }_x{w}_N\left(x,t\right)+{\partial }_x^3{w}_N\left(x,t\right)=0,\\ x=x_k,k=1,2,\cdots ,N,t>0\\ w_N\left(x,0\right)=g_0\left(x\right),x=x_k,k=0,1,\cdots ,N,\\ w_N\left(0,t\right)=g_1\left(t\right),t\ge 0\end{array}$ (7)

式(7)等价于

$\{\begin{array}{l}{\partial }_t{u}_k\left(t\right)=-u_k\left(t\right){\displaystyle \underset{j=0}{\overset{N}{\sum }}{u}_j}\left(t\right)d_{k,j}^{\left(\beta \right)}\\ -{\displaystyle \underset{j=0}{\overset{N}{\sum }}{u}_j}\left(t\right){\tilde{d}}_{k,j}^{\left(\beta \right)},u_j\left(t\right)=u\left(x_j,t\right),k=1,2,\cdots ,N,\\ u_k\left(0\right)=g_0\left(x_k\right),k=0,1,\cdots ,N,\\ u_0\left(t\right)=g_1\left(t\right).\end{array}$ (8)

令 $X\left(t\right)={\left(u_1\left(t\right),u_2\left(t\right),\cdots ,u_N\left(t\right)\right)}^{\text{T}},X_0={\left(g_0\left(x_1\right),g_0\left(x_2\right),\cdots ,g_0\left(x_N\right)\right)}^{\text{T}}$，及

$A=\left(\begin{array}{cccc}{d}_{1,1}^{\left(\beta \right)}& d_{1,2}^{\left(\beta \right)}& \cdots & d_{1,N}^{\left(\beta \right)}\\ d_{2,1}^{\left(\beta \right)}& d_{2,2}^{\left(\beta \right)}& \cdots & d_{2,N}^{\left(\beta \right)}\\ ⋮& ⋮& \ddots & ⋮\\ d_{N,1}^{\left(\beta \right)}& d_{N,2}^{\left(\beta \right)}& \cdots & d_{N,N}^{\left(\beta \right)}\end{array}\right),B=\left(\begin{array}{cccc}{\tilde{d}}_{1,1}^{\left(\beta \right)}& {\tilde{d}}_{1,2}^{\left(\beta \right)}& \cdots & {\tilde{d}}_{1,N}^{\left(\beta \right)}\\ {\tilde{d}}_{2,1}^{\left(\beta \right)}& {\tilde{d}}_{2,2}^{\left(\beta \right)}& \cdots & {\tilde{d}}_{2,N}^{\left(\beta \right)}\\ ⋮& ⋮& \ddots & ⋮\\ {\tilde{d}}_{N,1}^{\left(\beta \right)}& {\tilde{d}}_{N,2}^{\left(\beta \right)}& \cdots & {\tilde{d}}_{N,N}^{\left(\beta \right)}\end{array}\right),C=\left(\begin{array}{ccc}{d}_{1,0}^{\left(\beta \right)}& \cdots & 0\\ 0& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & d_{N,0}^{\left(\beta \right)}\end{array}\right),F=\left(\begin{array}{c}{\tilde{d}}_{1,0}^{\left(\beta \right)}\\ {\tilde{d}}_{1,0}^{\left(\beta \right)}\\ ⋮\\ {\tilde{d}}_{N,0}^{(\; \beta \; )}\end{array}\right)$

式(8)的矩阵形式为

$\{\begin{array}{l}\frac{\text{d}X\left(t\right)}{\text{d}t}=-BX\left(t\right)-X\left(t\right)\cdot *\left(AX\left(t\right)\right)-g_1\left(t\right)CX\left(t\right)-g_1\left(t\right)F,t>0\\ X\left(0\right)=X_0,\end{array}$ (9)

其中符号“ $\cdot *$ ”表示对应位置元素相乘。

3.2. 数值结果

用格式(9)求解式(6)。在时间方向用步长为 $\tau$ 的Crank-Nicolson格式离散式(9)，得

$\{\begin{array}{l}\left(I+\frac{1}{2}\tau \left(B+g_1\left(t+\tau \right)C\right)\right)X\left(t+\tau \right)\\ =\left(I-\frac{1}{2}\tau \left(B+g_1\left(t\right)C\right)\right)X\left(t\right)-\frac{1}{2}\tau \left(X\left(t+\tau \right)\cdot *\left(DX\left(t+\tau \right)\right)+X\left(t\right)\cdot *\left(DX\left(t\right)\right)\right)\\ -\frac{1}{2}\tau \left(g_1\left(t\right)+g_1\left(t+\tau \right)\right),t=0,\tau ,2\tau ,\cdots T-\tau ,\\ X\left(0\right)=X_0.\end{array}$ (10)

上式中I是N阶单位矩阵。显然(10)式是关于 $X\left(t\right)$ 的非线性方程，通常用解非线性方程的迭代方法求其近似解。实际计算时一般用的是Newton迭代方法，但计算迭代矩阵比较麻烦。为简单起见，用如下的不动点迭代方法：对时间方向 $t+\tau$ 处构造迭代格式：

$\{\begin{array}{l}\left(I+\frac{1}{2}\tau \left(B+g_1\left(t+\tau \right)C\right)\right)X_{n+1}\left(t+\tau \right)\\ =\left(I-\frac{1}{2}\tau \left(B+g_1\left(t\right)C\right)\right)X_n\left(t\right)-\frac{1}{2}\tau \left(X_{n+1}\left(t+\tau \right)\cdot *\left(D{X}_{n+1}\left(t+\tau \right)\right)+X_n\left(t\right)\cdot *\left(D{X}_n\left(t\right)\right)\right)\\ -\frac{1}{2}\tau \left(g_1\left(t\right)+g_1\left(t+\tau \right)\right),t=0,\tau ,2\tau ,\cdots T-\tau ,\\ X\left(0\right)=X_0.\end{array}$ (11)

在迭代过程中给出条件：对任给的 $\epsilon >0$，当 $\left|X_{n+1}\left(t+\tau \right)-X_n\left(t+\tau \right)\right|\le \epsilon$ 时终止迭代，可得 $X\left(t\right)$ 在 $t+\tau$ 处的值。

在(1)式中取 $\gamma =\nu =1$，相应的KdV方程有如下孤波解 [1] [10]：

$u=12{\kappa }^2{\mathrm{sech}}^2\left(\kappa \left(x-4{\kappa }^{2}t-x_0\right)\right),$ (12)

上式中的 $\kappa$ 和 $x_0$ 是任意参数。用如下的无穷范数

$E_{N,\tau }\left(t\right)=\underset{0\le j\le N}{\max }\left|u\left(x_j,t\right)-w_N^{\left(\beta \right)}\left(x_j,t\right)\right|$

计算数值解与理论解之间的误差。图(a)给出 $t=1$ 时不同时间步长 $\tau$ 对应的误差 $E_{N,\tau }\left(t\right)$ 的常用对数 ${\log }_{10}{E}_{N,\tau }\left(t\right)$ 随N的变化关系。可以发现，数值误差 $E_{N,\tau }\left(1\right)$ 只须空间节点数 $N\ge 8$，时间步长 $\tau \le 0.001$ 时，误差即达到 $O\left({10}^{-14}\right)$ 量级，说明所提算法格式有谱精度；图(b)与图(a)相比较，就是理论解(12)中的参数 $\kappa$ 由0.3变为0.5，插值函数中的伸缩因子 $\beta$ 由1变为1.5。可以看出，算法格式(9)对(10)式中的参数 $\kappa$ 有较强的适应性，当 $\kappa$ 较大时插值基函数中的 $\beta$ 也适当的变大，空间方向仍可以达到谱精度；而已有的文献通常对较小的 $\kappa$ 逼近精确解(10)，所以算法格式(9)对参数 $\kappa$ 是稳定的，这是所提算法的一个优点。图(c)给出时间t分别取1和50时的误差变化，当 $t\ge 50$ 时，需要增加节点个数，所提算法格式仍然有效。图(d)是选取插值函数中的伸缩因子 $\beta$ 不同的值时的误差，表明对孤波解而言，当选择较大的 $\beta$ 时，数值误差 $E_{N,\tau }\left(1\right)$ 会变的小一些， $\beta$ 的值通过试验得到的，理论上还没有解决如何选取最好的 $\beta$ 值。

(a) $\beta =1,\kappa =0.3,x_0=-20$ (b) $\beta =1.5,\kappa =0.5,x_0=-20$

(c) $\beta =1,\kappa =0.3,x_0=-20,\tau =0.001$ (d) $\kappa =0.3,x_0=-20,\tau =0.01$

4. 结论

为了避免通常的等距节点Lagrange插值多项式当次数较高时( $N\ge 10$ )出现Runge现象，而采用Gauss型节点构造Lagrange插值多项式逼近问题(1)解。考虑到问题(1)的孤波解的性态，当 $\left|x\right|\to +\infty$ 时，其迅速衰减为零，用含有因子 ${\text{e}}^{-\beta x/2}$ 的插值函数可以使逼近函数与问题(1)的理论解更好的吻合，从而得到较高精度的数值解。特别，所给算法也适合于(10)中较大的参数 $\kappa$ 及变量时间t值，表明算法格式(9)是强健的。

基金项目

国家自然科学基金项目(11371123)；SRTP (N.2019198)。