有限差分–配点法求解SchrO^¨dinger方程

Finite Difference-Collocation Method for the SchrO^¨dinger Equation

摘要: 本文将重心插值配点法结合Crank-Nicolson差分格式求解线性Schrödinger方程。首先，对方程中的空间方向上采用重心插值Chebyshev配点格式进行离散，时间方向采用Crank-Nicolson差分格式，从而导出对应的代数方程组。最后，数值算例验证该计算格式具有高精度性，并且满足质量和能量守恒性。

Abstract: In this paper, the barycentric interpolation collocation method combined with Crank-Nicolson dif-ference scheme is used to solve the linear Schrödinger equation. Firstly, the space direction is dis-cretized by the Chebyshev collocation method of barycentric interpolation, and the time direction is discretized by the Crank-Nicolson difference scheme. Finally, numerical examples show that the numerical scheme has high precision and satisfies the conservation of mass and energy.

1. 引言

1926年，奥地利物理学家薛定谔提出了著名的Schrödinger方程 [1]，用以解决微观世界的不确定性，揭示了微观物理世界物质运动的基本规律，就像牛顿定律在经典力学中所起的作用一样，是原子物理学中处理一切非相对论问题的有力工具。Schrödinger方程的应用非常广泛，在原子、分子、固体物理、核物理、化学等领域中都大放异彩。可是Schrödinger方程通常难以找到解析解，必须通过数值解来分析。因此在过去的几十年里，对Schrödinger方程的数值计算研究已经引起了学者们的广泛关注。

近几年来，专家学者们针对Schrödinger方程的数值解进行了系统性的研究，例如：文献 [2] 主要介绍了通过分步傅里叶变换法对Schrödinger方程进行数值仿真。文献 [3] 提出了带有五次项的差分格式，并对数值格式的理论进行分析。文献 [4] 提出了一类非自共轭的Schrödinger方程的三层的差分格式，并证明了差分格式解的存在唯一性以及先验估计，最后数值实验证明了其差分格式在空间和时间上都达到了二阶收敛精度。在文献 [5] 中主要对带有波动算子的Schrödinger方程进行了研究，并提出几类守恒的差分格式。文献 [6] 首次提出了局部间断有限元法求解Schrödinger方程，并且证明了该方法的稳定性。文献 [7] 用有限差分显格式求解Schrödinger方程，并给出最大范数意义下的误差分析。文献 [8] 提出用有限差分显格式结合人工边界条件求解无界区域的Schrödinger方程，并给出算法的稳定性和收敛性分析。文献 [9] 用伪谱傅里叶方法结合分裂步求解光纤中的Schrödinger方程。文献 [10] 使用了一种hr自适应方法来求解带有三次Schrödinger方程。文献 [11] 使用了三次B样条Galerkin方法求解耦合Schrödinger方程的数值解。文献 [12] 利用积分因子法、对角阵技巧结合FFT算法求解高维带有阻尼薛定谔方程。

相较于显格式的有条件稳定，隐格式往往具有无条件稳定的特点，例如Crank-Nicolson格式。Patel等 [13] 对于偏微分方程提出一种无条件稳定的紧致有限差分格式。文献 [14] 结合不同的差分格式对Schrödinger方程进行数值求解。文献 [15] 用四阶精度差分法解定态薛定谔方程，并用于求解几个常见势场中的定态方程。利用隐式格式在进行数值模拟时，可以取较大的时间步长，但是每一时间层需要求解线性或者非线性的代数方程组。

上述方法都是依赖网格剖分来求解Schrödinger方程，本文利用一种新型的无网格计算方法——重心插值配点格式 [16] 来求解这类问题，通过选取Chebyshev节点能有效克服数值不稳定现象，具有格式简单、精度高、程序简便、节点适应性好等特点。不同于以单元为基础的传统数值方法，重心插值配点格式是无网格配点法的一种，其没有对计算区域进行网格划分，而是采用节点的近似来构造近似函数，能有效避免网格重构带来的累积误差，其计算结果的精度更高。最近，文献 [17] 利用重心插值配点格式求解微分方程初边值问题，很多研究学者将该方法推广到求解各类微分方程中，比如平面弹性问题 [18]、Fredholm积分方程 [19]、热传导方程 [20]、Allen-Cahn方程 [21] [22] [23] [24]、二维时间分数阶电报方程 [25]、Burgers方程 [26]、Cahn-Hilliard方程 [27] 和Black-Scholes方程 [28] 等微分方程。如今，重心插值配点格式广泛地被应用于弹性力学、微波技术及流体力学等多个方面。

