在很多科学研究和工程计算中, 我们通常需要求解一些非线性的抛物问题和特征值问题, 如Allen-Cahn方程, Cahn-Hilliard方程和传输特征值问题 ^[1-6]. 然而, 这些非线性问题线性化后最终都归结为反复地求解一些二/四阶问题. 关于二/四阶问题的理论分析和数值计算已有很多成果, 如有限元法^[7-10]和谱方法^[11-14].

众所周知, 谱方法是一种具有谱精度的高阶数值方法, 但其对计算区域要求比较严格, 通常为乘积形式的矩形区域或长方体区域, 因此将该方法延伸到更多区域上的二/四阶问题的计算引起了越来越多学者的关注. 文献[15]提出了圆域上传输特征值问题基于降维格式的谱Galerkin逼近.文献[16]提出了圆域上四阶问题基于降维格式的一种有效的谱方法及误差分析. 文献[17]提出了球域上四阶问题的一种新的谱方法及误差分析. 文献[18]提出了球域上Maxwell传输特征值问题的一种有效的谱Galerkin逼近及误差分析. 文献[19]提出了圆域上四阶Steklov问题的一种有效的谱方法及误差估计. 然而, 这些数值方法主要都是基于径向变系数, 不能直接地应用于带有一般变系数或者一些特殊区域(如椭圆域)上的二/四阶问题.

据我们所知, 很少有关于椭圆域上二/四阶变系数问题的谱方法的报道, 其主要原因是极坐标变换引入了奇异变系数, 且径向$ r $和切向$ \theta $是耦合的, 不能对原问题进行降维求解, 这对理论分析和算法设计都带来了一定的困难. 因此, 本文的目的是提出椭圆域上二/四阶变系数问题的一种有效的谱Galerkin逼近. 首先, 我们将原问题化为极坐标下的等价形式, 并建立其弱形式及相应的离散格式. 其次, 针对二阶情形, 我们证明了弱解和逼近解的存在唯一性及它们之间的误差估计. 另外, 根据极条件和勒让得多项式的正交性, 我们构造了一组有效的径向基函数, 并在$ \theta $方向作截断的傅立叶展开, 推导了离散格式等价的矩阵形式. 最后, 我们给出了一些的数值算例, 数值结果表明了我们算法的收敛性和谱精度.

本文剩余部分组织如下: 在第2节, 我们推导了原问题在极坐标下的弱形式及离散格式. 在第3节, 我们证明了解的存在唯一性及误差估计. 在第4节, 我们详细描述了离散格式的有效实现. 在第5节, 我们给出了一些数值算例. 最后, 在第6节我们给出了一些结论性评注.

在这一节, 我们将推导椭圆域上二/四阶变系数问题在极坐标下的弱形式及其离散格式.

作为一个模型, 我们考虑下面的二阶问题:

$ \begin{align} -&\Delta u(x, y)+\alpha(x, y) u(x, y)=f(x, y), (x, y)\in\Omega , \end{align} $

(2.1)

$ \begin{align} & u(x, y)=0, (x, y)\in \partial \Omega, \end{align} $

(2.2)

其中$ \Omega =\big\{(x, y)\in \mathbb{R}^{2}: \frac{x^2}{a^2}+\frac{y^2}{b^2}<1 \big\} $, $ \alpha(x, y) $为非负有界函数. 令$ H^{s}(\Omega) $和$ H_0^{s}(\Omega) $表示通常的$ s $阶Sobolev空间, $ \|\cdot\|_s $和$ |\cdot|_s $分别表示其范数和半范数. 由格林公式可知(2.1)-(2.2)的弱形式为: 找$ u\in H^1_{0}(\Omega) $, 使得

$ \begin{align} & a(u, v)=F(v), \forall{v}\in H^{1}_{0}(\Omega) , \end{align} $

(2.3)

其中

$ \begin{align*} &a(u, v)=\int_{\Omega}\nabla{u}\nabla{\bar{v}}+\alpha u\bar{v}dxdy, F(v)=\int_{\Omega}f{\bar{v}}dxdy. \end{align*} $

令

$ \begin{align} x=\frac{1}{2}a(t+1)\cos\theta, y=\frac{1}{2}b(t+1)\sin\theta. \end{align} $

(2.4)

显然, (2.4)是一个从区域$ D=\big\{(t, \theta):t\in[-1, 1), \theta\in[0, 2\pi)\big\} $到区域$ \Omega $上的映射. 记$ \hat{u}(t, \theta)=u(x, y) $, 则由(2.4)及复合函数求导法则可得梯度算子在极坐标下的等价形式为:

$ \begin{align} & \nabla u=\hat{\nabla} \hat{u}=2\bigg[\frac{1}{a}(\cos\theta\partial_t\hat{u}-\frac{1}{t+1}\sin\theta\partial_\theta \hat{u}) , \frac{1}{b}(\sin\theta\partial_t\hat{u}+\frac{1}{t+1}\cos\theta\partial_\theta\hat{u})\bigg]. \end{align} $

(2.5)

为了使(2.3)适定, $ \hat{u} $需要满足下面的本质极条件:

$ \begin{align} &\partial_\theta\hat{u}(-1, \theta)=0. \end{align} $

(2.6)

定义Sobolev空间:

$ \begin{align*} \mathcal{H}^{1}_{0}(D)=\big\{&\hat{u}:{\int_{D}}\big(|\cos\theta\partial_t\hat{u}-\frac{1}{t+1}\sin\theta\partial_\theta\hat{u}|^2+|\sin\theta\partial_t\hat{u}+\frac{1}{t+1}\cos\theta\partial_\theta\hat{u}|^2\big)\\ &(t+1)dtd\theta<\infty, \hat{u}(1, \theta)=0, \partial_\theta\hat{u} (-1, \theta)=0, \hat{u}(t, \theta)=\hat{u}(t, \theta+2\pi)\big\}, \end{align*} $

