圆域上变系数二阶椭圆方程有效的谱方法及其在奇异非线性问题中的应用

刘忠敏,安 静,陈 悦

(贵州师范大学 数学科学学院，贵州 贵阳550025)

二阶椭圆方程问题是一类重要而又基础的问题，在工程和科学领域中很多复杂的非线性问题最终都归结为求解二阶椭圆方程，如非线性Schrödinger方程[1-3]、Allen-Cahn方程[4]和Navier-stokes方程[5-7]。关于二阶椭圆方程的理论分析和数值计算已有很多成果[8-13]。这些数值方法主要基于低阶的有限元法，对于一些正则性较低的拟线性二阶椭圆方程，有限元法确实是一种较好的数值方法，尤其是与自适应网格相结合的有限元法。然而，对于正则性较高的二阶椭圆方程，提出一些基于高阶多项式逼近的有限元法也是有意义的，特别是对于一些带有曲边或曲面边界的特殊区域(如圆域、球域等)上的二阶变系数和非线性的椭圆方程。我们在用有限元求解这些特殊区域上的二阶椭圆问题时,通常以直代曲去逼近这些边界，从而会引入一些额外的误差。谱方法是一种具有谱精度的高阶数值方法[14-16]。目前还未有将谱方法应用于求解圆域上二阶非线性椭圆方程的研究。

因此，本文提出了圆域上二阶变系数椭圆方程的一种有效的谱方法。通过利用极坐标变换，将原问题转化为极坐标下的一种等价形式。根据极条件、边界条件以及方向的周期性，引入了适当的Sobolev空间，建立了极坐标系下二阶变系数椭圆方程的一种弱形式及其离散格式。然后，利用Lax-Milgram引理证明了弱解的存在唯一性。再由非一致带权Sobolev空间中投影算子的逼近性质和傅里叶基函数的逼近性质，证明了逼近解的误差估计。此外,还将提出的算法延伸到奇异非线性二阶椭圆方程的计算，并给出了数值算例，数值结果表明本文算法是收敛的和高精度的。

1 极坐标下的弱形式及其离散格式

考虑如下二阶变系数椭圆方程：

(1)

其中：Ω={(x,y):x2+y2≤R2};α(x,y)为非负有界函数;∂Ω为Ω的边界。

Δu(x,y)=Lσ(t,θ)=

(2)

由于极坐标变换引入了极点奇性，为了(2)式的适定性，σ(t,θ)需要满足以下本质极条件：

∂θσ(-1,θ)=0。

令f(t,θ)=g(x,y),β(t,θ)=α(x,y),则(1)式在极坐标下的等价形式为

(3)

令ω=t+1表示一个权函数。为了推导(3)式的弱形式，需要引入一些如下的带权Sobolev空间：

相应的内积和范数分别为

和

|∂θσ|2dtdθ<∞,

σ(1,θ)=∂θσ(-1,θ)=0,

其内积和范数分别为

A(σ,δ)=F(δ),

(4)

其中

定义有限元空间

χMN=span{σmN(t)eimθ:σmN(t)∈QmN,|m|≤M},

其中QmN={q∈pN:mq(-1)=q(1)=0},pN为N次多项式空间。则(4)式的离散格式为：取σMN∈χMN,使得

A(σMN,δMN)=F(δMN),∀δMN∈χMN。

(5)

2 弱解的存在唯一性及逼近解的误差估计

我们先证明弱解的存在唯一性。为了叙述方便，用a≼b表示存在一个正的常数c,使得a≤cb。

|A(σ,δ)|≼‖σ‖1,*‖δ‖1,*,

证明由边界条件和Cauchy-Schwarz不等式可得

则有

则由β(t,θ)的有界性和Cauchy-Schwarz不等式有

‖σ‖1,*‖δ‖1,*。

另一方面，由β(t,θ)的非负性有

A(σ,σ)=

证毕。

‖δ‖1,*,

引理2设σ和σMN分别是(4)和(5)式的解，则有

证明由(4)和(5)式有

(6)

A(σMN,δMN)=F(δMN),∀δMN∈χMN。

(7)

在(6)式中取δ=δMN可得

A(σ,δMN)=F(δMN)。

(8)

将(8)和(7)式作差，由A(σ,δ)的双线性性质可得

A(σ-σMN,δMN)=0。

(9)

由引理1有

A(σ-σMN,σ-σMN-δMN+δMN)=

A(σ-σMN,σ-δMN)≼

‖σ-σMN‖1,*‖σ-δMN‖1,*。

由δMN的任意性有

证毕。

ρ(x+2π)=ρ(x)},

其内积和范数分别为

‖ρ(θ)-ρM(θ)‖k≼Mk-s|ρ(θ)|s。

其内积和范数分别为

由于σ(t,θ)在θ方向上以2π为周期，则由傅里叶基函数展开和截断可得