因此，本文希望在上述工作的研究基础上，主要利用重心插值配点法设计线性Schrödinger方程的高精度数值算法格式，通过数值算例验证格式的高精度和有效性。与前人工作比较，我们的格式用较少的节点个数就可以达到高精度，以进一步推动求解非线性Schrödinger方程的高精度数值算法的发展。

2. 线性Schrödinger方程的重心插值配点格式

2.1. 线性Schrödinger方程

线性Schrödinger方程的初值、边值问题考虑如下：

$\{\begin{array}{l}i\frac{\partial u\left(x,t\right)}{\partial t}+\Delta u\left(x,t\right)+v\left(x\right)u\left(x,t\right)=0,\left(x,t\right)\in \Omega \times \left(0,T\right]\\ u\left(x,0\right)=u_0\left(x\right),x\in \Omega \end{array}$ (1)

边值条件为

$u\left(x,t\right)=\varphi \left(x,t\right),x\in \partial \Omega ,t\in \left(0,T\right]$ (2)

其中i是一个虚数，满足 $i^2=-1$， $\Omega$ 是 $R^d\left(d=1,2\right)$ 的有界域， $u\left(x,t\right)$ 是需要我们求解的未知波函数， $v\left(x\right)$ 是势函数。

具有狄利克雷边界的Schrödinger方程保持了很多物理量。本文也将会验证格式的质量守恒和能量守恒：

$M\left(t\right)={\displaystyle \underset{\Omega }{\iint }{\left|u\left(x,t\right)\right|}^2\text{d}\Omega }=M\left(0\right),t\in R^{+}$ (3)

以及

$E\left(t\right)=-\frac{1}{2}{\displaystyle \underset{\Omega }{\iint }{\left|\nabla u\left(x,t\right)\right|}^2\text{d}\Omega }+\frac{v\left(x\right)}{2}{\displaystyle \underset{\Omega }{\iint }{\left|u\left(x,t\right)\right|}^2\text{d}\Omega }=E\left(0\right),t\in R^{+}$ (4)

其中， $M\left(t\right)$ 和 $E\left(t\right)$ 分别代表不同时间的质量值和能量值。

2.2. 重心Lagrange插值

设有插值节点 $x_j$ 和一组对应的实数 $f_j\left(j=0,1,\cdots ,m\right)$,采用多项式插值，则在次数不超过m的多项式空间可以找到唯一的插值多项式 $p\left(x\right)$，满足 $p\left(x_j\right)=f_j,j=0,1,\cdots ,m$，写成Lagrange插值公式，即

$p\left(x\right)={\displaystyle \underset{j=0}{\overset{m}{\sum }}{\gamma }_j}\left(x\right)f_j$ (5)

其中 ${\gamma }_j\left(x\right)={\displaystyle \underset{i=0,i\ne j}{\overset{m}{\prod }}\left(x-x_i\right)}/{\displaystyle \underset{i=0,i\ne j}{\overset{m}{\prod }}\left(x_j-x_i\right)},j=0,1,\cdots ,m$，为Lagrange插值基函数，满足性质

${\gamma }_j\left(x_i\right)={\delta }_{ji}=\{\begin{array}{l}1,j=i,\\ 0,j\ne i,\end{array}\left(i,j=0,1,\cdots ,m\right)$ (6)

以及

${\displaystyle \underset{j=0}{\overset{m}{\sum }}{\gamma }_j}\left(x\right)=1$ (7)

令 $l\left(x\right)=\left(x-x_0\right)\left(x-x_1\right)\cdots \left(x-x_m\right)$，定义重心权： $w_j=\frac{1}{{\displaystyle \underset{i=0,i\ne j}{\overset{m}{\prod }}\left(x_j-x_i\right)}},j=0,1,\cdots ,m$

则插值基函数可以表示为

${\gamma }_j\left(x\right)=l\left(x\right)\frac{w_j}{x-x_j},j=0,1,\cdots ,m$ (8)

将(8)代入(5)，得

$p\left(x\right)=l\left(x\right){\displaystyle \underset{j=0}{\overset{m}{\sum }}\frac{w_j}{x-x_j}}{f}_j$ (9)

由(7)~(9)可得重心Lagrange插值公式

$p\left(x\right)=\frac{{\displaystyle \underset{j=0}{\overset{m}{\sum }}\frac{w_j}{x-x_j}}{f}_j}{{\displaystyle \underset{j=0}{\overset{m}{\sum }}\frac{w_j}{x-x_j}}}$ (10)

重心Lagrange插值对于等距节点是病态的，通过选用服从密度比为 ${\left(1-x^2\right)}^{-\frac{1}{2}}$ 的节点族，可使其具有非常好的数值稳定性。满足这样条件的最简单的节点分布是Chebyshev点族，因此本文采用第二类Chebyshev点族： $x_j=\cos \left(\frac{j}{m}\pi \right),j=0,1,\cdots ,m$，保证重心Lagrange插值的数值稳定性。

2.3. Schrödinger方程的半离散数值解格式

下面对一维Schrödinger方程的半离散格式进行推导，二维Schrödinger方程的推导类似。

对于Schrödinger方程的空间区域 $\left[a,b\right]$ 进行离散成为 $m+\text{1}$ 个第二类Chebyshev节点，即

$a=x_{\text{0}}<x_{\text{1}}<\cdots <x_m=b$

则区域上的 $m+\text{1}$ 个Chebyshev节点为 $\left(x_i,t\right),i=\text{0},\text{1},\cdots ,m$，将未知函数 $u\left(x,t\right)$ 在空间域节点 $x_{\text{0}},x_{\text{1}},\cdots ,x_m$ 上的值记作

$u\left(x_i,t\right)=u_i\left(t\right),i=\text{0},\text{1},\cdots ,m$ (11)

然后，将未知函数 $u\left(x,t\right)$ 表示为在节点 $x_{\text{0}},x_{\text{1}},\cdots ,x_m$ 上的重心型插值近似函数，即

$u\left(x,t\right)={\displaystyle \underset{j=\text{0}}{\overset{m}{\sum }}{L}_j\left(x\right)u_j\left(t\right)}$ (12)

其中， $L_j\left(x\right)$ 为此重心型插值的插值基函数。将式(12)代入方程，并且让方程(12)在节点 $x_{\text{0}},x_{\text{1}},\cdots ,x_m$ 上成立，可得常微分方程组，即

$i\left({\displaystyle \underset{j=\text{0}}{\overset{m}{\sum }}{L}_j\left(x_i\right)\frac{\text{d}{u}_j\left(t\right)}{\text{d}t}}\right)+{\displaystyle \underset{j=\text{0}}{\overset{m}{\sum }}{L^{″}}_j\left(x_i\right)u_j\left(t\right)}+V\left({\displaystyle \underset{j=\text{0}}{\overset{m}{\sum }}{L}_j\left(x_i\right)u_j\left(t\right)}\right),i=\text{0},\text{1},\cdots ,m$ (13)

注意到 $L_j\left(x_i\right)={\delta }_{ij},{L^{″}}_j\left(x_i\right)=\frac{{\text{d}}^2{L}_j\left(x_i\right)}{\text{d}{x}^2}=C_{ij}^{\left(2\right)},V\left(u\right)=v\left(x\right)u$

其中， ${\delta }_{ij}=\{\begin{array}{l}1,i=j\\ 0,i\ne j\end{array}$， $C_{ij}^{\left(p\right)},p=1,2$。为节点 $x_0,x_1,\cdots ,x_m$ 上重心型插值p阶微分矩阵的元素。则方程(13)可以记作下面的矩阵形式，即

$i\left(\begin{array}{c}\frac{\text{d}{u}_0\left(t\right)}{\text{d}t}\\ ⋮\\ \frac{\text{d}{u}_m\left(t\right)}{\text{d}t}\end{array}\right)+\left(\begin{array}{ccc}{C}_{\text{00}}^{\left(2\right)}& \cdots & C_{\text{0}m}^{\left(2\right)}\\ ⋮& \ddots & ⋮\\ C_{m\text{0}}^{\left(2\right)}& \cdots & C_{mm}^{\left(2\right)}\end{array}\right)\left(\begin{array}{c}{u}_{\text{0}}\left(t\right)\\ ⋮\\ u_m\left(t\right)\end{array}\right)+\left(\begin{array}{c}V\left(u_{\text{0}}\left(x,t\right)\right)\\ ⋮\\ V\left(u_m\left(x,t\right)\right)\end{array}\right)=0$ (14)

若令 $U={\left[u_0\left(t\right),u_1\left(t\right),\cdots ,u_m\left(t\right)\right]}^{\text{T}}$

则方程(14)可以改写为以下矩阵形式：

$i\frac{\text{d}U}{\text{d}t}+\left(C^{\left(2\right)}\otimes I_m\right)U+V\left(U\right)=0$ (15)

上式中：符号 $\otimes$ 表示矩阵的Kronecker积， $C^{\left(2\right)}$ 为关于节点 $x_0,x_1,\cdots ,x_m$ 的重心型插值2阶微分矩阵； $I_m$ 分别为m阶单位矩阵。

类似一维问题的推导，二维Schrödinger方程的半离散格式写成矩阵形式如下：

$i\frac{\text{d}U}{\text{d}t}+\left(C^{\left(2\right)}\otimes I_m\right)U+\left(I_n\otimes D^{\left(2\right)}\right)U+V\left(U\right)=0$ (16)

其中 $C^{\left(2\right)}$ 为关于节点 $x_0,x_1,\cdots ,x_m$ 的重心型插值2阶微分矩阵， $D^{\left(2\right)}$ 为关于节点 $y_0,y_1,\cdots ,y_n$ 的重心型插值2阶微分矩阵。

2.4. Schrödinger方程的有限差分配点全离散格式

为了得到非线性Schrödinger方程的全离散格式，进一步对时间方向上用等距剖分成 $n+1$ 个节点，故节点为 $0=t_{\text{0}}<t_{\text{1}}<\cdots <t_n=T$。

将未知函数 $u_i\left(t\right)$ 在时间域节点 $t_{\text{0}},t_1,\cdots ,t_n$ 上的值记作 $u_i\left(t_k\right)=u\left(x_i,t_k\right)=u_i^k,i=\text{0},1,\cdots ,m;k=0,1,\cdots ,n$

在时间方向上采用Crank-Nicolson差分格式，则可以改写出以下矩阵形式：

$\frac{i}{\tau }\left(U_i^{k+1}-U_i^k\right)+\frac{1}{2}\left(C^{\left(2\right)}\otimes I_m\right)\left(U_i^{k+1}+U_i^k\right)+\frac{\text{1}}{2}\left(V_i^{k+1}+V_i^k\right)=0$ (17)

其中 $I_m$ 为m阶单位矩阵， $V\left(u\right)=v\left(x\right)u$ 从而得到Schrödinger方程的Crank-Nicolson差分–配点全离散数值格式。

那么，方程(17)可以化简为如下格式：

$\left[iI+\frac{\tau }{2}\left(C^{\left(2\right)}\otimes I_m\right)+\frac{\tau }{2}diag\left(v\left(x\right)\right)\right]U_i^{k+1}=F_i^k$ (18)

其中， $F^k=\left[iI-\frac{\tau }{2}\left(C^{\left(2\right)}\otimes I_m\right)-\frac{\tau }{2}diag\left(v\left(x\right)\right)\right]U_i^k$。

最后进行迭代求解得出相应数值解。

同样类似一维方程的推导，二维Schrödinger方程的全离散格式写成矩阵形式如下：

$\left[iI+\frac{\tau }{2}\left(C^{\left(2\right)}\otimes I_m\right)+\frac{\tau }{2}\left(I_n\otimes D^{\left(2\right)}\right)+\frac{\tau }{2}diag\left(v\left(x\right)\right)\right]U_i^{k+1}=F_i^k$ (19)

其中， $F^k=\left[iI-\frac{\tau }{2}\left(C^{\left(2\right)}\otimes I_m\right)-\frac{\tau }{2}\left(I_n\otimes D^{\left(2\right)}\right)-\frac{\tau }{2}diag\left(v\left(x\right)\right)\right]U_i^k$。

3. 数值算例

为验证利用重心插值配点格式求解Schrödinger方程的有效性和合理性，下面我们通过四个算例对一维、二维Schrödinger方程运用重心插值配点格式，对此格式的高精度和质量能量守恒进行验证。

为了方便，定义如下的最大绝对误差和相对误差符号：

最大绝对误差： $E_{\infty }=\underset{0\le i\le N}{{\displaystyle \max }}\left|U_i-u_i\right|$

最大相对误差： $E{r}_{\infty }=\frac{\underset{0\le i\le N}{{\displaystyle \max }}\left|U_i-u_i\right|}{\underset{0\le i\le N}{{\displaystyle \max }}\left|u_i\right|}$

其中， $u_i$ 表示点 $x_i$ 处的精确解， $U_i$ 表示点 $x_i$ 处的数值解。

3.1. 算例一

对于一维Schrödinger方程的边界条件为 $u\left(0,t\right)=u\left(1,t\right)=0,0\le t\le 1$ 时函数 $v\left(x\right)=0$。

初值条件为 $u\left(x,0\right)=\sin \left(\pi x\right),0\le x\le 1$ 此问题的带右端项的精确解为 $u\left(x,t\right)=\cos t\sin \left(\pi x\right),0\le x\le 1$。

表1表明，对于一维Schrödinger方程，推导得出其在时间上是二阶精度的。

Table 1. The time convergence order of one-dimensional Schrödinger equation

表1. 一维Schrödinger方程的时间收敛阶

固定时间节点，改变空间节点，表2给出了二阶有限差分格式和重心插值配点格式所得的最大绝对误差和最大相对误差，经过比较可以看出，重心插值配点格式在空间上可以用较少的点便能达到更高的精度。比如，当 $\tau ={10}^{-4}$ 时，重心Lagrange插值只需要在空间上取8个点就可以达到10⁻⁵的精度，但二阶有限差分格式需要取256个点才能达到这个精度。

下表中M + 1是为空间的节点数，N取10,000等分。

Table 2. Maximum absolute error and maximum relative error of second-order finite difference scheme and barycentric interpolation collocation scheme ( τ = 10 − 4 )

表2. 二阶有限差分格式和重心插值配点格式的最大绝对误差和最大相对误差( $\tau ={10}^{-4}$)

通过下面的图1和图2我们可以看出，重心插值配点格式的数值解和解析解的等势线图像基本一致，这进一步说明该数值格式的稳定性。

下图为M = 16，N = 1000时的图像：

Figure 1. Analytical solution of example one

图1. 算例一的解析解

Figure 2. Example one numerical solution of barycentric interpolation collocation scheme

图2. 算例一重心插值配点法格式数值解

为了进一步说明该格式的优越性，画出两种格式的绝对误差图。由图3和图4可知，当时间剖分一致时，重心插值配点格式在空间上取16个点的时候，绝对误差精度就达到10⁻⁹，而二阶有限差分格式在空间上取128个点的时候，精度只有10⁻⁵。

Figure 3. Example one absolute error figure of second-order finite difference scheme

图3. 算例一二阶有限差分格式绝对误差图

Figure 4. Example one absolute error figure of barycentric interpolation collocation scheme

图4. 算例一重心插值配点格式绝对误差图

3.2. 算例二

Schrödinger方程的边界条件为 $u\left(0,t\right)=u\left(1,t\right)=0,0\le t\le 1$ 时函数 $v\left(x\right)={\cos }^{2}x$。

初值条件为 $u\left(x,0\right)=\sin x,0\le x\le 2\pi ,0\le t\le 1$ 此问题的精确解 $u\left(x,t\right)=\exp \left(-\frac{3it}{2}\right)\sin x,\text{0}\le x\le 2\pi ,0\le t\le 1$。

进一步针对Schrödinger方程的离散质量和能量守恒性进行验证。选取空间步长为 $\tau ={10}^{-4}$，对应的计算区间为 $\left[0,1\right]\times \left[0,t\right]$，实验结果以图形的方式给出。图5给出重心Lagrange插值格式对质量守恒量的保持情况，图6给出重心Lagrange插值格式对能量守恒量的保持情况。从图像可以看出我们的配点格式很好地满足离散质量和能量守恒性。

(a) 总质量 (b) 质量相对误差

Figure 5. Example two mass conservation

图5. 算例二质量守恒

(a) 总能量 (b) 能量相对误差

Figure 6. Example two energy conservation

图6. 算例二能量守恒

3.3. 算例三

对于二维Schrödinger方程满足 $\Omega ={\left[-1,1\right]}^2\times \left(0,T\right]$，有

$\{\begin{array}{l}i\frac{\partial u}{\partial t}+\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}=f\\ u\left(x,y,0\right)=\sin \pi x\sin \pi y\\ u\left(-1,y,t\right)=0,u\left(1,y,t\right)=0\\ u\left(x,-1,t\right)=0,u\left(x,1,t\right)=0\end{array}$

其中

$\{\begin{array}{l}u=\sin \pi x\sin \pi y\cos t\\ f=i\sin t\sin \pi x\sin \pi y-{\pi }^2\cos t\sin \pi x\sin \pi y\end{array}$

下表中M + 1是为空间的节点数，N取1000等分。

Table 3. Maximum absolute error and maximum relative error of second-order finite difference scheme and barycentric interpolation collocation scheme ( τ = 10 − 3 )

表3. 二阶有限差分格式和重心插值配点格式的最大绝对误差和最大相对误差( $\tau ={10}^{-3}$)

固定时间节点，改变空间节点，表3给出了二维Schrödinger方程采用二阶有限差分格式和重心插值配点格式所得的最大绝对误差和最大相对误差。与一维Schrödinger方程类似，重心插值配点格式在空间上可以用较少的点便能达到更高的精度。比如，当 $\tau ={10}^{-3}$ 时，重心Lagrange插值只需要在空间上取8个点就可以达到10⁻⁴的精度，但二阶有限差分格式需要取40个点才能达到这个精度。

表4可知，对于二维Schrödinger方程，推导得出其在时间上是二阶精度的。

Table 4. The time convergence order of two-dimensional Schrödinger equation

表4. 二维Schrödinger方程的时间收敛阶

通过下面的图7和图8我们可以看出，重心插值配点格式的数值解和解析解的等势线图像基本一致，这进一步说明该数值格式的稳定性。

下图为M = 16，N = 1000时的图像：

画出两种格式的绝对误差图。由图9和图10可知，当时间剖分一致时，重心插值配点格式在空间上取16个点的时候，绝对误差精度就达到10⁻⁷，而二阶有限差分格式在空间上取31个点的时候，精度只有10⁻⁴。

Figure 7. Analytical solution of example three

图7. 算例三的解析解

Figure 8. Example three numerical solution of barycentric interpolation collocation scheme

图8. 算例三重心插值配点法格式数值解

Figure 9. Example three absolute error figure of second-order finite difference scheme

图9. 算例三二阶有限差分格式绝对误差图

Figure 10. Example three absolute error figure of barycentric interpolation collocation scheme

图10. 算例三重心插值配点格式绝对误差图

3.4. 算例四

对于二维Schrödinger方程满足 $\Omega =\left[\text{0},\text{2}\pi \right]\times \left[\text{0},\text{2}\pi \right]$，有

$\{\begin{array}{l}i{u}_t=-\frac{1}{2}\Delta u+V\left(X\right)+f,\left(X,t\right)\in \Omega \times \left(0,T\right]\\ u\left(X,0\right)=\sin x\sin y\end{array}$

其中

$V\left(X\right)=1-{\sin }^{2}x{\sin }^{2}y$

该问题的精确解为 $u={\text{e}}^{-2ti}\sin x\sin y$，再进一步针对Schrödinger方程的离散质量和能量守恒性进行验证。图11给出重心Lagrange插值格式对质量守恒量的保持情况，图12给出重心Lagrange插值格式对能量守恒量的保持情况。从图像能再一次看出我们的配点格式很好地满足离散质量和能量守恒性。

(a) 总质量 (b) 质量相对误差

Figure 11. Example four mass conservation

图11. 算例四质量守恒

(a) 总能量 (b) 能量相对误差

Figure 12. Example four energy conservation

图12. 算例四能量守恒

4. 结语

本文利用Crank-Nicolson格式结合重心Lagrange插值配点法求解Schrödinger方程一维二维问题，与经典的二阶有限差分格式进行数值比较，重心Lagrange插值配点格式在空间方向用较少的点时就可以达到很高的精度，且误差结果与理论分析相吻合，进一步数值算例也验证配点格式满足方程的质量和能量守恒性。今后，可以将重心插值配点格式推广到更高维Schrödinger方程和时间分数阶Schrödinger方程问题，也可将重心插值配点格式与其他格式相结合，探讨Schrödinger方程的高精度数值格式。

致谢

由衷感谢华侨大学数学科学学院翁智峰老师对本论文的指导与帮助。

基金项目

国家级大学生创新创业训练计划项目(S202110385029)；华侨大学实验教学与管理改革课题(SY2021J12)。

NOTES

^*通讯作者。