其相应的内积和范数分别为:

$ \begin{align*} &\langle\hat{u}, {\hat{v}}\rangle_1={\int_D}(\cos\theta\partial_t\hat{u}-\frac{1}{t+1}\sin\theta\partial_\theta\hat{u})\cdot(\cos\theta\partial_t\bar{\hat{v}}-\frac{1}{t+1}\sin\theta\partial_\theta\bar{\hat{v}})(t+1)dtd\theta\\ &+{\int_D}(\sin\theta\partial_t\hat{u}+\frac{1}{t+1}\cos\theta\partial_\theta\hat{u})\cdot(\sin\theta\partial_t\bar{\hat{v}}+\frac{1}{t+1}\cos\theta\bar{\hat{v}})(t+1)dtd\theta, |||\hat{u}|||_1=\sqrt{\langle\hat{u}, \hat{u}\rangle_1}. \end{align*} $

由(2.5)可得弱形式(2.3)在极坐标下的等价形式为: 找$ \hat{u}\in\mathcal{H}^{1}_{0}(D) $, 使得

$ \begin{align} \mathcal{A}(\hat{u}, \hat{v})=\mathcal{F}(\hat{v}), \forall{\hat{v}}\in \mathcal{H}^{1}_{0}(D), \end{align} $

(2.7)

其中

$ \begin{align*} \mathcal{A}(\hat{u}, \hat{v})=&4{\int_D}\bigg[\frac{1}{a^2}(\cos\theta\partial_t\hat{u}-\frac{1}{t+1}\sin\theta\partial_\theta\hat{u})(\cos\theta\partial_t\bar{\hat{v}}-\frac{1}{t+1}\sin\theta\partial_\theta\bar{\hat{v}})\\ &+\frac{1}{b^2}(\sin\theta\partial_t\hat{u}+\frac{1}{t+1}\cos\theta\partial_\theta\hat{u})(\sin\theta\partial_t\bar{\hat{v}}+\frac{1}{t+1}\cos\theta\bar{\hat{v}})\bigg](t+1)dtd\theta+{\int_D}\hat{\alpha}\hat{u}\bar{\hat{v}}(t+1)dtd\theta, \\ \mathcal{F}(\hat{v})=&{\int_D}\hat{f}(t, \theta)\bar{\hat{v}}(t, \theta)(t+1)dtd\theta, \hat{\alpha}(t, \theta)=\alpha(x, y), \hat{f}(t, \theta)=f(x, y). \end{align*} $

由于$ \hat{u}(t, \theta) $在$ \theta $方向上是周期的, 则可由Fourier基函数展开如下:

$ \begin{align} \hat{u}(t, \theta)=\sum^{\infty}_{|m|=0}\hat{u}_{m}(t)e^{im\theta}. \end{align} $

(2.8)

将(2.8)代入(2.6)可得

$ \begin{align} \partial_\theta\hat{u}(-1, \theta)=\sum^{\infty}_{|m|=0}im \hat{u}_{m}(-1)e^{im\theta}=0. \end{align} $

(2.9)

由Fourier基函数的正交性, (2.9)可进一步简化为

$ \begin{align} m\hat{u}_{m}(-1)=0. \end{align} $

(2.10)

令$ P_{N} $为$ N $次多项式空间, 定义$ \mathcal{H}^{1}_{0}(D) $的一个逼近空间如下:

$ \begin{align*} X_{NM}=span\bigg\{\hat{u}_{mN}(t)e^{im\theta}:\hat{u}_{mN}(t)\in X_{N}^{m}, |m|\leq M\bigg\}, \end{align*} $

其中$ X_{N}^{m}=\{ p\in P_{N}: mp(-1)=p(1)=0\}. $则(2.7)式的一种离散格式为: 找$ u_{NM}\in X_{NM} $使得

$ \begin{align} \mathcal{A}(u_{NM}, v_{NM})=\mathcal{F}(v_{NM}), \forall v_{NM}\in X_{NM}. \end{align} $

(2.11)

我们考虑如下具有简支板边界条件的四阶问题:

$ \begin{align} &\Delta^2\sigma-\nabla(\alpha\nabla \sigma)+\beta \sigma=f, (x, y)\in\Omega, \end{align} $

(2.12)

$ \begin{align} &\sigma=\Delta \sigma=0, (x, y)\in\partial\Omega, \end{align} $

(2.13)

其中$ \alpha, \beta $为非负有界函数. 由格林公式可知, (2.12)-(2.13)的弱形式为: 找$ \sigma\in H_{0}^{1}(\Omega)\cap H^{2}(\Omega) $, 使得

$ \begin{align} G(\sigma, \varsigma)=F(\varsigma), \forall \varsigma\in H_{0}^{1}(\Omega)\cap H^{2}(\Omega), \end{align} $

(2.14)

其中$ G(\sigma, \varsigma)= \int_{\Omega}(\Delta \sigma\cdot\Delta\bar{\varsigma}+\alpha \nabla \sigma\nabla \bar{\varsigma}+\beta \sigma \bar{\varsigma})dxdy. $令$ \hat{\sigma}(t, \theta)=\sigma(x, y), \hat{\alpha}(t, \theta)=\alpha(x, y), \hat{\beta}(t, \theta)=\beta(x, y). $类似地, 利用复合函数求导法则，我们能够得到拉普拉斯算子在极坐变下的等价形式为:

