贵州师范大学数学科学学院，贵州 贵阳

收稿日期：2022年6月4日；录用日期：2022年6月29日；发布日期：2022年7月6日

摘要

本文针对周期边界条件下变系数非线性二阶问题提出了一种有效的Fourier谱方法。首先，根据边界条件引入了适当的Sobolev空间及其逼近空间，建立了变系数非线性二阶问题的弱形式和相应的离散格式。基于这非线性的离散格式，我们建立了一种线性迭代算法，并给出了该算法相应的Matlab程序设计。最后，我们给出了数值算例，数值结果表明我们提出的算法是收敛的和高精度的。

关键词

二阶非线性问题，周期边界条件，Fourier谱方法，程序设计，数值实验

Efficient Fourier Spectral Approximation for Nonlinear Second-Order Problems with Variable Coefficients

Tingting Jiang

School of Mathematical Sciences, Guizhou Normal University, Guiyang Guizhou

Received: Jun. 4^th, 2022; accepted: Jun. 29^th, 2022; published: Jul. 6^th, 2022

ABSTRACT

In this paper, an efficient Fourier spectral method is proposed for nonlinear second-order problems with variable coefficients under periodic boundary conditions. Firstly, an appropriate Sobolev space and its approximation space are introduced according to the boundary conditions, and the weak form and the corresponding discrete scheme of the nonlinear second-order problem with variable coefficients are established. Based on the nonlinear discrete scheme, we establish a linear iterative algorithm and its Matlab program design. Finally, we give a numerical example, and the numerical results show that our proposed algorithm is convergent and highly accurate.

Keywords:Second-Order Nonlinear Problems, Periodic Boundary Conditions, Fourier Spectral Method, Program Design, Numerical Experiments

Copyright © 2022 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

很多科学和工程问题最终都归结为求解非线性偏微分方程，如等离子体物理、非线性光学，激光脉冲中的自聚焦、热脉冲在晶体中的传播以及在极低温下玻色–爱因斯坦凝聚的动力学可由非线性薛定谔方程来描述 [1] [2] [3] [4] [5]；在材料科学和流体动力学中的许多复杂的运动界面问题可由Allen-Cahn和Cahn-Hilliard方程来描述 [6] [7] [8] [9] [10]。因此，提出一种有效求解非线性偏微分方程的高精度数值方法是非常有意义的。

到目前为止，已有很多数值方法求解非线性偏微分方程 [11] - [18]。但它们主要都是基于有限元方法和有限差分方法，要获得高精度的数值解需要花费很多计算时间和内存容量。在 [19] 中，非线性特征值问题在矩形网格一定条件下，证明了光谱的高精度。据我们所知，很少有关于周期边界条件下变系数非线性二阶问题的Fourier谱方法的报道。因此，本文的目的是针对周期边界条件下变系数非线性二阶问题提出了一种有效的Fourier谱方法。首先，根据边界条件引入了适当的Sobolev空间及其逼近空间，建立了变系数非线性二阶问题的弱形式和相应的离散格式。基于这非线性的离散格式，我们建立了一种线性迭代算法，并给出了该算法相应的Matlab程序设计。最后，我们给出了数值算例，数值结果表明我们提出的算法是收敛的和高精度的。

本文剩余部分安排如下：在第二节中，我们推导了变系数非线性二阶问题的弱形式及其离散格式。在第三节中，我们详细描述了算法的实现过程及其程序设计。在第四节中，我们给出了一些数值算例。

2. 弱形式及其离散格式

作为一个模型，我们考虑如下的变系数非线性二阶问题：

$-\Delta u+V\left(x,y\right)u+h\left(u^2\right)u=f,\left(x,y\right)\in \Omega ,$ (2.1)

$u\left(x,y\right)=u\left(x+2\pi ,y\right),u\left(x,y\right)=u\left(x,y+2\pi \right),$ (2.2)

其中 $V\left(x,y\right)$ 为非负变系数， $h\left(u^2\right)$ 为非线性项， $\Omega =\left(0,2\pi \right)\times \left(0,2\pi \right)$ 为计算区域。

下面我们将推导(2.1)~(2.2)的弱形式及其离散格式，用 $H^s\left(\Omega \right)$ 表示s阶Sobolev空间， ${\Vert \cdot \Vert }_s$ 表示 $H^s\left(\Omega \right)$ 中的范数。特别地，我们有

$H^0\left(\Omega \right)=L^2\left(\Omega \right)=\left\{u:{\displaystyle {\int }_{\Omega }{\left|u\right|}^2\text{d}x\text{d}y}<\infty \right\},$

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