‖σ-σMN‖1,*≼M1-s+N1-q。

证明由于

(1+t)|∂t(σM-δMN)|2dtdθ+

则有

从而有

由引理3可得

由引理4可得

由引理2可得

M1-s+N1-q。

证毕。

3 算法的有效实现

我们对(5)式的算法进行详细描述。令

Ψmi(t)=Li(t)-Li+2(t),(i=0,1，…,N-2),

其中Li(t)为i次Legendre多项式。则逼近空间χMN可表示为

χMN=span{Ψmi(t)eimθ:-M≤m≤M,

i=0,1,…,N-1-sign(|m|)}。

令

eimθe-inθdtdθ,

我们将寻找

(10)

将表达式(10)代入(5)式，让δMN取遍χMN中的所有基函数，则离散格式(5)可化为以下线性系统：

(A+B+D)H=f。

其中：

H=[h-M,0,h-M,1,…,h-M,N-2,…,h00,h01,…,

h0,N-1,…,hM,0,hM,1,…,hM,N-2]T,

A=(akjnm),B=(bkjnm),

D=(dkjnm),f=(fkn)。

由Fourier基函数和Legendre多项式的正交性质可知，矩阵A、B都是稀疏的分块带状矩阵，矩阵D的稀疏性依赖于变系数α(x,y)的性质。

4 奇异非线性二阶椭圆方程

在前文提出的算法基础上，给出非线性二阶椭圆方程问题的算法。我们考虑如下非线性二阶椭圆方程问题：

(11)

类似(3)式的推导，(11)式在极坐标下的等价形式为

(12)

(13)

其中

则(13)式相应的离散格式为：取σMN∈χMN，使得

A(σMN,δMN)=F(δMN),∀δMN∈χMN。

(14)

为了求解非线性方程(14)，建立如下迭代格式：

因此，我们把非线性格式(14)转化成变系数的迭代格式，从而可用第3节提出的算法有效地求解。

5 数值算例

为了表明算法的收敛性和高精度，在MATLAB 2016b平台上进行一系列数值实验。设uMN(x,y)为精确解u(x,y)的逼近解，则由极坐标变换和变量替换有

定义逼近解uMN(x,y)与精确解u(x,y)间的误差

e(u(x,y),uMN(x,y))=

‖u(x,y)-uMN(x,y)‖L∞(Ω)=

‖σ(t,θ)-σMN(t,θ)‖L∞(D)。

例1我们考虑二阶变系数椭圆方程(1)，取α(x,y)=ex+y+1，u(x,y)=(x2+y2-1)cosxy，Ω={(x,y):x2+y2≤1}，显然u(x,y)满足方程(1)的边界条件，再将精确解u(x,y)代入方程(1)可得到方程(1)右端的g(x,y)。我们用前文提出的算法求解问题(1)。对于不同的N和M，表1列出了逼近解与精确解之间的误差。为进一步表明本文提出算法的收敛性和高精度，图1和图2给出了精确解与N=25、M=16时的逼近解对比图与不同的N和M时逼近解与精确解间的误差图像。

表1 对于不同的M和N，逼近解与精确解间的误差e(u(x,y),uMN(x,y))

图1 精确解和N=25、M=16时的逼近解对比图

图2 当N=30、M=20和N=50、M=25时的逼近解与精确解之间的误差图像

从表1可以观察到，当N≥15且M≥16时，逼近解达到了10-14左右精度。另外，从图1和图2也可观察到本文算法是收敛的和高精度的。

例2考虑奇异非线性二阶椭圆方程(11)。类似地，取Ω={(x,y):x2+y2≤1}，α(x,y)=ex+y+1,u(x,y)=(x2+y2-1)ex+y,将u(x,y)代入方程(11)可得到右端的g(x,y)。对于不同的N和M，表2列出了逼近解与精确解间的误差。图3和图4分别给出了精确解与N=35、M=15时逼近解的对比图和不同的N和M时逼近解与精确解间的误差图像。

表2 对于不同的M和N，逼近解与精确解间的误差e(u(x,y),uMN(x,y))

图3 精确解和N=35、M=15时的逼近解对比图

图4 当N=15、M=12和N=30、M=14时逼近解与精确解之间的误差图像

从表2可以观察到，当N≥15且M≥12时，逼近解达到了10-13左右精度。另外，从图3和图4也观察到本文算法是收敛的和高精度的。

6 结语

本文提出了圆域上二阶变系数椭圆方程的一种高精度数值方法，并证明了弱解的存在唯一性和逼近解的误差估计。将该算法应用于奇异非线性二阶椭圆方程的计算，数值结果表明本文算法是收敛的和高精度的。但本文仅讨论了二阶椭圆方程，没有对一些更复杂的高阶非线性问题进行分析和讨论，这将是未来的研究目标。