$ \begin{align} \Delta \sigma=&\mathcal{L}\hat{\sigma}=2(\frac{1}{a^2}+\frac{1}{b^2})({\partial}^2_t\hat{\sigma}+\frac{1}{t+1}\partial_t\hat{\sigma}+\frac{1}{(t+1)^2}{\partial}^2_\theta\hat{\sigma})\\ &+2(\frac{1}{a^2}-\frac{1}{b^2})[\cos2\theta({\partial}^2_t\hat{\sigma}-\frac{1}{t+1}\partial_t\hat{\sigma}-\frac{1}{(t+1)^2}{\partial}^2_\theta\hat{\sigma})+\frac{2}{t+1}\sin 2\theta(\frac{1}{t+1}\partial_\theta\hat{\sigma}-\partial_t\partial_\theta\hat{\sigma})]. \end{align} $

(2.15)

为了使(2.14)适定, $ \hat{\sigma} $需要满足下面的本质极条件

$ \begin{align} {\partial}^2_\theta\hat{\sigma}(-1, \theta)=0, \lim\limits_{t\rightarrow -1}(\partial_t\hat{\sigma}+\frac{1}{t+1}{\partial}^2_\theta\hat{\sigma})=0. \end{align} $

(2.16)

定义Sobolev空间:

$ \begin{align*} \mathcal{H}_{*}^{2}(D)=\bigg\{&\hat{\sigma}:\int_{D}(t+1)|\mathcal{L}\hat{\sigma}|^2dtd\theta<\infty, \hat{\sigma}(1, \theta)={\partial}^2_\theta\hat{\sigma}(-1, \theta)=0;\\ &\lim\limits_{t\rightarrow -1}(\partial_t\hat{\sigma}+\frac{1}{t+1}{\partial}^2_\theta\hat{\sigma})=0 , \hat{\sigma}(t, \theta)=\hat{\sigma}(t, \theta+2\pi)\bigg\}, \end{align*} $

相应的内积和范数分别为:

$ \begin{align*} &\langle\hat{\sigma}, \hat{\varsigma}\rangle_2=\int_{D}\mathcal{L}\hat{\sigma}\cdot\mathcal{L}\bar{\hat{\varsigma}}(t+1)dtd\theta, |||\hat{\sigma}|||_2=\sqrt{\langle\hat{\sigma}, \hat{\sigma}\rangle_2}. \end{align*} $

由(2.5)和(2.15)可得弱形式(2.14)在极坐标下的等价形式为: 找$ \hat{\sigma}\in \mathcal{H}_{*}^{2}(D), $使得

$ \begin{align} \mathcal{G}(\hat{\sigma}, \hat{\varsigma})=\mathcal{F}(\hat{\varsigma}), \forall \hat{\varsigma}\in \mathcal{H}_{*}^{2}(D). \end{align} $

(2.17)

其中

$ \begin{align*} &\mathcal{G}(\hat{\sigma}, \hat{\varsigma})=\int_{D}[\mathcal{L}\hat{\sigma}\mathcal{L}\bar{\hat{\varsigma}}+\alpha(t, \theta)\hat{\nabla}\hat{\sigma}\hat{\nabla}\bar{\hat{\varsigma}}+\beta(t, \theta)\hat{\sigma}\bar{\hat{\varsigma}}](t+1)dtd\theta. \end{align*} $

类似地, 我们将$ \hat{\sigma} $在$ \theta $方向用Fourier基函数展开:

$ \begin{align} \hat{\sigma}(t, \theta)=\sum\limits_{|m|=0}^{\infty}\hat{\sigma}_{m}(t)e^{im\theta}. \end{align} $

(2.18)

将(2.18)代入(2.16)可得:

$ \begin{align*} m^2\hat{\sigma}_m(-1)=(1-m^2)\partial_t\hat{\sigma}_m(-1)=0. \end{align*} $

定义$ \mathcal{H}^{2}_{*}(D) $的一个逼近空间如下:

$ \begin{align*} \mathcal{X}_{NM}=span\bigg\{\hat{\sigma}_{mN}(t)e^{im\theta}:\hat{\sigma}_{mN}(t)\in \mathcal{X}_{N}^{m}, |m|\leq M\bigg\}, \end{align*} $

其中

$ \begin{align*} \mathcal{X}_{N}^{m}=\{ p\in P_{N}: m^2\hat{p}(-1)=(1-m^2)\partial_t\hat{p}(-1)=\hat{p}(1)=0\}. \end{align*} $

则(2.17)的一种离散格式为: 找$ \sigma_{NM}\in \mathcal{X}_{NM}, $使得

$ \begin{align} \mathcal{G}(\sigma_{NM}, \varsigma_{NM})=\mathcal{F}(\varsigma_{NM}), \forall \varsigma_{NM}\in \mathcal{X}_{NM}. \end{align} $

(2.19)

为了简洁起见, 我们仅给出二阶情况下解的存在唯一性证明及误差估计.

引理1    $ \mathcal{A}(\hat{u}, \hat{v}) $为定义在$ \mathcal{H}^{1}_{0}(D)\times\mathcal{H}^{1}_{0}(D) $上的连续强制的双线性形式, 即

$ \begin{align*} &|\mathcal{A}(\hat{u}, \hat{v})|\lesssim |||\hat{u}|||_1|||\hat{v}|||_1, \forall \hat{u}, \hat{v}\in \mathcal{H}^{1}_{0}(D); \\ &\mathcal{A}(\hat{u}, \hat{u})\gtrsim |||\hat{u}|||_1^2, \forall \hat{u} \in \mathcal{H}^{1}_{0}(D). \end{align*} $

