基于广义有限差分方法计算二维非线性耦合薛定谔方程

Calculating the Two-Dimensional Nonlinear Coupled Schrödinger Equation Based on the Generalized Finite Difference Method

摘要: 耦合非线性薛定谔方程组(Coupled Nonlinear Schrödinger Equations, CNLS)描述了许多量子系统的动力学行为，寻找一种有效的数值方法是科学和工程应用的关键。本文介绍了一种求解二维非线性耦合薛定谔方程的无网格方法，该方法称为广义有限差分法(GFDM)。该方法将Taylor多项式与最小二乘法相结合，在推导过程中，首先将CNLS非线性系统在时间方向使用线性化方法进行离散，推导出CNLS在时间上的线性化离散形式，再使用GFDM在空间上的计算。该方法为CNLS生成稀疏的非线性代数系统。最后给出了数值算例，表明了本文方法求解两个区域线性代数方程的有效性。数值结果表明，发展的数值方法是一种稳定、快速和精确的计算方法。

Abstract: Coupled nonlinear Schrödinger equations (CNLS) describe the dynamical behavior of many quantum systems, and finding an efficient numerical method is crucial for scientific and engineering applications. This paper introduces a meshless method for solving two-dimensional nonlinear coupled Schrödinger equations, called the generalized finite difference method (GFDM). This method combines Taylor polynomials with the least squares method. In the derivation, the CNLS nonlinear system is first discretized in the time direction using the Crank-Nicolson method, deriving the linearized discretized form of CNLS in time. Then, GFDM is used for spatial computation. This method generates a sparse nonlinear algebraic system for CNLS. Finally, numerical examples are given to demonstrate the effectiveness of the proposed method in solving linear algebraic equations in two regions. Numerical results show that the developed numerical method is a stable, fast, and accurate computational approach.

1. 引言

耦合非线性薛定谔方程组(Coupled Nonlinear Schrödinger Equations, CNLS)在量子物理、非线性光学、晶体物理、玻色–爱因斯坦凝聚和水波动力学等物理领域有着重要的应用[1]-[3]，近几十年来得到了广泛的研究。本章主要考虑下列耦合非线性薛定谔方程组[4]：

$\begin{array}{l}\text{i}\frac{\partial u\left(x,t\right)}{\partial t}+\Delta u\left(x,t\right)+\left(b_{11}{\left|u\left(x,t\right)\right|}^2+b_{12}{\left|v\left(x,t\right)\right|}^2\right)u\left(x,t\right)=f_1\left(x,t\right),x\in \Omega ,t\in \left(0,T\right],\\ \text{i}\frac{\partial v\left(x,t\right)}{\partial t}+\Delta v\left(x,t\right)+\left(b_{21}{\left|u\left(x,t\right)\right|}^2+b_{22}{\left|v\left(x,t\right)\right|}^2\right)v\left(x,t\right)=f_2\left(x,t\right),x\in \Omega ,t\in \left(0,T\right],\\ u\left(x,t\right)=\overline{u}\left(x,t\right),v\left(x,t\right)=\overline{v}\left(x,t\right),x\in \partial \Omega ,t\in \left(0,T\right],\\ u\left(x,0\right)=u^0\left(x\right),v\left(x,0\right)=v^0\left(x\right),x\in \Omega .\end{array}$(1.1)

其中$\text{i}=\sqrt{-1}$ ，$u$ 和$v$ 是复函数，表示两个耦合的场分量(如偏振模、能级态、原子种类等)，$x=\left[x,y\right]$ 是矢量坐标，$b_{11},b_{12},b_{21}$ 和$b_{22}$ 是非线性项强度的任意实数，表示自相互作用与交叉相互作用强度。在玻色–爱因斯坦凝聚体的系统中，该方程组模拟了两组分凝聚体[5]，其中$u\left(x,t\right)$ 和$v\left(x,t\right)$ 对应于组分的波函数。该方程组可以模拟双折射介质中的耦合模式[6]或多模光纤中的脉冲传播[7]，其中$u$ 和$v$ 表示沿着每个双折射轴或光纤模式的电场强度。$f_1,f_2$ 是源项，$\Omega$ 是边界为$\partial \Omega$ 的有界域，$\overline{u},\overline{v},u^0,v^0$ 是已知连续函数。

近年来，一种称为广义有限差分法(GFDM)的无网格技术在工程和科学界得到了广泛的应用[8] [9]，因为它可以很容易地扩展到不规则形状区域上的高维模型。本文的目的是发展一种用于求解不规则形状区域上的二维CNLS的GFDM。GFDM将Taylor多项式与最小二乘法相结合，该方法为CNLS生成稀疏的非线性代数系统。

2. 时间离散化方法

现在，我们在方程(1.1)中构造一个关于时间方向的耦合CNLS的半离散格式。取N为一个正整数，$\tau =T/N$ 为时间步长，使用网格点${\left\{t_n=n\tau \right\}}_{n=0}^N$ 将时域$\left[0,T\right]$ 划分为N个子区间。为简单起见，记

$t_{n+\frac{1}{2}}=\frac{t_n+t_{n+1}}{2},u\left(x,t_n\right)=u^n\left(x\right),v\left(x,t_n\right)=v^n\left(x\right),f\left(x,t_n\right)=f^n\left(x\right),0\le n\le N.$

考虑方程(1.1)在点$\left(x,t_{n+\frac{1}{2}}\right)$ 处，我们得到

