二维薛定谔方程的全离散有限元两层网格方法

王建云，田智鲲，张 丹

（1.湖南工业大学 理学院，湖南 株洲 412007；2.湖南工程学院 理学院，湖南 湘潭 411104）

1 研究背景

量子力学是描述微观世界结构、运动与变化规律的物理学科。它的发现与建立是20世纪人类文明发展的一个重大飞跃，而且引发了一系列划时代的科学发现与技术发明。

薛定谔（Schrödinger）方程是量子力学最基本的方程，由奥地利物理学家薛定谔于1926年提出，揭示了原子世界中物质运动的基本规律。例如，最简单的一维无界域上依赖于时间的薛定谔方程，就反映了一个在守恒场中无旋转的粒子的运动状态，而其解的平方就表示该粒子在给定时刻出现在给定位置上的概率密度。许多物理问题，如高能物理、非线性媒体中的激光束扫描、浓缩问题等的数学模型也都可归结为（非线性）薛定谔方程。

薛定谔方程也是量子力学的一个基本假定，不能用比它更根本的假定来证明，其正确性只能靠实践来检验，而且在实际中复杂的薛定谔方程不易求得其精确解。因此，关于其数值解的研究越来越受到人们的重视和关注。

本文主要考虑如下二维依赖于时间的线性薛定谔方程的初边值问题：

有限元方法是求解偏微分方程数值解的一种重要方法，它的基本思想是根据变分原理，利用有限元空间上的离散解近似无穷维空间上的连续解。目前，许多数值求解方法被应用到线性或非线性薛定谔方程中，其中有限元方法为常见的方法之一，具体参见文献[1-10]。

文献[1]利用Galerkin有限元方法将空间离散，用Crank-Nicolson方法将时间方向离散得到有限元全离散格式，证明了数值解在全离散格式下的存在性和唯一性，并给出了在时间和空间步长满足一定条件下的最优阶L2误差估计。文献[2]利用人工边界条件，构造和分析了无界域二维依赖于时间的薛定谔方程的一种全离散格式，并证明了该格式具有无条件稳定和收敛。文献[7]讨论了二维薛定谔方程的有限元超收敛问题，利用椭圆投影算子得到数值解的H1误差具有超收敛性，并运用插值后处理技术得到整体超收敛结果。

两层网格方法是求解偏微分方程的一种高效数值计算方法，其基本思想是通过构造两种不同尺度（基于粗网格和细网格）的有限元空间。首先在粗网格上求解原来的非对称或非线性问题；然后利用粗网格上的数值解，将原问题用合适的方式进行对称化或线性化；再在细网格上求解相应的对称化或线性化问题。与直接在细网格上数值求解原问题相比，两层网格方法减少了计算量，大幅提高了算法效率。该方法最初由许进超提出，他在文献[11-13]中，通过求解非对称不定以及非线性椭圆问题，引入粗细两个子空间进行离散，构造了椭圆型问题有限元的一系列两层网格算法。几乎在同一时期，黄云清等在文献[14]中，研究了非线性奇异两点边值问题的多层迭代校正算法的收敛性误差估计与逼近解的渐近展式。金继承等在文献[15]中，首次将两层网格方法运用到求解一类耦合的偏微分方程组，构造了解耦的有限元两层网格算法，并进行了误差估计。

后来两层网格有限元方法被应用到线性和非线性薛定谔方程，具体见文献[16-19]。文献[18]针对一类二维依赖于时间的线性薛定谔方程，其空间方向是在三角形网格上采用线性元进行离散，而时间方向采用向后欧拉方法进行离散，得到一种全离散的两层网格有限元算法，并且进行了理论分析。本文基于已有研究，主要考虑空间方向在矩形网格上利用双线性元进行离散，并通过数值算例验证了两层网格方法的高效性，且数值解在相同网格尺寸下比文献[18]中数值解的误差更小。

2 全离散有限元格式

记Lp(Ω)为标准的Banach空间，具有范数，Wm,p为定义在区域Ω上的标准Sobolev空间，其范数定义为，并且当p=2时，记Wm,p=Hm，相应的范数简记为m=m,2，=0.2。

对应的L2范数为