证   由于$ \alpha $是非负有界函数, 利用Cauchy-Schwarz不等式可得:

$ \begin{align*} |\mathcal{A}(\hat{u}, \hat{v})|=&4{\int_D}\bigg[\frac{1}{a^2}(\cos\theta\partial_t\hat{u}-\frac{1}{t+1}\sin\theta\partial_\theta\hat{u})(\cos\theta\partial_t\bar{\hat{v}}-\frac{1}{t+1}\sin\theta\partial_\theta\bar{\hat{v}})\\ &+\frac{1}{b^2}(\sin\theta\partial_t\hat{u}+\frac{1}{t+1}\cos\theta\partial_\theta\hat{u})(\sin\theta\partial_t\bar{\hat{v}}+\frac{1}{t+1}\cos\theta\bar{\hat{v}})\bigg](t+1)dtd\theta\\ &+{\int_D}\hat{\alpha}\hat{u}\bar{\hat{v}}(t+1)dtd\theta=\frac{2}{ab}\bigg|{\int_\Omega}\nabla u\nabla\bar{ v}+\alpha u\bar{v}dxdy\bigg|\lesssim \|u\|_1\|v\|_1. \end{align*} $

由Poincaré不等式有

$ \begin{align*} |\mathcal{A}(\hat{u}, \hat{v})|\lesssim &\|u\|_1\|v\|_1\lesssim |u|_1|v|_1 =({\int_\Omega}|\nabla u|^2dxdy)^\frac{1}{2}({\int_\Omega}|\nabla v|^2dxdy)^\frac{1}{2}\\ \lesssim&\bigg\{{\int_D}(|\cos\theta\partial_t\hat{u}-\frac{1}{t+1}\sin\theta\partial_\theta\hat{u}|^2+|\sin\theta\partial_t\hat{u}+\frac{1}{t+1}\cos\theta\partial_\theta\hat{u}|^2)(t+1)dtd\theta\bigg\}^\frac{1}{2}\\ &\bigg\{{\int_D}(|\cos\theta\partial_t\hat{v}-\frac{1}{t+1}\sin\theta\partial_\theta\hat{v}|^2+|\sin\theta\partial_t\hat{v}+\frac{1}{t+1}\cos\theta\partial_\theta\hat{v}|^2)(t+1)dtd\theta\bigg\}^\frac{1}{2}\\ =&|||u|||_1|||v|||_1. \end{align*} $

另一方面, 我们有

$ \begin{align*} \mathcal{A}(\hat{u}, \hat{u})&\geq4{\int_D}\bigg(\frac{1}{a^2}|\cos\theta\partial_t\hat{u}-\frac{1}{t+1}\sin\theta\partial_\theta\hat{u}|^2+\frac{1}{b^2}|\sin\theta\partial_t\hat{u}+\frac{1}{t+1}\cos\theta\partial_\theta\hat{u}|^2\bigg)(t+1)dtd\theta\\ &\gtrsim{\int_\Omega}(|\cos\theta\partial_t\hat{u}-\frac{1}{t+1}\sin\theta\partial_\theta\hat{u}|^2+|\sin\theta\partial_t\hat{u}+\frac{1}{t+1}\cos\theta\partial_t\hat{u}|^2)(t+1)dtd\theta=|||\hat{u}|||_{1}^{2}. \end{align*} $

证毕.

定理1   若$ f\in L^{2}(\Omega) $, 则弱形式(2.7)和离散格式(2.11)分别存在唯一解$ \hat{u}\in\mathcal{H}^{1}_{0}(D) $和$ u_{NM}\in X_{NM} $.

证明: 由于$ f\in L^{2}(\Omega) $, 则由Cauchy-Schwarz不等式及Poincaré不等式有

$ \begin{align*} |\mathcal{F}(\hat{v})|&=|{\int_D}\hat{f}(r, \theta)\bar{\hat{v}}(t+1)dtd\theta|=\frac{2}{ab}|{\int_\Omega}f\bar{v}dxdy|\\ &\leqslant\frac{2}{ab}({\int_\Omega}|f|^2dxdy)^\frac{1}{2}({\int_\Omega}|v|^2dxdy)^\frac{1}{2}\lesssim ({\int_\Omega}|v|^2dxdy)^\frac{1}{2}\\ &\lesssim ({\int_\Omega}|\nabla v|^2dxdy)^\frac{1}{2}\lesssim|||\hat{u}|||_{1}. \end{align*} $

即$ \mathcal{F}(\hat{v}) $为定义在$ \mathcal{H}^{1}_{0}(D) $上的有界线性泛函. 由引理1及Lax-Milgram定理可知弱形式(2.7)存在唯一解$ \hat{u}\in\mathcal{H}^{1}_{0}(D) $. 类似地, 我们可证离散格式(2.11) 存在唯一解$ u_{NM}\in X_{NM} $. 证毕.

下面将证明逼近解的误差估计. 首先, 我们证明下面的Céa引理.

引理2   设$ \hat{u} $和$ u_{NM} $分别为弱形式(2.7)和离散格式(2.11)的解, 则有下列不等式成立:

$ \begin{align*} |||\hat{u}-u_{NM}|||_1\lesssim \inf\limits_{v_{NM} \in X_{NM}}|||\hat{u}-v_{NM}|||_1. \end{align*} $

证明: 在(2.7)中取$ \hat{v}=\hat{v}_{NM} $, 结合(2.11)可得

$ \begin{align} \mathcal{A}(\hat{u}-u_{NM}, v_{NM})=0, \forall v_{NM}\in X_{NM}. \end{align} $