$\begin{array}{c}\text{i}\frac{\partial u\left(x,t_{n+\frac{1}{2}}\right)}{\partial t}+\Delta u\left(x,t_{n+\frac{1}{2}}\right)+\left(b_{11}{\left|u\left(x,t_{n+\frac{1}{2}}\right)\right|}^2+b_{12}{\left|v\left(x,t_{n+\frac{1}{2}}\right)\right|}^2\right)u\left(x,t_{n+\frac{1}{2}}\right)=f_1\left(x,t_{n+\frac{1}{2}}\right),\\ \text{i}\frac{\partial v\left(x,t_{n+\frac{1}{2}}\right)}{\partial t}+\Delta v\left(x,t_{n+\frac{1}{2}}\right)+\left(b_{21}{\left|u\left(x,t_{n+\frac{1}{2}}\right)\right|}^2+b_{22}{\left|v\left(x,t_{n+\frac{1}{2}}\right)\right|}^2\right)v\left(x,t_{n+\frac{1}{2}}\right)=f_2\left(x,t_{n+\frac{1}{2}}\right),\end{array}$

这里$0\le n\le N-1$ ，然后我们有以下离散化：

$\begin{array}{l}\text{i}\frac{u\left(x,t_{n+1}\right)-u\left(x,t_n\right)}{\tau }+\frac{1}{2}\left[\Delta u\left(x,t_{n+1}\right)+\Delta u\left(x,t_n\right)\right]+\left(b_{11}{\left|u\left(x,t_{n+\frac{1}{2}}\right)\right|}^2+b_{12}{\left|v\left(x,t_{n+\frac{1}{2}}\right)\right|}^2\right)\\ \cdot \frac{1}{2}\left[u\left(x,t_{n+1}\right)+u\left(x,t_n\right)\right]=\frac{1}{2}\left[f_1\left(x,t_{n+1}\right)+f_1\left(x,t_n\right)\right]+O\left({\tau }^2\right),\\ \text{i}\frac{v\left(x,t_{n+1}\right)-v\left(x,t_n\right)}{\tau }+\frac{1}{2}\left[\Delta v\left(x,t_{n+1}\right)+\Delta v\left(x,t_n\right)\right]+\left(b_{21}{\left|u\left(x,t_{n+\frac{1}{2}}\right)\right|}^2+b_{22}{\left|v\left(x,t_{n+\frac{1}{2}}\right)\right|}^2\right)\\ \cdot \frac{1}{2}\left[v\left(x,t_{n+1}\right)+v\left(x,t_n\right)\right]=\frac{1}{2}\left[f_2\left(x,t_{n+1}\right)+f_2\left(x,t_n\right)\right]+O\left({\tau }^2\right).\end{array}$ (2.1)

通过Taylor展开式很容易得到$u\left(x,t\right),v\left(x,t\right)$ 在$t=t_{n+\frac{1}{2}}$ 处的三阶近似，对于$n\ge 2$ ，

$\begin{array}{l}u\left(x,t_{n+\frac{1}{2}}\right)=\frac{3}{8}u\left(x,t_{n-2}\right)-\frac{5}{4}u\left(x,t_{n-1}\right)+\frac{15}{8}u\left(x,t_n\right)+O\left({\tau }^3\right),\\ v\left(x,t_{n+\frac{1}{2}}\right)=\frac{3}{8}v\left(x,t_{n-2}\right)-\frac{5}{4}v\left(x,t_{n-1}\right)+\frac{15}{8}v\left(x,t_n\right)+O\left({\tau }^3\right),\end{array}$

在$n=1$ 时，

$\begin{array}{l}u\left(x,t_{\frac{3}{2}}\right)=\frac{3}{2}u\left(x,t_1\right)-\frac{1}{2}u\left(x,t_0\right)+O\left({\tau }^2\right),\\ v\left(x,t_{\frac{3}{2}}\right)=\frac{3}{2}v\left(x,t_1\right)-\frac{1}{2}v\left(x,t_0\right)+O\left({\tau }^2\right).\end{array}$

为了保证数值格式的精确性，我们考虑$n=0$ 时非线性项的近似。在等式(1.1)中，令$t\to 0$ ，可以得到

$\begin{array}{l}\frac{\partial u\left(x,0\right)}{\partial t}=\text{i}\Delta u^0\left(x\right)+\text{i}\left(b_{11}{\left|u^0\left(x\right)\right|}^2+b_{12}{\left|v^0\left(x\right)\right|}^2\right)u^0\left(x\right)-\text{i}{f}_1^0\left(x\right),x\in \Omega ,\\ \frac{\partial v\left(x,0\right)}{\partial t}=\text{i}\Delta v^0\left(x\right)+\text{i}\left(b_{21}{\left|u^0\left(x\right)\right|}^2+b_{22}{\left|v^0\left(x\right)\right|}^2\right)v^0\left(x\right)-\text{i}{f}_2^0\left(x\right),x\in \Omega .\end{array}$

为了简单起见，记