将区域Ω拟一致剖分为矩形网格Γh，其中网格步长0

记τ=T/N为区间[0,T]上的时间步长，tj=jτ，j=0,1,…,N为区间上的时间节点，其中N为某个正整数。

为简单起见，将函数ω(x,y,tj)简记为ωj。另外，对任意函数序列ωj(x,y),j=1,2…，利用向后欧拉方法定义其关于时间的差商，为

那么，问题（1）的全离散有限元解

可定义为满足

式中Phu0(x,y)为u0(x,y)在Sh空间的椭圆投影[18]。

因此，类似文献[18]，可得如下定理1。

定理1设u(x,y,t)为问题（1）的精确解，函数序列Uhn(x,y)为满足问题（2）的全离散有限元解，则有误差估计

其中C为与h无关的任意常数。

在实际编程计算时，需要将问题（2）按实部和虚部分开，则可以写成如下等价的耦合方程组：

式 中Phu0R(x,y)、Phu0I(x,y)分 别为u0R(x,y)、u0I(x,y)在Sh空间的椭圆投影。

3 两层网格算法

将Ω拟一致剖分为矩形网格ΓH，其网格尺寸H>>h，而SH为定义在ΓH上的另一个有限元空间。利用向后Euler方法在时间方向离散，构造问题（1）的一种全离散两网格有限元算法，先在粗网格上对原问题进行求解，然后在细网格上计算时利用粗网格上已算得的部分数值作为已知数据代入，将薛定谔方程耦合的实部和虚部进行解耦，从而简化为在细网格上求解两个泊松方程。这样，由于粗网格尺寸H>>h，因此该算法比直接在细网格上求解原问题要简便很多，将大大减少数值计算的时间。

算法具体步骤如下。

步骤1求函数序列满足

步骤2求函数序列满足

因此，类似文献[18]，可得如下定理2。

定理2设u(x,y,t)为问题（1）的精确解，函数序列uhn(x,y)为两层网格算法中满足问题（5）的两层网格有限元解，则有误差估计

其中C为与h和H无关的任意常数。

4 数值实验

考虑如下二维依赖于时间的线性薛定谔方程

式中：Ω=[0,1]×[0,1]；T=1；右端项函数f(x,y,t)选择满足如下精确解[18]

采用空间步长为h的均匀网格，对区域Ω进行拟一致矩形网格剖分Γh，利用双线性有限元进行数值求解。取时间步长为τ=10-3，分别取细网格步长h=1/4,1/16,1/64，而ΓH为粗网格步长的矩形剖分。Uhn为细网格Γh上，时间方向利用向后欧拉格式，空间方向利用有限元方法，得到的标准有限元解；为粗网格ΓH和细网格Γh上，利用两层网格算法得到的两层网格有限元解。分别计算精确解与标准有限元解的误差及与两层网格有限元解的误差。不同时刻的误差结果及CPU运行时间，如表1～3所示。

表1 数值解在t=0.2时的误差及CPU运行时间Table1 Errors of numerical solutions and CPU running time with t=0.2

表2 数值解在t=0.5时的误差及CPU运行时间Table2 Errors of numerical solutions and CPU running time with t=0.5

表3 数值解在t=1.0时的误差及CPU运行时间Table3 Errors of numerical solutions and CPU running time with t=1.0

图1和图2分别给出了本文与文献[18]中的两层网格有限元解在不同时刻的误差比较，从图中可以看出，本文的两层网格有限元解在每一时刻都要更加精确。

图1 16×16网格上两层网格解在不同时刻的误差Fig.1 Errors of two-grid solutions at different time levels on 16×16 meshes

图2 64×64网格上两层网格解在不同时刻的误差Fig.2 Errors of two-grid solutions at different time levels on 64×64 meshes

5 结语

本文研究了两层网格方法在二维线性薛定谔方程中的应用，对空间区域进行拟一致矩形网格剖分，构建了一种向后欧拉全离散有限元格式的两层网格算法，解释了两层网格算法在H1求解线性薛定谔方程中的思想和原理，给出了数值解在范数下的误差阶，最后通过一个数值算例验证了该算法的可行性和高效性。