(3.1)

由(3.1)和引理1, 我们有

$ \begin{align*} |||\hat{u}-u_{NM}|||_{1}^{2}&\lesssim \mathcal{A}(\hat{u}-u_{NM}, \hat{u}-u_{NM})\\ &=\mathcal{A}(\hat{u}-u_{NM}, \hat{u}-v_{NM}+v_{NM}-u_{NM})\\ &=\mathcal{A}(\hat{u}-u_{NM}, \hat{u}-v_{NM})\\ &\lesssim|||\hat{u}-u_{NM}|||_1|||\hat{u}-v_{NM}|||_1. \end{align*} $

由$ v_{NM} $的任意性可知结论成立. 证毕.

定义如下的Sobolev空间:

$ \begin{align*} H_p^{s}(0, 2\pi)=\{\phi\in H^{s}(0, 2\pi):\phi(\theta+2\pi)=\phi(\theta)\}, \end{align*} $

其内积和范数分别为：

$ \begin{align*} &(\phi, \xi)_s=\int^{2\pi}_{0}\sum^{s}_{i=0}\partial_\theta^k\phi\cdot\partial_\theta^k\bar{\xi}d\theta, ||\phi||_s=\sqrt{(\phi, \phi)_s}. \end{align*} $

引理3^[14]  对$ \forall \phi(\theta)\in H^{\kappa}_{p}(0, 2\pi) $及$ 0\leq\ell\leq \kappa $, 下列不等式成立:

$ \begin{align*} ||\phi(\theta)-\phi_M(\theta)||_\ell\lesssim M^{\ell-\kappa}|\phi(\theta)|_\kappa, \end{align*} $

其中

$ \begin{align*} \phi_{M}(\theta)=\sum\limits_{|m|=0}^{M}\phi_{m}e^{im\theta}, \phi_{m}=\int_0^{2\pi}\phi(\theta) e^{-im\theta}d\theta. \end{align*} $

令$ \omega^{k, l}(t)=(1-t)^k(1+t)^l $为定义在$ I=(-1, 1) $上的一个权函数. 定义正交投影$ \pi_{N, \omega^{-1, -1}}: L_{\omega^{-1, -1}}^2(I)\rightarrow P_N^0=\{p\in P_N: p(\pm1)=0\} $, 使得

$ \begin{align*} (\psi-\pi_{N, \omega^{-1, -1}}\psi, \varrho_N)_{\omega^{-1, -1}}=0, \forall \varrho_N\in P_N^0. \end{align*} $

定义非一致带权Sobolev空间:

$ \begin{align*} H_{\omega^{-1, -1}, *}^{s}(I)=\{\psi:\partial_t^j\psi\in L_{\omega^{-1+j, -1+j}}^{2}(I), 0\leq j\leq s \}, \end{align*} $

其内积和范数分别为:

$ \begin{align*} (\psi, \phi)_{s, \omega^{-1, -1}, *}=\sum^{s}_{j=0}(\partial_t^j\psi, \partial_t^j\phi)_{\omega^{-1+j, -l+j}}, \|\psi\|_{s, \omega^{-1, -1}, *}=\sqrt{(\psi, \psi)_{s, \omega^{-1, -1}, *}} . \end{align*} $

引理4^[13]   对于任意$ \psi\in H_{\omega^{-1, -1}, *}^{s}(I), $下列估计式成立:

$ \begin{align*} \|\partial_t^\mathcal{\mu}(\pi_{N, \omega^{-1, -1}}\psi-\psi)\|_{\omega^{-1+\mathcal{\mu}, -1+\mathcal{\mu}}}\lesssim N^{\mathcal{\mu}-s}\|\partial_t^s\psi\|_{\omega^{-1+s, -1+s}}, 0\leq \mathcal{\mu}\leq s. \end{align*} $

定义Sobolev空间:

$ \begin{align*} &H_{0, m}^{1}(I)=\{\hat{u}_{m}:\int_{I}(t+1)|{\partial_t\hat{u}_m}|^2+\frac{m^2}{t+1}{|\hat{u}_{m}|}^2dt<\infty, m^2\hat{u}_{m}(-1)=\hat{u}_{m}(1)=0\};\\ &\mathcal{H}_{0, m}^s(I):=\{\hat{u}_{m}\in H_{0, m}^{1}(I):\partial_t^k\hat{u}_m\in L_{\omega^{-1+k, -1+k}}^2, 1\leq k\leq s\}, \end{align*} $

其相应的内积和范数分别为:

$ \begin{align*} &(\hat{u}_m, \hat{v}_m)_{1, m}=\int_I(t+1)\partial_t\hat{u}_{m}\partial_t\bar{\hat{ v}}_{m}+\frac{m^2}{t+1}\hat{u}_{m} \bar{\hat{v}}_{m}dt, \|\hat{u}_m\|_{1, m}= \sqrt{(\hat{u}_m, \hat{u}_m)_{1, m}};\\ &[\hat{u}_m, \hat{v}_m]_{s, m}=(\hat{u}_m, \hat{v}_m)_{1, m}+\sum\limits_{k=1}^{s}(\partial_t^k\hat{u}_m, \partial_t^k\hat{v}_m)_{\omega^{-1+k, -1+k}}, \||\hat{u}_m\||_{s, m}=[\hat{u}_m, \hat{u}_m]_{s, m}^{\frac{1}{2}}. \end{align*} $