$\left(u,v\right)={\displaystyle {\int }_{\Omega }u\bar{v}\text{d}x\text{d}y},\Vert u\Vert ={\left({\displaystyle {\int }_{\Omega }{\left|u\right|}^2\text{d}x\text{d}y}\right)}^{\frac{1}{2}}$

定义Sobolev空间：

$H_p^1\left(\Omega \right)=\left\{u\in H^1\left(\Omega \right):u\left(x,y\right)=u\left(x+2\pi ,y\right),u\left(x,y\right)=u\left(x,y+2\pi \right)\right\},$

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

${\left(u,v\right)}_{1,\Omega }={\displaystyle \underset{\left|\alpha \right|=0}{\overset{1}{\sum }}{\displaystyle {\int }_{\Omega }{D}^{\alpha }u{D}^{\alpha }\bar{v}\text{d}x\text{d}y}},$

${\Vert u\Vert }_{1,\Omega }={\left({\displaystyle \underset{\left|\alpha \right|=0}{\overset{1}{\sum }}{\Vert D^{\alpha }u\Vert }^2}\right)}^{\frac{1}{2}},$

其中 $D^{\alpha }=\frac{{\partial }^{\left|\alpha \right|}}{\partial x^{{\alpha }_1}\partial y^{{\alpha }_2}}$，$\alpha =\left({\alpha }_1,{\alpha }_2\right)$，$\left|\alpha \right|={\alpha }_1+{\alpha }_2$。

由格林公式及周期边界条件可知，问题(2.1)~(2.2)的弱形式为：找 $u\in H_p^1\left(\Omega \right)$，使得

$a\left(u,v\right)=F\left(v\right),\forall v\in H_p^1\left(\Omega \right),$ (2.3)

其中

$a\left(u,v\right)={\displaystyle {\int }_{\Omega }\nabla u\nabla \bar{v}\text{d}x\text{d}y}+{\displaystyle {\int }_{\Omega }Vu\bar{v}\text{d}x\text{d}y}+{\displaystyle {\int }_{\Omega }h\left(u^2\right)u\bar{v}\text{d}x\text{d}y},$

$F\left(v\right)={\displaystyle {\int }_{\Omega }f\bar{v}\text{d}x\text{d}y}.$

定义逼近空间：

$X_M\left(\Omega \right)=\text{Span}\left\{{\text{e}}^{itx}{\text{e}}^{iqy}:\left|t\right|=0,1,\cdots ,M;\left|q\right|=0,1,\cdots ,M\right\}.$

则弱形式(2.3)相应的离散格式为：找 $u_M\in X_M\left(\Omega \right)$，使得

$a\left(u_M,v_M\right)=F\left(v_M\right),\forall v_M\in X_M\left(\Omega \right).$ (2.4)

3. 算法的实现及其程序设计

3.1. 算法的实现

在这一节，我们将详细描述算法的实现过程，并给出该算法相应的程序设计。由于离散格式(2.4)是非线性的，我们将通过Picard迭代法进行求解。我们将(2.4)相应的线性问题的解作为迭代初值，即：找 $u_M^0\in X_M\left(\Omega \right)$，使得

${\left(\nabla u_M^0,\nabla v_M\right)}_N+{\left(V{u}_M^0,v_M\right)}_N={\left(f,v_M\right)}_N,\forall v_M\in X_M\left(\Omega \right),$ (3.1)

其中 ${\left(u,v\right)}_N$ 表示u与v的离散内积，则可建立(2.4)的一种Picard迭代格式：找 $u_M^m\in X_M\left(\Omega \right)$，使得

${\left(\nabla u_M^m,\nabla v_M\right)}_N+{\left(V{u}_M^m,v_M\right)}_N+{\left(h\left({\left(u_M^{m-1}\right)}^2\right)u_M^m,v_M\right)}_N={\left(f,v_M\right)}_N,\forall v_M\in X_M\left(\Omega \right),$ (3.2)

其中 $u_M^0$ 由(3.1)求出。下面我们将分别建立(3.1)和(3.2)的矩阵形式，令

$u_M^0={\displaystyle \underset{\left|t\right|=0}{\overset{M}{\sum }}{\displaystyle \underset{\left|q\right|=0}{\overset{M}{\sum }}{u}_{tq}^0{\text{e}}^{itx}{\text{e}}^{iqy}}},$ (3.3)

$u_M^m={\displaystyle \underset{\left|t\right|=0}{\overset{M}{\sum }}{\displaystyle \underset{\left|q\right|=0}{\overset{M}{\sum }}{u}_{tq}^m{\text{e}}^{itx}{\text{e}}^{iqy}}}.$ (3.4)