${\widehat{u}}^n\left(x\right)=\{\begin{array}{l}{u}^0\left(x\right)+\frac{1}{2}[\text{i}\Delta u^0\left(x\right)+\text{i}\left(b_{11}{\left|u^0\left(x\right)\right|}^2+b_{12}{\left|v^0\left(x\right)\right|}^2\right)u^0\left(x\right)-\text{i}{f}_1^0\left(x\right)],n=0,\\ \frac{3}{2}u\left(x,t_1\right)-\frac{1}{2}u\left(x,t_0\right),n=1,\\ \frac{3}{8}u\left(x,t_{n-2}\right)-\frac{5}{4}u\left(x,t_{n-1}\right)+\frac{15}{8}u\left(x,t_n\right),n\ge 2.\end{array}$

${\widehat{v}}^n\left(x\right)=\{\begin{array}{l}{v}^0\left(x\right)+\frac{1}{2}[\text{i}\Delta v^0\left(x\right)+\text{i}\left(b_{21}{\left|u^0\left(x\right)\right|}^2+b_{22}{\left|v^0\left(x\right)\right|}^2\right)v^0\left(x\right)-\text{i}{f}_2^0\left(x\right)],n=0,\\ \frac{3}{2}v\left(x,t_1\right)-\frac{1}{2}v\left(x,t_0\right),n=1,\\ \frac{3}{8}v\left(x,t_{n-2}\right)-\frac{5}{4}v\left(x,t_{n-1}\right)+\frac{15}{8}v\left(x,t_n\right),n\ge 2.\end{array}$

将${\widehat{u}}^n\left(x\right),{\widehat{v}}^n\left(x\right)$ 代入方程(2.1)，并分别用$u^n\left(x\right),v^n\left(x\right)$ 逼近方程(2.1)中的$u\left(x,t_n\right),v\left(x,t_n\right)$ ，得到

$\begin{array}{l}2\text{i}\frac{u^{n+1}\left(x\right)-u^n\left(x\right)}{\tau }+\Delta u^{n+1}\left(x\right)+\Delta u^n\left(x\right)+\left(b_{11}{\left|{\widehat{u}}^n\left(x\right)\right|}^2+b_{12}{\left|{\widehat{v}}^n\left(x\right)\right|}^2\right)\left[u^{n+1}\left(x\right)+u^n\left(x\right)\right]\\ =f_1^{n+1}\left(x\right)+f_1^n\left(x\right),0\le n\le N-1\\ 2\text{i}\frac{v^{n+1}\left(x\right)-v^n\left(x\right)}{\tau }+\Delta v^{n+1}\left(x\right)+\Delta v^n\left(x\right)+\left(b_{21}{\left|{\widehat{u}}^n\left(x\right)\right|}^2+b_{22}{\left|{\widehat{v}}^n\left(x\right)\right|}^2\right)\left[v^{n+1}\left(x\right)+v^n\left(x\right)\right]\\ =f_2^{n+1}\left(x\right)+f_2^n\left(x\right),0\le n\le N-1\end{array}$(2.2)

重新排列方程(2.2)中的各项，我们得到一个线性化的时间离散格式，如下所示：

$\begin{array}{l}{k}_1^n{u}^{n+1}\left(x\right)+\Delta u^{n+1}\left(x\right)={\Phi }_1^n\left(x\right),0\le n\le N-1,\\ k_2^n{v}^{n+1}\left(x\right)+\Delta v^{n+1}\left(x\right)={\Phi }_2^n\left(x\right),0\le n\le N-1,\end{array}$ (2.3)

其中

$\begin{array}{l}{k}_1^n=\frac{2\text{i}}{\tau }+\left(b_{11}{\left|{\widehat{u}}^n\left(x\right)\right|}^2+b_{12}{\left|{\widehat{v}}^n\left(x\right)\right|}^2\right),k_2^n=\frac{2\text{i}}{\tau }+\left(b_{21}{\left|{\widehat{u}}^n\left(x\right)\right|}^2+b_{22}{\left|{\widehat{v}}^n\left(x\right)\right|}^2\right),\\ {\Phi }_1^n\left(x\right)=\left[\frac{2\text{i}}{\tau }-\left(b_{11}{\left|{\widehat{u}}^n\left(x\right)\right|}^2+b_{12}{\left|{\widehat{v}}^n\left(x\right)\right|}^2\right)\right]u^n\left(x\right)-\Delta u^n\left(x\right)+f_1^{n+1}\left(x\right)+f_1^n\left(x\right),\\ {\Phi }_2^n\left(x\right)=\left[\frac{2\text{i}}{\tau }-\left(b_{21}{\left|{\widehat{u}}^n\left(x\right)\right|}^2+b_{22}{\left|{\widehat{v}}^n\left(x\right)\right|}^2\right)\right]v^n\left(x\right)-\Delta v^n\left(x\right)+f_2^{n+1}\left(x\right)+f_2^n\left(x\right).\end{array}$

3. 空间近似的GFDM

为了简单起见，在半离散格式(2.3)中忽略时间上标，令$u\left(x\right),v\left(x\right)$ 表示微分方程当前时刻的未知量，${\Phi }_1\left(x\right),{\Phi }_2\left(x\right)$ 表示已知量，包含时间层中已知的函数数值解以及源项，$k_1,k_2$ 表示非线性项的线性近似。在时间离散化之后，在每个时间层$t_n\left(n=1,\cdots ,N\right)$ 处得到二阶线性微分方程

$\begin{array}{l}{k}_{1}u\left(x\right)+\Delta u\left(x\right)={\Phi }_1\left(x\right),\\ k_{2}v\left(x\right)+\Delta v\left(x\right)={\Phi }_2\left(x\right),\end{array}$ (3.1)

以及边界条件

$u\left(x\right)=\overline{u}\left(x\right),v\left(x\right)=\overline{v}\left(x\right),x\in \partial \Omega .$

GFDM是求解数学物理方程的一种有效方法[10] [11]。在本节中，我们应用GFDM方法求解方程(3.1)中的二阶线性系统的解。接下来我们介绍GFDM在二维情形($\Omega \subset R^2$ )下的实现。取计算区域内一点$x_0=\left(x_0,y_0\right)$ 作为中心节点，距离中心节点$x_0$ 最近的$m$ 个点作为其支持节点(如图1所示)，设$u_0,v_0$ 为中心节点$x_0$ 处的函数值，$u_i,v_i$ 为$m$ 个支持节点$x_i=\left(x_i,y_i\right),i=1,\cdots ,m$ 处的函数值。

考虑到上面的微分方程组，在每个局部子域中，以$x_0$ 为中心，对支持节点$x_i$ 进行Taylor级数展开：

Figure 1. Diagram of central node and supporting nodes

图1. 中心节点与支持节点示意图

$\begin{array}{l}u\left(x_i\right)={\displaystyle \sum _{j=0}^{\infty }\frac{1}{j!}}{\left(h_i\frac{\partial }{\partial x}+k_i\frac{\partial }{\partial y}\right)}^{j}u\left(x_0\right),i=1,2,\cdots ,m,\\ v\left(x_i\right)={\displaystyle \sum _{j=0}^{\infty }\frac{1}{j!}}{\left(h_i\frac{\partial }{\partial x}+k_i\frac{\partial }{\partial y}\right)}^{j}v\left(x_0\right),i=1,2,\cdots ,m,\end{array}$

其中$h_i=x_i-x_0,k_i=y_i-y_0$ ，通过截断四阶导数之后的部分，可以定义如下的残差函数：

$\begin{array}{c}R\left(u\right)={\displaystyle \sum _{i=1}^m[(u_0-u_i+h_i\frac{\partial u_0}{\partial x}+k_i\frac{\partial u_0}{\partial y}+\frac{h_i^2}{2}\frac{{\partial }^2{u}_0}{\partial x^2}+\frac{k_i^2}{2}\frac{{\partial }^2{u}_0}{\partial y^2}+h_i{k}_i\frac{{\partial }^2{u}_0}{\partial x\partial y}}\\ +\frac{h_i^3}{6}\frac{{\partial }^3{u}_0}{\partial x^3}+\frac{k_i^3}{6}\frac{{\partial }^3{u}_0}{\partial y^3}+\frac{h_i^2{k}_i}{2}\frac{{\partial }^3{u}_0}{\partial x^2\partial y}+\frac{h_i{k}_i^2}{2}\frac{{\partial }^3{u}_0}{\partial x\partial y^2}+\frac{h_i^4}{24}\frac{{\partial }^4{u}_0}{\partial x^4}\\ +\frac{k_i^4}{24}\frac{{\partial }^4{u}_0}{\partial y^4}+\frac{h_i^3{k}_i}{6}\frac{{\partial }^4{u}_0}{\partial x^3\partial y}+\frac{h_i{k}_i^3}{6}\frac{{\partial }^4{u}_0}{\partial x\partial y^3}+\frac{h_i^2{k}_i^2}{4}\frac{{\partial }^4{u}_0}{\partial x^2\partial y^2}){\omega }_i^2],\end{array}$

$\begin{array}{c}R\left(v\right)={\displaystyle \sum _{i=1}^m[(v_0-v_i+h_i\frac{\partial v_0}{\partial x}+k_i\frac{\partial v_0}{\partial y}+\frac{h_i^2}{2}\frac{{\partial }^2{v}_0}{\partial x^2}+\frac{k_i^2}{2}\frac{{\partial }^2{v}_0}{\partial y^2}+h_i{k}_i\frac{{\partial }^2{v}_0}{\partial x\partial y}}\\ +\frac{h_i^3}{6}\frac{{\partial }^3{v}_0}{\partial x^3}+\frac{k_i^3}{6}\frac{{\partial }^3{v}_0}{\partial y^3}+\frac{h_i^2{k}_i}{2}\frac{{\partial }^3{v}_0}{\partial x^2\partial y}+\frac{h_i{k}_i^2}{2}\frac{{\partial }^3{v}_0}{\partial x\partial y^2}+\frac{h_i^4}{24}\frac{{\partial }^4{v}_0}{\partial x^4}\\ +\frac{k_i^4}{24}\frac{{\partial }^4{v}_0}{\partial y^4}+\frac{h_i^3{k}_i}{6}\frac{{\partial }^4{v}_0}{\partial x^3\partial y}+\frac{h_i{k}_i^3}{6}\frac{{\partial }^4{v}_0}{\partial x\partial y^3}+\frac{h_i^2{k}_i^2}{4}\frac{{\partial }^4{v}_0}{\partial x^2\partial y^2}){\omega }_i^2],\end{array}$

其中${\omega }_i$ 为加权函数[12]，在本文中为

${\omega }_{\text{i}}=1-6{\left(\frac{d_i}{d_m}\right)}^2+8{\left(\frac{d_i}{d_m}\right)}^3-3{\left(\frac{d_i}{d_m}\right)}^4,d_i\le d_m,i=1,\cdots ,m,$

这里$d_i=\Vert x-x_0\Vert$ 为中心节点到支持节点的范数，$d_m=\max \left\{d_i,i=1,\cdots ,m\right\}$ 。加权函数来表示不同节点处所取近似值的比重，区域内的支持节点离中心节点越近，这个节点的计算权重就越大，这样所得到的近似值就越精确。

为了使残差最小化，此处对$R\left(u\right),R\left(v\right)$ 中的未知变量

$\begin{array}{l}{D}_u={\left[\frac{\partial u_0}{\partial x},\frac{\partial u_0}{\partial y},\frac{{\partial }^2{u}_0}{\partial x^2},\frac{{\partial }^2{u}_0}{\partial y^2},\frac{{\partial }^2{u}_0}{\partial x\partial y},\frac{{\partial }^3{u}_0}{\partial x^3},\frac{{\partial }^3{u}_0}{\partial y^3},\frac{{\partial }^3{u}_0}{\partial x^2\partial y},\frac{{\partial }^3{u}_0}{\partial x\partial y^2},\frac{{\partial }^4{u}_0}{\partial x^4},\frac{{\partial }^4{u}_0}{\partial y^4},\frac{{\partial }^4{u}_0}{\partial x^3\partial y},\frac{{\partial }^4{u}_0}{\partial x\partial y^3},\frac{{\partial }^4{u}_0}{\partial x^2\partial y^2}\right]}_{14\times 1}^T,\\ D_v={\left[\frac{\partial v_0}{\partial x},\frac{\partial v_0}{\partial y},\frac{{\partial }^2{v}_0}{\partial x^2},\frac{{\partial }^2{v}_0}{\partial y^2},\frac{{\partial }^2{v}_0}{\partial x\partial y},\frac{{\partial }^3{v}_0}{\partial x^3},\frac{{\partial }^3{v}_0}{\partial y^3},\frac{{\partial }^3{v}_0}{\partial x^2\partial y},\frac{{\partial }^3{v}_0}{\partial x\partial y^2},\frac{{\partial }^4{v}_0}{\partial x^4},\frac{{\partial }^4{v}_0}{\partial y^4},\frac{{\partial }^4{v}_0}{\partial x^3\partial y},\frac{{\partial }^4{v}_0}{\partial x\partial y^3},\frac{{\partial }^4{v}_0}{\partial x^2\partial y^2}\right]}_{14\times 1}^T,\end{array}$

求极小值，即令

$\frac{\partial R\left(u\right)}{\partial D_u}=0,\frac{\partial R\left(v\right)}{\partial D_v}=0.$

这样便可得到关于未知变量$D_u,D_v$ 的两个线性方程组

$A{D}_u=b_u,A{D}_v=b_v,$

此处线性方程组的系数矩阵$A$ ，右端项向量$b_u,b_v$ 分别为

$A=P^{T}WP,$

$b_u=BU,b_v=BV,$

其中

$\begin{array}{l}W=diag\left\{{\omega }_1^2,{\omega }_2^2,\cdots ,{\omega }_m^2\right\},\\ p_i={\left[h_i,k_i,\frac{h_i^2}{2},\frac{k_i^2}{2},h_i{k}_i,\frac{h_i^3}{6},\frac{k_i^3}{6},\frac{h_i^2{k}_i}{2},\frac{h_i{k}_i^2}{2},\frac{h_i^4}{24},\frac{k_i^4}{24},\frac{h_i^3{k}_i}{6},\frac{h_i{k}_i^3}{6},\frac{h_i^2{k}_i^2}{4}\right]}_{1\times 14},i=1,2,\cdots ,m,\\ P=\left[\begin{array}{c}{p}_1\\ p_2\\ ⋮\\ p_m\end{array}\right]=\left[\begin{array}{cccc}{h}_1& k_1& \cdots & \frac{h_1^2{k}_1^2}{4}\\ h_2& k_2& \cdots & \frac{h_2^2{k}_2^2}{4}\\ ⋮& ⋮& \ddots & ⋮\\ h_m& k_m& \cdots & \frac{h_m^2{k}_m^2}{4}\end{array}\right],\\ B={\left[-{\displaystyle \sum _{i=1}^m{\omega }_i^2{p}_i^{\text{T}},{\omega }_1^2{p}_1^{\text{T}},{\omega }_2^2{p}_2^{\text{T}},\cdots ,{\omega }_m^2{p}_m^{\text{T}}}\right]}_{1\times m+1},\\ U={\left[u_0,u_1,u_2,\cdots ,u_m\right]}_{1\times m+1},\\ V={\left[v_0,v_1,v_2,\cdots ,v_m\right]}_{1\times m+1}.\end{array}$

由此未知量的各阶偏导数$D_u,D_v$ 可被表示成如下形式：

$D_u=A^{-1}{b}_u=A^{-1}BU=EU,D_v=A^{-1}{b}_v=A^{-1}BV=EV.$

上述运算过程可以将$x_0$ 处函数值$u_0,v_0$ 的各阶偏导数表示为其支持点$u_i,v_i$ 函数值的线性组合，对于二维耦合非线性薛定谔方程时间离散之后的格式(3.1)，我们可以得到

$\frac{{\partial }^2{u}_0}{\partial x^2}=e_{3,1}^0{u}_0+{\displaystyle \sum _{i=1}^m{e}_{3,i+1}^0{u}_i},\frac{{\partial }^2{u}_0}{\partial y^2}=e_{4,1}^0{u}_0+{\displaystyle \sum _{i=1}^m{e}_{4,i+1}^0{u}_i},$

以及

$\frac{{\partial }^2{v}_0}{\partial x^2}=e_{3,1}^0{v}_0+{\displaystyle \sum _{i=1}^m{e}_{3,i+1}^0{v}_i},\frac{{\partial }^2{v}_0}{\partial y^2}=e_{4,1}^0{v}_0+{\displaystyle \sum _{i=1}^m{e}_{4,i+1}^0{v}_i},$

其中$e_{p,q}^0$ 是矩阵$E$ 中的元素。在整个计算区域内重复如上的过程，最终实现全域离散。值得注意的是，支持节点的选取不依赖网格，这意味着整个离散过程可以在不规则区域上进行。

我们对方程组(3.1)分开求解$u\left(x\right),v\left(x\right)$ ，在整个计算区域内离散布点，令$M_{inp},M_{bp}$ 分别表示计算区域内部以及边界∂Ω上的离散点的数量。首先对于方程

$\begin{array}{l}{k}_{1}u\left(x\right)+\Delta u\left(x\right)={\Phi }_1\left(x\right),\\ u\left(x\right)=\overline{u}\left(x\right),x\in \partial \Omega ,\end{array}$

令区域内每个节点满足此方程，其中边界点满足$M_{bp}$ 条边界条件，内点满足$M_{inp}$ 条控制方程。最终得到了有$M_{inp}+M_{bp}$ 条方程的代数方程组，其矩阵表示为

$\left[\begin{array}{cc}{I}_{M_{bp}\times M_{bp}}& 0\\ 0& k_1{I}_{M_{inp}\times M_{inp}}\end{array}\right]u=\left[\begin{array}{c}u\left(x_1\right)\\ ⋮\\ u\left(x_{M_{bp}}\right)\\ {\Phi }_1\left(x_{M_{bp}+1}\right)\\ ⋮\\ {\Phi }_1\left(x_{M_{bp}+M_{inp}}\right)\end{array}\right],$ (3.2)

同样的，对于方程

$\begin{array}{l}{k}_{2}v\left(x\right)+\Delta v\left(x\right)={\Phi }_2\left(x\right),\\ v\left(x\right)=\overline{v}\left(x\right),x\in \partial \Omega ,\end{array}$

有$M_{inp}+M_{bp}$ 条方程的代数方程组，矩阵表示为

$\left[\begin{array}{cc}{I}_{M_{bp}\times M_{bp}}& 0\\ 0& k_2{I}_{M_{inp}\times M_{inp}}\end{array}\right]v=\left[\begin{array}{c}v\left(x_1\right)\\ ⋮\\ v\left(x_{M_{bp}}\right)\\ {\Phi }_2\left(x_{M_{bp}+1}\right)\\ ⋮\\ {\Phi }_2\left(x_{M_{bp}+M_{inp}}\right)\end{array}\right],$ (3.3)

这样最终得到了两个大型且易于计算的稀疏系数矩阵，如果改变计算区域的几何形状，仅更改区域内点的分布，对于不同的边界条件，也仅需使其满足对应的约束条件，这使得所提方法能够处理许多复杂的计算区域以及应对不同的边界条件问题。另外，需要注意的是，方程(3.2)和(3.3)中的系数矩阵是稀疏的，这样可以提高计算速度。

4. 数值算例

在本节中，我们提供了三个数值例子来测试所提出的GFDM求解任意区域的CNLS的效率。在所有的数值例子中，我们在使用GFDM进行了4阶泰勒展开。为了测试GFDM求解CNLS的准确性和稳定性。首先，定义误差范数：

$L^{\infty }=\max \left|{\tilde{u}}_k-u_k\right|,k=1,2,\cdots ,M,$

$RE=\sqrt{\frac{{\displaystyle \sum _{k=1}^M{\left|{\tilde{u}}_k-u_k\right|}^2}}{{\displaystyle \sum _{k=1}^M{\left|u_k\right|}^2}}},$

以及空间收敛阶$p_1$ 和时间收敛阶$p_2$ ：

$\begin{array}{l}{p}_1=\frac{\log \left[E\left(h_1,\tau \right)/E\left(h_2,\tau \right)\right]}{\log \left(h_1/h_2\right)},\\ p_2=\frac{\log \left[E\left(h,{\tau }_1\right)/E\left(h,{\tau }_2\right)\right]}{\log \left({\tau }_1/{\tau }_2\right)},\end{array}$

其中${\tilde{u}}_k,u_k$ 分别表示离散点$x_i$ 处的数值解和精确解。$M$ 为计算区域内离散点数，$E$ 表示误差范数。

4.1. 收敛性分析

本算例考虑方形计算区域，$\Omega =0\le x,y\le 1$ ，边界为$\partial \Omega$ ，将计算域以步长$h$ 随机均匀布点(如图2所示)，由方程(1.1)中$b_{11}=b_{21}=2,b_{12}=b_{22}=1$ 所定义，即控制方程为：

$\begin{array}{l}\text{i}\frac{\partial u\left(x,t\right)}{\partial t}+\Delta u\left(x,t\right)+\left(2{\left|u\left(x,t\right)\right|}^2+{\left|v\left(x,t\right)\right|}^2\right)u\left(x,t\right)=f_1\left(x,t\right),x\in \Omega ,t\in \left(0,T\right],\\ \text{i}\frac{\partial v\left(x,t\right)}{\partial t}+\Delta v\left(x,t\right)+\left(2{\left|u\left(x,t\right)\right|}^2+{\left|v\left(x,t\right)\right|}^2\right)v\left(x,t\right)=f_2\left(x,t\right),x\in \Omega ,t\in \left(0,T\right],\end{array}$

其边界条件为Dirichlet 条件。该问题的精确解为：

$\begin{array}{l}u\left(x,t\right)=\mathrm{sech}\left(x+y-4t\right){\text{e}}^{\text{i}\left(2x+2y-3t\right)},\\ v\left(x,t\right)=\mathrm{sech}\left(x+y+2t\right){\text{e}}^{\text{i}\left(x+y-t\right)},\end{array}$

考虑初值

$\begin{array}{l}{u}^0\left(x\right)=\mathrm{sech}\left(x+y\right){\text{e}}^{\text{i}\left(2x+2y\right)},\\ v^0\left(x\right)=\mathrm{sech}\left(x+y\right){\text{e}}^{\text{i}\left(x+y\right)}.\end{array}$

Figure 2. Schematic diagram of a square computational domain with randomly uniformly distributed points

图2. 方形计算域随机均匀布点示意图

首先，我们检验数值格式的空间精度，表1与表2给出了该数值格式在$\tau =1/200,T=1$ 时的空间误差和收敛阶。然后通过空间布点距离期望$h=1/100$ 来检验该数值格式的时间精度，表3与表4给出了数值结果。表1与表2中的数据表明数值格式的空间收敛阶为$O\left(h^4\right)$ ，与离散格式的推导分析是一致的。另外，表3与表4中的数据表明数值格式的时间收敛阶为$O\left({\tau }^2\right)$ ，这意味着我们对非线性项的线性近似求解对时间方向上的推导分析是一致的。

Table 1. $L^{\infty },RE$ and convergence order of function $u$ ($\tau =1/200,T=1$ )

表1. 函数$u$ 的$L^{\infty },RE$ 和收敛阶($\tau =1/200,T=1$ )

$h$	$L^{\infty }$	收敛阶$p_1$	$RE$	收敛阶$p_2$
1/60	$8.41\times {10}^{-5}$	——	$3.46\times {10}^{-4}$	——
1/80	$2.64\times {10}^{-5}$	4.02	$1.14\times {10}^{-4}$	3.85
1/100	$1.09\times {10}^{-5}$	3.97	$4.86\times {10}^{-5}$	3.83
1/120	$5.51\times {10}^{-6}$	3.74	$2.26\times {10}^{-5}$	4.19
1/140	$3.14\times {10}^{-6}$	3.64	$1.28\times {10}^{-5}$	3.71

Table 2. $L^{\infty },RE$ and convergence order of function $v$ ($\tau =1/200,T=1$ )

表2. 函数$v$ 的$L^{\infty },RE$ 和收敛阶($\tau =1/200,T=1$ )

$h$	$L^{\infty }$	收敛阶$p_1$	$RE$	收敛阶$p_2$
1/60	$9.56\times {10}^{-6}$	——	$4.04\times {10}^{-5}$	——
1/80	$2.92\times {10}^{-6}$	4.13	$1.29\times {10}^{-5}$	3.97
1/100	$1.18\times {10}^{-6}$	4.07	$5.28\times {10}^{-6}$	3.99
1/120	$5.64\times {10}^{-7}$	4.03	$2.63\times {10}^{-6}$	3.82
1/140	$3.06\times {10}^{-7}$	3.97	$1.46\times {10}^{-6}$	3.80

Table 3. $L^{\infty },RE$ and convergence order of function $u$ ($h=1/100,T=1$ )

表3. 函数$u$ 的$L^{\infty },RE$ 和收敛阶($h=1/100,T=1$ )

$\tau$	$L^{\infty }$	收敛阶$p_1$	$RE$	收敛阶$p_2$
1/100	$4.38\times {10}^{-5}$	——	$1.94\times {10}^{-4}$	——
1/120	$3.10\times {10}^{-5}$	1.89	$1.38\times {10}^{-4}$	1.89
1/140	$2.23\times {10}^{-5}$	2.13	$1.02\times {10}^{-4}$	1.97
1/160	$1.71\times {10}^{-5}$	1.98	$7.59\times {10}^{-5}$	2.18
1/180	$1.35\times {10}^{-5}$	2.01	$5.96\times {10}^{-5}$	2.06

Table 4. $L^{\infty },RE$ and convergence order of function $v$ ($h=1/100,T=1$ )

表4. 函数$v$ 的$L^{\infty },RE$ 和收敛阶($h=1/100,T=1$ )

$\tau$	$L^{\infty }$	收敛阶$p_1$	$RE$	收敛阶$p_2$
1/100	$4.51\times {10}^{-6}$	——	$2.04\times {10}^{-5}$	——
1/120	$3.20\times {10}^{-6}$	1.88	$1.43\times {10}^{-5}$	1.96
1/140	$2.39\times {10}^{-6}$	1.89	$1.07\times {10}^{-5}$	1.89
1/160	$1.82\times {10}^{-6}$	2.03	$8.05\times {10}^{-6}$	2.10
1/180	$1.45\times {10}^{-6}$	1.95	$6.40\times {10}^{-6}$	1.95

4.2. 时间步误差计算

本算例考虑圆形计算区域(以点$\left(1,1\right)$ 为圆心，半径为1的封闭圆和计算区域被定义在$0\le x,y\le 1$ 内的五角星状不规则区域(如图3所示)，由方程(1.1)中$b_{11}=b_{21}=1,$ $b_{12}=b_{22}=1$ 所定义，即控制方程为：

$\begin{array}{l}\text{i}\frac{\partial u\left(x,t\right)}{\partial t}+\Delta u\left(x,t\right)+\left({\left|u\left(x,t\right)\right|}^2+{\left|v\left(x,t\right)\right|}^2\right)u\left(x,t\right)=f_1\left(x,t\right),x\in \Omega ,t\in \left(0,T\right],\\ \text{i}\frac{\partial v\left(x,t\right)}{\partial t}+\Delta v\left(x,t\right)+\left({\left|u\left(x,t\right)\right|}^2+{\left|v\left(x,t\right)\right|}^2\right)v\left(x,t\right)=f_2\left(x,t\right),x\in \Omega ,t\in \left(0,T\right],\end{array}$

其边界条件为Dirichlet 条件。该问题的精确解为：

$\begin{array}{l}u\left(x,t\right)=\sin \left(4x+y+2\right){\text{e}}^{\text{i}\left(2x+2y-3t+2\right)},\\ v\left(x,t\right)=\cos \left(2x+2y+1\right){\text{e}}^{\text{i}\left(x+y-t\right)},\end{array}$

考虑初值

$\begin{array}{l}{u}^0\left(x\right)=\sin \left(4x+y+2\right){\text{e}}^{\text{i}\left(2x+2y+2\right)},\\ v^0\left(x\right)=\cos \left(2x+2y+1\right){\text{e}}^{\text{i}\left(x+y\right)}.\end{array}$

(a) (b)

Figure 3. Schematic diagram of randomly uniformly distributed points in circular and five-pointed star-shaped computational domains

图3. 圆形计算域和五角星状计算域随机均匀布点示意图

我们在圆形计算域上分布了4384个配置节点，其中包括126个边界节点和4258个内部节点。我们将时间区间从$t=0$ 到$t=1$ ，并且时间步长被设置为$\tau =0.01$ 。数值实验得到函数$u,v$ 的每个时间步$L^{\infty },RE$ 误差如图4所示。从图中结果看出，在圆形计算区域上，所使用方法能保持较高的精度。在五角星状不规则区域上随机分布了上分布了9560个配置节点，其中包括258个边界节点和9302个内部节点。我们将时间区间从$t=0$ 到$t=1$ ，并且时间步长被设置为$\tau =0.01$ 。数值实验得到函数$u,v$ 的每个时间步$L^{\infty },RE$ 误差如图5所示。从图中结果看出，在五角星状不规则计算区域上，所使用方法同样能保持较高的精度。以上可以表明线性化格式的GFDM在解决问题上提供了有效的解决方案。

(a) (b)

Figure 4. $L^{\infty }$ and $RE$ of function $u$ , $v$ (circular computational domain)

图4. 函数$u$ ，$v$ 的$L^{\infty },RE$ (圆形计算域)

(a) (b)

Figure 5. $L^{\infty }$ and $RE$ of function $u$ , $v$ (five-pointed star-shaped computational domain)

图5. 函数$u$ ，$v$ 的$L^{\infty },RE$ (五角星状计算域)

4.3. 亮孤子传播数值实验

亮孤子的传播行为是非线性科学中研究的重要方向之一[12]，尤其在二维空间中，其演化特性更为复杂和丰富。亮孤子作为一种局域化波动结构，能够在具有聚焦型非线性特性的介质中传播时保持稳定的形态，不发生扩散或塌缩。与线性波不同，亮孤子的形成与维持依赖于非线性效应与色散效应之间的精确平衡。在二维空间中，由于自由度的增加，这一平衡更加脆弱，孤子的稳定性与动力学演化更具有挑战性，因此成为实验与理论研究的热点。

本节的目标是利用我们的方法模拟两个空间上偏移的二维矢量亮孤子在二维空间中的传播过程，全面展示其从初始形成到稳定传播的演化全过程。实验中，我们引入具有正非线性系数的二维非线性薛定谔方程作为基础模型，通过调控初始波包的幅度与宽度，构造出符合稳定传播条件的亮孤子。随后，通过数值模拟观察亮孤子在自由传播过程中是否能够保持结构稳定，避免扩散或塌缩。为了构造两个孤子的传播模拟，由方程(1.1)中$b_{11}=b_{21}=1,b_{12}=b_{22}=1$ 所定义，即控制方程为：

$\begin{array}{l}\text{i}\frac{\partial u\left(x,t\right)}{\partial t}+\Delta u\left(x,t\right)+\left({\left|u\left(x,t\right)\right|}^2+{\left|v\left(x,t\right)\right|}^2\right)u\left(x,t\right)=f_1\left(x,t\right),x\in \Omega ,t\in \left(0,T\right],\\ \text{i}\frac{\partial v\left(x,t\right)}{\partial t}+\Delta v\left(x,t\right)+\left({\left|u\left(x,t\right)\right|}^2+{\left|v\left(x,t\right)\right|}^2\right)v\left(x,t\right)=f_2\left(x,t\right),x\in \Omega ,t\in \left(0,T\right],\end{array}$

其边界条件为Dirichlet条件。该问题的精确解为：

$\begin{array}{l}u\left(x,t\right)=\mathrm{sech}\left(5x+5t+2\right)\mathrm{sech}\left(5y+5t+2\right){\text{e}}^{\text{i}\left(2x+2y-3t\right)},\\ v\left(x,t\right)=\mathrm{sech}\left(5x-5t-2\right)\mathrm{sech}\left(5y-5t-2\right){\text{e}}^{\text{i}\left(x+y-t\right)},\end{array}$

(a) t = 0 (b) t = 0.1 (c) t = 0.3

(d) t = 0.4 (e) t = 0.55 (f) t = 0.7

(g) t = 0.8 (h) t = 0.9 (i) t = 1

Figure 6. Profiles of the numerical solution at different times

图6. 数值解在不同时刻的剖面图

(a) t = 0 (b) t = 0.1 (c) t = 0.3

(d) t = 0.4 (e) t = 0.55 (f) t = 0.7

(g) t = 0.8 (h) t = 0.9 (i) t = 1

Figure 7. Contour plots of numerical solution at different times

图7. 数值解在不同时刻的等高线图

考虑初值

$\begin{array}{l}{u}^0\left(x\right)=\mathrm{sech}\left(5x+2\right)\mathrm{sech}\left(5y+2\right){\text{e}}^{\text{i}\left(2x+2y\right)},\\ v^0\left(x\right)=\mathrm{sech}\left(5x-2\right)\mathrm{sech}\left(5y-2\right){\text{e}}^{\text{i}\left(x+y\right)}.\end{array}$

在我们的计算中，计算区域被定义为$-2\le x,y\le 2$ 的矩形区域，在空间中分布160,000个配置节点，包括1596个边界点以及158,404个内点，该解被定义为在该区间之外为零。该解从$t=0$ 模拟到$t=1$ ，时间步长设为$\tau =0.01$ 。图6和图7显示了两个孤子的传播过程。两个振幅不同的孤子从$t=0$ 开始向对方移动。大约在$t=0.55$ 时刻，两个孤子开始碰撞，它们逐渐变形，最后合并成一个波包，这个波包碰撞后很快分裂成两个孤子。在$0\le t\le 1$ 范围内，函数$u$ 的解析解与数值解的最大误差为$3.6654\times {10}^{-4}$ ，函数$v$ 的解析解与数值解的最大误差为$1.9861\times {10}^{-3}$ 。在我们实验模拟中，亮孤子在二维传播过程中，其中心位置、峰值幅度以及空间宽度基本保持不变。在适当的初始条件与非线性参数设定下，亮孤子在二维空间中可以实现稳定的自维持传播。

5. 结论

本文提出了一种求解CNLS的GFDM与线性化格式相结合的数值方法，并进行数值计算验证了所提出数值格式的收敛性。误差结果表明，GFDM对求解不规则区域中的非线性耦合薛定谔方程组是准确有效的。此外，为进一步增强处理此类问题的精度和效率，本文基于GFDM的无网格方法进行了数值模拟，模拟效果表明所提出的格式在演化孤子传播中能保持稳定性和可靠性，这对CNLS的应用研究具有重要意义。

NOTES

^*通讯作者。