引理5   存在一个算子$ \pi_N^{1, m}:H_{0, m}^1(I)\rightarrow X_N(m)=\{p\in P_{N}:mp(-1)=p(1)=0\} $使得$ \pi_N^{1, m}\hat{u}_m(\pm1)=\hat{u}_m(\pm1) $, 且对任意$ \hat{u}_m\in \mathcal{H}_{0, m}^s(I) $, $ s\geq 1 $, 下列估计式成立：

$ \begin{align*} \|\partial_t^k(\pi_{N}^{1, m}\hat{u}_m-\hat{u}_m)\|_{\omega^{-1+k, -1+k}}\lesssim N^{k-s}\|\partial_t^s\hat{u}_m\|_{\omega^{-1+s, -1+s}}, 0\leq k\leq s. \end{align*} $

证   令$ \hat{u}_{m}^*=\frac{1-t}{2}\hat{u}_m(-1) $, 则有$ (\hat{u}_m-\hat{u}_{m}^*)(\pm1)=0. $对于$ \forall \hat{u}_{m} \in \mathcal{H}_{0, m}^s(I) $, 我们有$ \hat{u}_{m}-\hat{u}_{m}^*\in H_{\omega^{-1, -1}, *}^s(I) $. 事实上, 由Hardy不等式(参见文献[14] 中的B.8.8), 我们有

$ \begin{align} &\int_I\omega^{-1, -1}(\hat{u}_m-\hat{u}_{m}^*)^2dt\lesssim\int_I[\partial_t(\hat{u}_m-\hat{u}_{m}^*)]^2dt. \end{align} $

(3.2)

由Cauchy-Schwarz不等式可得

$ \begin{align*} \int_I(\partial_t\hat{u}_{m}^*)^2dt&=\int_I\frac{1}{4}[-\hat{u}_m(-1)]^2dt= \frac{1}{2}(\int_I\partial_t\hat{u}_mdt)^2\leq\int_I(\partial_t\hat{u}_m)^2dt. \end{align*} $

因此, 我们有

$ \begin{align} \int_I[\partial_t(\hat{u}_m-\hat{u}_{m}^*)]^2dt&\lesssim\int_I(\partial_t\hat{u}_m)^2dt+\int_I(\partial_t\hat{u}_{m}^*)^2dt \lesssim \int_I(\partial_t\hat{u}_m)^2dt. \end{align} $

(3.3)

对于$ k\geq2 $, 我们有

$ \begin{align} \int_I\omega^{-1+k, -1+k}[\partial_t^k(\hat{u}_m-\hat{u}_{m}^*)]^2dt&=\int_I\omega^{-1+k, -1+k}(\partial_t^k\hat{u}_m)^2dt. \end{align} $

(3.4)

由(3.2)-(3.4)可得$ \hat{u}_{m}-\hat{u}_{m}^*\in H_{\omega^{-1, -1}, *}^s(I) $. 定义投影算子

$ \begin{align*} \pi_N^{1, m}\hat{u}_{m}=\pi_{N, \omega^{-1, -1}}(\hat{u}_{m}-\hat{u}_{m}^*)+\hat{u}_{m}^* \in X_N(m), \forall \hat{u}_{m}\in \mathcal{H}_{0, m}^s(I). \end{align*} $

由引理4可得

$ \begin{align*} \|\partial_t^k(\pi_N^{1, m}\hat{u}_{m}-\hat{u}_{m})\|_{\omega^{-1+k, -1+k}}&=\|\partial_t^k(\pi_{N, \omega^{-1, -1}}(\hat{u}_m-\hat{u}_{m}^*)-(\hat{u}_m-\hat{u}_{m}^*))\|_{\omega^{-1+k, -1+k}}\\ &\lesssim N^{k-s}\|\partial_t^s(\hat{u}_m-\hat{u}_{m}^*)\|_{\omega^{-1+s, -1+s}}. \end{align*} $

结合(3.3)-(3.4)可知结论成立. 证毕.

记

$ \begin{align*} &u_M(t, \theta)=\sum^{M}_{|m|=0}\hat{u}_m(t)e^{im\theta}, v_{NM}(t, \theta)=\sum^{M}_{|m|=0}\hat{v}_{mN}(t)e^{im\theta}. \end{align*} $

定理2    设$ \hat{u} $和$ u_{NM} $分别为弱形式(2.7)和离散格式(2.11) 的解. 若对$ \forall t\in I, \hat{u}(t, \theta)\in H_{p}^\kappa(0, 2\pi), $且$ \hat{u}_m\in \mathcal{H}_{0, m}^s(I) $, 则有下列不等式成立:

$ \begin{align*} |||\hat{u}-u_{NM}|||_1\lesssim M^{1-\kappa}+N^{1-s}. \end{align*} $

证    由于

$ \begin{align*} |||\hat{u}-v_{NM}|||_{1}^{2}=&\int_{D}|\partial_t(\hat{u}-v_{NM})|^2(t+1)+\frac{1}{t+1}|\partial_\theta(\hat{u}-v_{NM})|^2dtd\theta\\ \leq& 2\int_{D}(1+t)|\partial_t(\hat{u}-u_M)|^2+(1+t)|\partial_t(u_M-v_{NM})|^2dtd\theta\\ &+2\int_{D}\frac{1}{t+1}|\partial_\theta(\hat{u}-u_M)|^2+\frac{1}{t+1}|\partial_\theta(u_M-v_{NM})|^2dtd\theta. \end{align*} $

而