$U^m=\left(\begin{array}{ccccc}{u}_{-M,-M}^m& \cdots & u_{-M,0}^m& \cdots & u_{-M,M}^m\\ ⋮& \ddots & ⋮& \ddots & ⋮\\ u_{0,-M}^m& \cdots & u_{0,0}^m& \cdots & u_{0,M}^m\\ ⋮& \ddots & ⋮& \ddots & ⋮\\ u_{M,-M}^m& \cdots & u_{M,0}^m& \cdots & u_{M,M}^m\end{array}\right).$

我们用 ${\bar{U}}^m$ 表示由 $U^m$ 的列构成的长度为 ${\left(2M+1\right)}^2$ 的列向量。令 ${\phi }_t\left(x\right)={\text{e}}^{itx}$，${\phi }_q\left(x\right)={\text{e}}^{iqx}$，用 ${\zeta }_{\mu },w_{\mu }\left(\mu =0,1,\cdots ,N-1\right)$ 分别表示傅里叶积分的高斯点和权，将(3.3)代入(3.1)，并取 $v_M={\text{e}}^{ikx}{\text{e}}^{ily}$，$\left(\left|k\right|=0,1,\cdots ,M;\left|l\right|=0,1,\cdots ,M\right)$，则有

$\begin{array}{c}{\left(\nabla u_M^0,\nabla v_M\right)}_N={\displaystyle \underset{\left|t\right|=0}{\overset{M}{\sum }}{\displaystyle \underset{\left|q\right|=0}{\overset{M}{\sum }}{u}_{tq}^0}}{\left(\nabla \left({\text{e}}^{itx}{\text{e}}^{iqy}\right),\nabla \left({\text{e}}^{-ikx}{\text{e}}^{-ily}\right)\right)}_N\\ ={\displaystyle \underset{\left|t\right|=0}{\overset{M}{\sum }}{\displaystyle \underset{\left|q\right|=0}{\overset{M}{\sum }}{u}_{tq}^0}}{\displaystyle \underset{\mu ,\sigma =0}{\overset{N-1}{\sum }}{\partial }_x{\phi }_t}\left({\zeta }_{\mu }\right){\phi }_q\left({\zeta }_{\sigma }\right){\partial }_x{\phi }_k\left({\zeta }_{\mu }\right){\phi }_l\left({\zeta }_{\sigma }\right)w_{\mu }{w}_{\sigma }\\ +{\displaystyle \underset{\left|t\right|=0}{\overset{M}{\sum }}{\displaystyle \underset{\left|q\right|=0}{\overset{M}{\sum }}{u}_{tq}^0}}{\displaystyle \underset{\mu ,\sigma =0}{\overset{N-1}{\sum }}{\phi }_t}\left({\zeta }_{\mu }\right){\partial }_y{\phi }_q\left({\zeta }_{\sigma }\right){\phi }_k\left({\zeta }_{\mu }\right){\partial }_y{\phi }_l\left({\zeta }_{\sigma }\right)w_{\mu }{w}_{\sigma },\end{array}$

$\begin{array}{c}{\left(V{u}_M^0,v_M\right)}_N={\displaystyle \underset{\left|t\right|=0}{\overset{M}{\sum }}{\displaystyle \underset{\left|q\right|=0}{\overset{M}{\sum }}{u}_{tq}^0}}{\left(V{\text{e}}^{itx}{\text{e}}^{iqy},{\text{e}}^{-ikx}{\text{e}}^{-ily}\right)}_N\\ ={\displaystyle \underset{\left|t\right|=0}{\overset{M}{\sum }}{\displaystyle \underset{\left|q\right|=0}{\overset{M}{\sum }}{u}_{tq}^0}}{\displaystyle \underset{\mu ,\sigma =0}{\overset{N-1}{\sum }}V}\left({\zeta }_{\mu },{\zeta }_{\sigma }\right){\phi }_t\left({\zeta }_{\mu }\right){\phi }_q\left({\zeta }_{\sigma }\right){\phi }_k\left({\zeta }_{\mu }\right){\phi }_l\left({\zeta }_{\sigma }\right)w_{\mu }{w}_{\sigma },\end{array}$

${\left(f,v_M\right)}_N={\displaystyle \underset{\mu ,\sigma =0}{\overset{N-1}{\sum }}f}\left({\zeta }_{\mu },{\zeta }_{\sigma }\right){\phi }_k\left({\zeta }_{\mu }\right){\phi }_l\left({\zeta }_{\sigma }\right)w_{\mu }{w}_{\sigma },$

则(2.3)可以写为下列的矩阵形式

