全直线区域上的对角化Chebyshev有理谱方法

赵云阁，余旭洪

（上海理工大学 理学院，上海 200093）

科学与工程中的许多问题都被归结为无界区域问题，例如，无界区域的流体问题、量子力学中的非线性波动方程等。随着科技的发展，无界区域问题越来越受到人们的重视。因此，如何高效地求解此类问题一直是研究的热点，同时由于定义区域的无界特性也使得此类问题成为数值分析和实际应用的难点。在实践中人们采用的第一类方法是人工边界方法，即将无界区域问题转化为有界区域问题，但该类方法往往会带来人工误差，甚至导致数值不稳定；第二类方法是直接使用Hermite或Laguerre函数逼近无界问题[1-3]，但该类方法往往由于逼近函数的Gauss点过于集中在原点附近，而使得远离原点区域的逼近效果不够理想；第三类方法是通过变量变换将无界区域问题变换成有界区域上的奇异问题，再利用Jacobi谱方法等进行数值求解[4]，但该类方法所得奇异问题一般形式过于复杂，不利于理论分析与数值计算；第四类方法即有理谱方法[5-7]，其作为一类求解无界区域问题的行之有效的方法已越来越被人们重视。

通常应用Chebyshev有理谱方法求解无界区域问题时，虽然问题离散后所得代数方程组的系数矩阵往往是带状稀疏矩阵，但条件数却往往以基函数个数的幂次方快速增长，从而导致病态矩阵问题，给实际计算带来了困难。因此，人们希望获得一类对角化的系数矩阵来解决条件数的快速增长问题。受文献[8-10]的启发，本文将建立全直线区域上二阶微分方程的对角化Chebyshev有理谱方法，并构造相应的Fourier型基函数。特别地，通过Fourier型基函数，能将原问题真解展开成Fourier级数的形式，并将数值解表示成级数的截断形式。在数值实验中可以看出，新算法优化了计算过程，减少了计算量，并且简单易行。

1 预备知识及引理

对任意的正整数r，在区间 Λ =(−∞,+∞)上按照通常方式定义Sobolev空间Hr(Λ),其內积和范数定义为

当r=0时，记且和为了叙述方便，分别记和

记Tk(y)为定义在上由如下三递推公式生成的k次Chebyshev多项式：

Chebyshev多项式具有正交性：

其中，δk,l为 Kronecker 数，c0=2且当k≥1时，ck=1。此外，Chebyshev多项式还满足递推关系：

且有Tk(±1)=(±1)k和Tk′(±1)=(±1)k−1k2。

定义全直线区域 Λ上的k次Chebyshev有理函数

由式（5）和式（2）可知Chebyshev有理函数具有正交性：

且有

由式（1）,(3）和（4）可得递推关系：

引理1任意的有如下等式成立：

证明由式（8），（9）和（10）可以得到

且

结合式（14）和式（15），很容易得到结论式（12）成立。

2 对角化Chebyshev有理谱方法

建立全直线区域上带渐近边界条件的二阶椭圆型微分方程及非线性热传导方程的对角化Chebyshev有理谱方法，并构造相应的Fourier型Sobolev正交基函数，使得方程的真解和数值解能展开成Fourier级数的形式。

2.1 二阶椭圆型边值问题的对角化有理谱方法

考虑定义在全直线区域上的二阶模型问题

为了获得谱格式式（18）的对角化逼近格式，需要构造在Sobolev内积 Aλ(·,·)下的正交基函数

引理 2设是满足

的Sobolev正交函数，则有

证明首先利用数学归纳法证明式（20）成立。由Schmidt正交化思想，令

因此，由式（12）和式（6）可得

成立。事实上，由式（17），（23），（12），（6）以及归纳法假设，当且时，有

从而由式（22）可得，当k≥5时，

另一方面，由式（19）和式（20）可得

因此，有

这样即得结论式（21）成立。

定理1设u(x)和分别是方程（16）和有理谱格式式（18）的真解和数值解，那么，和分别可以展开成基函数的级数形式和级数截断形式，即

2.2 非线性热传导方程边值问题的对角化有理谱方法

考虑非线性热传导方程

这里，µ >0且f(x,t)是已知函数。

建立方程（24）的混合差分−有理谱格式，其中，时间方向使用中心差分格式，空间方向使用对角化Chebyshev有理谱格式。

其中

其中

为了获得Chebyshev有理谱格式式（27）的对角化逼近格式，需要构造在Sobolev内积下的正交基函数采用与引理2同样的证明方法，可得引理3。

引理 3设ψk∈Pk(Λ)是满足ψk−Rk∈ Pk−1(Λ)和的Sobolev正交函数，则有

定理2设是有理谱格式式（27）的解，那么，有

有理谱格式（27）是隐式格式，在实际计算中需要应用迭代法计算解的系数。

3 数值模拟

应用对角化Chebyshev有理谱方法求解二阶微分方程（16）和非线性热传导方程（24），并验证本文所提方法的快速与高精度特性。

首先，考察二阶方程（16）在λ =1 时，针对测试函数具有不同衰减速度时的收敛性。

图1 指数衰减时谱格式（18）的误差Fig.1 Errors of the scheme (18) with exponential decay function

针对震荡且无穷远处代数衰减的测试函数u(x)=(x2+1)−2sin(2x),图 2 绘 制了谱格式式（18）取自然对数后的离散L2和H1误差与 lg10N之间的关系。从图中可以发现，谱格式式（18）对衰减速度较慢的测试函数仍然具有代数收敛率。

图2 代数衰减时谱格式（18）的误差Fig.2 Errors of the scheme (18) with algebraic decay function

其次，考察非线性热传导方程问题（24）针对震荡且无穷远处指数衰减的测试函数u(x,t)=e−x2sin(2x+t)且 µ =1 时的收敛性。图3绘制了有限差分−有理谱格式式（27）在τ=0.1,0.01,0.001,0.000 1时的离散L2误差与N之间关系。从图3中可以发现，谱格式式（27）保持了几何收敛率；且随着时间步长τ的 减小，数值误差也越来越小。图4绘制了在时间 0 ⩽t⩽100取τ=0.001 时的长时间计算的L2误差，可以发现，本文所建立的对角化算法具有长时间计算稳定性。

图3 谱格式（27）的 L 2误差Fig.3 L 2 errors of the scheme (27)

最后，为了体现对角化Chebyshev有理谱方法的最主要优点，比较新方法与经典Chebyshev有理谱方法以及Hermite谱方法的刚度矩阵条件数以及计算速度的差别。

在表2和表3中分别列出了应用对角化Chebyshev有理谱方法和经典Chebyshev有理谱方法求解方程（16）时的L∞误差以及CPU耗时（单位：s），从表中可以看出，对角化方法明显地节约了计算时间。

表1 经典Chebyshev有理谱方法(CRSM)和Hermite谱方法(HSM)的条件数Tab.1 Condition numbers of the classical Chebyshev and Hermite spectral methods

表2 方程（16）的对角化Chebyshev有理谱方法耗时Tab.2 Time consumption of the diagonalized Chebyshev rational spectral method for (16)

表3 方程（16）的经典Chebyshev有理谱方法耗时Tab.3 Time consumption of the classical Chebyshev rational spectral method for (16)