$ \begin{align*} &|\partial_t(u_M-v_{NM})|^2\leq[\sum\limits_{|m|=0}^{M}|\partial_t(\hat{u}_m-\hat{v}_{mN})|]^2 \leq(2M+1)\sum\limits_{|m|=0}^{M}|\partial_t(\hat{u}_m-\hat{v}_{mN})|^2.\\ &|\partial_\theta(\hat{u}_M-v_{MN})|^2=|\sum\limits_{|m|=1}^{M}(\hat{u}_m-\hat{v}_{mN})e^{im\theta}\cdot im|^2 \leq 2M\sum\limits_{|m|=1}^{M}m^2|\hat{u}_m-\hat{v}_{mN}|^2. \end{align*} $

则有

$ \begin{align*} |||\hat{u}-v_{NM}|||_1^2\leq& 4\int_{D}|\partial_t(\hat{u}-u_M)|^2dtd\theta+4\pi(2M+1)\sum\limits_{|m|=0}^{M}\int_{I}(1+t)|\partial_t(\hat{u}_m-\hat{v}_{mN})|^2dt\\ &+2\int_{D}\frac{1}{t+1}|\partial_\theta(\hat{u}-u_M)|^2dtd\theta+16\pi M\sum\limits_{|m|=1}^{M}m^2\int_{I}\frac{1}{1-t^2}|\hat{u}_m-\hat{v}_{mN}|^2dt. \end{align*} $

由引理3和引理5可得

$ \begin{align*} \inf\limits_{v_{NM} \in X_{NM}}|||\hat{u}-v_{NM}|||_1^2\lesssim& \int_{I}dt\int_{0}^{2\pi}|\partial_t(\hat{u}-u_M)|^2d\theta+M\sum\limits_{|m|=0}^{M}\int_{I}|\partial_t(\hat{u}_m-\pi_{N}^{1, m}\hat{u}_m)|^2dt\\ &+\int_{D}\frac{1}{t+1}|\partial_\theta(\hat{u}-u_M)|^2dtd\theta+M\sum\limits_{|m|=1}^{M}m^2\int_{I}\omega^{-1, -1}|\hat{u}_m-\pi_{N}^{1, m}\hat{u}_m|^2dt\\ \lesssim& M^{-2\kappa}\int_{D}\bigg|\frac{\partial^{\kappa}\partial_t\hat{u}}{\partial_\theta^\kappa}\bigg|^2d\theta dt +N^{2(1 -s)}\|\partial_t^s\hat{u}_m\|_{\omega^{-1+s, -1+s}}^2\\ &+M^{2(1-\kappa)}\int_{D}\frac{1}{t+1}\bigg|\frac{\partial^{\kappa}\hat{u}}{\partial_\theta^{\kappa}}\bigg|^2d\theta dt +N^{-2s}\|\partial_t^s\hat{u}_m\|_{\omega^{-1+s, -1+s}}^2\\ \lesssim& M^{2(1-\kappa)}+N^{2(1-s)}. \end{align*} $

由引理2可知结论成立. 证毕.

在这一节, 我们将建立离散格式(2.11)和(2.19)等价的矩阵形式. 首先, 我们分别构造逼近空间$ X_{NM} $和$ \mathcal{X}_{NM} $中的一组基函数.

设$ L_k(t) $为$ k $阶Legendre多项式. 令

$ \begin{align*} &\varphi_k(t)=L_k(t)-L_{k+2}(t), i=0, 1, \cdots, N-2, \varphi_{N-1}(t)=\frac{1}{2}(1-t). \end{align*} $

则有

$ \begin{align*} &X_{NM}=span\{\varphi_k(t) e^{im\theta}: k=0, 1, \cdots, N-1-sgn(|m|), |m|\leq M \}. \end{align*} $

我们将$ u_{NM} $展开如下:

$ \begin{align} u_{NM}=\sum\limits_{\lvert m\rvert=0}^{M}\sum\limits_{j=0}^{N-1-sign(|m|)}u_{jm}\varphi_{j}(t)e^{im\theta}. \end{align} $

(4.1)

将(4.1)代入离散格式(2.11), 让$ v_{NM} $取遍$ X_{NM} $中的一组基函数, 则可得到如下的线性系统:

$ \begin{align} \mathbf{A}\mathbf{u}=\mathbf{f}, \end{align} $

(4.2)

其中

$ \begin{align*} &\mathbf{A}=(a_{kjnm}), a_{kjnm}=\mathcal{A}(\varphi_{j}(t)e^{im\theta}, \varphi_{k}(t)e^{in\theta}), \mathbf{f}=(f_{kn}), f_{kn}=\int_{D}(t+1)\hat{f}(t, \theta)\varphi_{k}(t)e^{-in\theta}dtd\theta.\\ &\mathbf{u}=[u_{0, -M}, u_{1, -M}, \cdots, u_{{N-2}, -M}, \cdots, u_{0, 0}, u_{1, 0}, \cdots , u_{N-1, 0}, \cdots, u_{0, M}, u_{1, M}, \cdots, u_{N-2, M}]^T. \end{align*} $

注: 当$ \alpha $为非负常数时, 由Legendre多项式和Fourier基函数的正交性可知(4.2)中系数矩阵都是分块带状稀疏矩阵, 从而可用一些稀疏矩阵求解器有效地求解(4.2). 当$ \alpha $为非负函数时, 系数矩阵$ \mathbf{A} $有可能是满的, 这时我们可采用预条件迭代法对(4.2)进行求解.

构造一组基函数如下:

$ \begin{align*} &\tilde{\phi}_i(t)=(1-t^2)[L_i(t)-L_{i+2}(t)], i=0, 1, \cdots, N-4, \tilde{\phi}_{N-3}(t)=\frac{1}{4}(1-t^2)(t+1), \\ &\tilde{\phi}_{N-2}(t)=\frac{1}{4}(1-t)(t+3), \phi_{i}(t)=\tilde{\phi}_i(t), i=0, 1, \cdots, N-4, \\ & \phi_{N-3}(t)=\frac{1}{4}(1-t^2)(t+1), \phi_{N-2}(t)=\frac{1}{4}(1-t^2). \end{align*} $