$\left(A+B+C\right){\bar{U}}^0=F,$ (3.5)

其中

$a_{tqkl}={\displaystyle \underset{\mu ,\sigma =0}{\overset{N-1}{\sum }}{\partial }_x{\phi }_t}\left({\zeta }_{\mu }\right){\phi }_q\left({\zeta }_{\sigma }\right){\partial }_x{\phi }_k\left({\zeta }_{\mu }\right){\phi }_l\left({\zeta }_{\sigma }\right)w_{\mu }{w}_{\sigma },$

$b_{tqkl}={\displaystyle \underset{\mu ,\sigma =0}{\overset{N-1}{\sum }}{\phi }_t}\left({\zeta }_{\mu }\right){\partial }_y{\phi }_q\left({\zeta }_{\sigma }\right){\phi }_k\left({\zeta }_{\mu }\right){\partial }_y{\phi }_l\left({\zeta }_{\sigma }\right)w_{\mu }{w}_{\sigma },$

$c_{tqkl}={\displaystyle \underset{\mu ,\sigma =0}{\overset{N-1}{\sum }}V}\left({\zeta }_{\mu },{\zeta }_{\sigma }\right){\phi }_t\left({\zeta }_{\mu }\right){\phi }_q({\zeta }_{\sigma }){\phi }_k\left({\zeta }_{\mu }\right){\phi }_l\left({\zeta }_{\sigma }\right)w_{\mu }{w}_{\sigma },$

$f_{kl}={\displaystyle \underset{\mu ,\sigma =0}{\overset{N-1}{\sum }}f}\left({\zeta }_{\mu },{\zeta }_{\sigma }\right){\phi }_k\left({\zeta }_{\mu }\right){\phi }_l\left({\zeta }_{\sigma }\right)w_{\mu }{w}_{\sigma }$

$A={\left(a_{tqkl}\right)}_{t,q,k,l=0}^{2M+1},B={\left(b_{tqkl}\right)}_{t,q,k,l=0}^{2M+1},C={\left(c_{tqkl}\right)}_{t,q,k,l=0}^{2M+1},F={\left(f_{kl}\right)}_{k,l=0}^{2M+1}.$

类似地，将(3.4)代入(3.2)可得

$\begin{array}{c}{\left(\nabla u_M^m,\nabla v_M\right)}_N={\displaystyle \underset{\left|t\right|=0}{\overset{M}{\sum }}{\displaystyle \underset{\left|q\right|=0}{\overset{M}{\sum }}{u}_{tq}^m}}{\left(\nabla \left({\text{e}}^{itx}{\text{e}}^{iqy}\right),\nabla \left({\text{e}}^{-ikx}{\text{e}}^{-ily}\right)\right)}_N\\ ={\displaystyle \underset{\left|t\right|=0}{\overset{M}{\sum }}{\displaystyle \underset{\left|q\right|=0}{\overset{M}{\sum }}{u}_{tq}^m}}{\displaystyle \underset{\mu ,\sigma =0}{\overset{N-1}{\sum }}{\partial }_x{\phi }_t}\left({\zeta }_{\mu }\right){\phi }_q\left({\zeta }_{\sigma }\right){\partial }_x{\phi }_k\left({\zeta }_{\mu }\right){\phi }_l\left({\zeta }_{\sigma }\right)w_{\mu }{w}_{\sigma }\\ +{\displaystyle \underset{\left|t\right|=0}{\overset{M}{\sum }}{\displaystyle \underset{\left|q\right|=0}{\overset{M}{\sum }}{u}_{tq}^m}}{\displaystyle \underset{\mu ,\sigma =0}{\overset{N-1}{\sum }}{\phi }_t}\left({\zeta }_{\mu }\right){\partial }_y{\phi }_q\left({\zeta }_{\sigma }\right){\phi }_k\left({\zeta }_{\mu }\right){\partial }_y{\phi }_l\left({\zeta }_{\sigma }\right)w_{\mu }{w}_{\sigma },\end{array}$

$\begin{array}{c}{\left(V{u}_M^m,v_M\right)}_N={\displaystyle \underset{\left|t\right|=0}{\overset{M}{\sum }}{\displaystyle \underset{\left|q\right|=0}{\overset{M}{\sum }}{u}_{tq}^m}}{\left(V{\text{e}}^{itx}{\text{e}}^{iqy},{\text{e}}^{-ikx}{\text{e}}^{-ily}\right)}_N\\ ={\displaystyle \underset{\left|t\right|=0}{\overset{M}{\sum }}{\displaystyle \underset{\left|q\right|=0}{\overset{M}{\sum }}{u}_{tq}^m}}{\displaystyle \underset{\mu ,\sigma =0}{\overset{N-1}{\sum }}V}\left({\zeta }_{\mu },{\zeta }_{\sigma }\right){\phi }_t\left({\zeta }_{\mu }\right){\phi }_q\left({\zeta }_{\sigma }\right){\phi }_k\left({\zeta }_{\mu }\right){\phi }_l\left({\zeta }_{\sigma }\right)w_{\mu }{w}_{\sigma },\end{array}$