则有

$ \begin{align*} \mathcal{X}_{NM}&=span\{\tilde{\phi}_k(t), k=0, 1, \cdots, N-2;\phi_k(t) e^{im\theta}, \\ k&=0, \cdots, N-2-sign(|m|-1), 1\leq|m|\leq M \}. \end{align*} $

我们将$ \sigma_{NM} $展开如下:

$ \begin{align} \sigma_{NM}=\sum\limits_{|m|=1}^{M}\sum\limits_{j=0}^{N-2-sign(|m|-1)}\sigma_{jm}\phi_{j}(t)e^{im\theta}+\sum\limits_{j=0}^{N-2}\sigma_{j_0}\tilde{\phi}_{j}(t). \end{align} $

(4.3)

将(4.3)代入离散格式(2.19), 让$ \varsigma_{NM} $取遍$ \mathcal{X}_{NM} $中的一组基函数, 则可得到如下的线性系统:

$ \begin{align} \mathbf{B}\mathbf{\sigma}=\mathbf{g}, \end{align} $

(4.4)

其中

$ \begin{align*} &B=(b_{kjnm}), b_{kjnm}=\mathcal{G}(\phi_{jm}(t)e^{im\theta}, \phi_{kn}(t)e^{in\theta}), \phi_{jm}(t)=\tilde{\phi}_{j}(t), (m=0), \\ &\phi_{jm}(t)=\phi_{j}(t), (m\neq0); \phi_{jn}(t)=\tilde{\phi}_{j}(t), (n=0), \phi_{jn}(t)=\phi_{j}(t), (n\neq0), \\ & \mathbf{g}=(g_{kn}), g_{kn}=\int_{D}(t+1)\hat{f}(t, \theta)\varphi_{kn}(t)e^{-in\theta}dtd\theta, \\ &\sigma=[\sigma_{0, -M}, \cdots, \sigma_{{N-3}, -M};\cdots;u_{0, -1}, \cdots\sigma_{N-2, -1};u_{0, 0}, \cdots\sigma_{N-2, 0}; \\ &\qquad u_{0, 1}, \cdots\sigma_{N-2, 1}; \cdots; \sigma_{0, M}, \cdots, \sigma_{N-3, M}]^T. \end{align*} $

为了表明算法的收敛性和谱精度, 我们将给出一些数值算例, 并在MATLAB软件R2016b版本上进行编程计算. 定义精确解与逼近解在$ L^{\infty} $意义下的误差如下:

$ \begin{align*} &e(u, u_{NM})=\|\hat{u}(t, \theta)-u_{NM}(t, \theta)\|_{L^\infty(D)}, e(\sigma, \sigma_{NM})=\|\hat{\sigma}(t, \theta)-\sigma_{NM}(t, \theta)\|_{L^\infty(D)}. \end{align*} $

例1：我们考虑变系数二阶问题(2.1)-(2.2). 取$ a=3, b=1, \alpha(x, y)=3+sin(x+y) $和给定的$ f $, 使得精确解$ u(x, y)=\cos(x + y)sin[\pi(\frac{x^2}{9} + y^2 - 1)] $. 对于不同的$ N $和$ M={N}/{2} $, 我们在表 1中列出了精确解$ u $与逼近解$ u_{NM} $之间的误差$ e(u, u_{NM}) $. 从表 1我们可以观察到, 当$ N\geq30 $时, 逼近解与精确解之间的误差达到了大约$ 10^{-11} $的精度.为了更直观地表明我们算法的收敛性和谱精度, 我们分别在图 1和图 2中列出了精确解与$ N=30 $时逼近解的图像, 在图 3和 图 4中分别列出了$ N=30 $和$ N=40 $时精确解与逼近解之间的误差图像.从图 1-4我们可以观察到我们的算法是收敛的和谱精度的。

例2：我们考虑四阶问题(2.12)-(2.13). 取$ a=3, b=1, \alpha=1, \beta=3 $和给定的$ f $, 使得精确解$ \sigma(x, y)=(\frac{x^2}{9} + y^2 - 1)^2sin[\pi(\frac{x^2}{9} + y^2 - 1)] $. 对于不同的$ N $和$ M={N}/{2} $, 我们在表 2中列出了精确解$ \sigma $与逼近解$ \sigma_{NM} $之间的误差$ e(\sigma, \sigma_{NM}) $, 从表 2我们可以观察到, 当$ N\geq26 $时, 逼近解与精确解之间的误差达到了大约$ 10^{-15} $的精度. 类似地, 我们分别在图 5和图 6中列出了精确解与$ N=30 $时逼近解的图像, 在图 7和图 8中分别列出了$ N=30 $和$ N=40 $时精确解与逼近解之间的误差图像. 从图 5-8我们再次观察到我们的算法是收敛的和谱精度的.

本文提出了椭圆域上二/四阶方程的一种有效的谱Galerkin逼近.针对二阶问题, 我们从理论上证明了弱解和逼近解的存在唯一性及它们之间的误差估计.该方法不仅具有谱精度，而且还克服了有限元求解过程中区域剖分及其基函数构造的复杂性, 尤其是椭圆域上的一些四阶问题.

另外，本文提出的算法还可应用于一些非线性的抛物方程, 如Cahn-Hilliard方程, Schrödinger方程等, 这是我们将来的研究目标.