$\begin{array}{c}{\left(h\left({\left(u_M^{m-1}\right)}^2\right)u_M^m,v_M\right)}_N={\displaystyle \underset{\left|t\right|=0}{\overset{M}{\sum }}{\displaystyle \underset{\left|q\right|=0}{\overset{M}{\sum }}{u}_{tq}^m}}{\left(h\left({\left(u_M^{m-1}\right)}^2\right){\text{e}}^{itx}{\text{e}}^{iqy},{\text{e}}^{-ikx}{\text{e}}^{-ily}\right)}_N\\ ={\displaystyle \underset{\left|t\right|=0}{\overset{M}{\sum }}{\displaystyle \underset{\left|q\right|=0}{\overset{M}{\sum }}{u}_{tq}^m}}{\displaystyle \underset{\mu ,\sigma =0}{\overset{N-1}{\sum }}h}\left({\left(u_M^{m-1}\left({\zeta }_{\mu },{\zeta }_{\sigma }\right)\right)}^2\right){\phi }_t\left({\zeta }_{\mu }\right){\phi }_q\left({\zeta }_{\sigma }\right){\phi }_k\left({\zeta }_{\mu }\right){\phi }_l\left({\zeta }_{\sigma }\right)w_{\mu }{w}_{\sigma },\end{array}$

则(2.4)可以写为下列的矩阵形式

$\left(A+B+C+D_m\right){\bar{U}}^m=F,$ (3.6)

其中

$d_{tqkl}={\displaystyle \underset{\mu ,\sigma =0}{\overset{N-1}{\sum }}h}\left({\left(u_M^{m-1}\left({\zeta }_{\mu },{\zeta }_{\sigma }\right)\right)}^2\right){\phi }_t\left({\zeta }_{\mu }\right){\phi }_q\left({\zeta }_{\sigma }\right){\phi }_k\left({\zeta }_{\mu }\right){\phi }_l\left({\zeta }_{\sigma }\right)w_{\mu }{w}_{\sigma },D_m={\left(d_{tqkl}\right)}_{t,q,k,l=0}^{2M+1}.$

3.2. 算法的程序设计

4. 数值实验

为了表明算法的有效性，我们将进行数值实验。我们在MATLAB2018b平台上编程计算。定义逼近解与精确解之间的误差如下：

$e\left(u\left(x,y\right),u_M^K\left(x,y\right)\right)={\Vert u\left(x,y\right)-u_M^K\left(x,y\right)\Vert }_{L^2\left(\Omega \right)},$

其中K表示迭代次数。

例1 我们取 $u={\text{e}}^{\cos \left(x+y\right)}$，$V\left(x,y\right)=\frac{1}{2}\left(x^2+y^2\right)$，$h\left(u^2\right)={\left|u\right|}^2$，显然u满足周期边界条件，将 $u,V,h(\; u\; 2\; )$

代入(2.3)可算出f。我们取迭代次数 $K=50$，对于不同的M和N，我们在表1中列出了逼近解与精确解之间的误差结果。

表1. 当K = 50时，对于不同的M和N，逼近解与精确解之间的误差结果

为了进一步表明算法的收敛性和高精度，我们还在图1中分别画出了精确解和逼近解的图像，在图2中分别画出了精确解与 $N=30$，$M=10$ 和 $N=40$，$M=20$ 时的逼近解之间的误差图像。

图1. 精确解(左)与N = 40和M = 20时的逼近解(右)的图像

图2. 精确解与N = 30, M = 10 (左)和N = 40, M = 20 (右)时的逼近解之间的误差图像

从表1中可以看出，当迭代次数 $K=50$，$N\ge 35$，$M\ge 15$ 时，逼近解 $u_M^K\left(x,y\right)$ 达到了大约10⁻¹³的精度。另外，从图1，图2进一步观察到我们的算法是收敛的和高精度的。

文章引用

江婷婷. 变系数非线性二阶问题有效的Fourier谱逼近
Efficient Fourier Spectral Approximation for Nonlinear Second-Order Problems with Variable Coefficients[J]. 应用数学进展, 2022, 11(07): 4268-4277. https://doi.org/10.12677/AAM.2022.117453

参考文献