非線性Schr dinger方程的龍格庫塔配點格式

摘要： 采用4階龍格庫塔方法結合重心Lagrange插值配點法求解非線性Schr dinger方程。首先，在空間方向上采用重心Lagrange插值配點格式進行離散，在時間方向上采用4階龍格庫塔方法離散，從而得到非線性Schr dinger方程的龍格庫塔配點格式。其次，對全離散格式進行相容性分析。結果表明：龍格庫塔配點格式具有高精度，且滿足離散的質量和能量守恒。

關鍵詞： 非線性Schr dinger方程； 4階龍格庫塔方法； 重心Lagrange插值配點； 相容性分析

中圖分類號： O 241.82文獻標志碼： A"" 文章編號： 1000 5013（2024）04 0534 09

Runge-Kutta Collocation Scheme for Nonlinear Schr dinger Equation

YAO Mengli， TENG Yuhang， LAI Yiying， WENG Zhifeng

（School of Mathematical Sciences， Huaqiao University， Quanzhou 362021， China）

Abstract： The fourth order Runge-Kutta method and barycentric Lagrange interpolation collocation method are used to solve the nonlinear Schr dinger equation. Firstly， the barycentric Lagrange interpolation collocation scheme is discreted in the spatial direction， and the fourth-order Runge-Kutta method is discreted in the temporal direction. The Runge-Kutta collocation scheme of the nonlinear Schr dinger equation is obtained. Secondly， the consistency analysis of the fully discrete scheme is analyzed. The results show that the Runge-Kutta collocation scheme has the high accuracy and satisfies the conservation of discrete mass and energy.

Keywords： nonlinear Schr dinger equation； fourth-order Runge-Kutta method； barycentric Lagrange interpolation collocation； consistency analysis

在經典力學中，質點的狀態采用質點的坐標和速度進行描述，質點的運動方程就是牛頓運動方程。而在量子力學中，微觀粒子的狀態則采用波函數進行描述，且決定粒子狀態變化的方程不再是牛頓運動方程，而是Schr dinger方程 ［1］。1926年，奧地利物理學家薛定諤提出了Schr dinger方程，該方程是量子力學領域的基本方程，其在量子力學的重要性毫不遜色于牛頓運動定律在經典力學中的重要性.非線性Schr dinger（NLS）方程在量子力學、非線性光學、電磁學等離子理論和固體物理等領域中有著廣泛的應用。帶初邊值條件的NLS方程為

iut+ρΔu+v（x）u+βu2u=0，" x∈Ω， t∈［0，T］，u（x，0）=u0（x），" x∈Ω，u（x，t）=0，" x∈Ω， t∈［0，T］。（1）

式（1）中：i為虛數單位；ρ，β為實常數；Ω∈Rd（d=1，2）是有界區域；Δ為Laplace算子，復值函數u（x，t）為波函數；實值函數v（x）為外場的勢能；u0（x）為足夠光滑的函數。

帶初邊值條件的NLS方程的計算結果很好地反映量子力學效應和微觀系統性質，且能很好地描述微觀粒子的狀態隨時間的變化情況。NLS方程滿足質量守恒，即

M（t）=∫Ω|u（x，t）|2dx=M（0），（2）

以及滿足能量守恒，即

E（t）=∫Ω（ρ|SymbolQC@u（x，t）|2-v（x）|u（x，t）|2-β2|u（x，t）|4）dx=E（0）。（3）

近年來，對NLS方程的數值研究引起了廣大學者的關注。Hu等［2］研究了帶5次項的NLS方程4階緊致差分格式，并證明最大范數下的無條件穩定性和收斂性。Guo等［3］ 利用兩級有限差分格式結合吸收邊界方法求NLS方程。Feng 等［4］構造NLS方程的高階守恒SAV-Gauss配置有限元格式。Wang等［5］提出NLS方程的兩層網格有限元格式，并對其進行超收斂性分析。Fu等［6］研究二維NLS方程的顯式高階保指數差分龍格庫塔格式。Hu等［7］提出帶波動算子的二維NLS方程的守恒型差分格式。Zhai等［8］利用算子分裂法求解空間分數階NLS方程，并進行嚴格的誤差分析和數值模擬。Deng等［9］ 提出NLS方程的二階SAV格式，并給出嚴格的誤差分析。

以上求解NLS方程的數值方法都是基于網格剖分，從而建立逼近格式。近年來，一種新型的無網格方法（重心插值配點法）引起學者關注。重心插值配點法具有計算格式簡單、精度高、程序實施方便和節點適應性好等特點，使用Chebyshev節點的重心Lagrange插值公式還可以有效克服Runge現象。目前，重心插值配點法已經被推廣到求解Sine-Gordon方程［10］、Burgers方程［11-13］、粘彈性波方程［14］、Allen-Cahn方程［15-18］、非線性對流擴散最優控制問題［19-20］和分數階電報方程［21］等?；谇叭说墓ぷ?，本文主要采用4階龍格庫塔和重心Lagrange插值配點格式相結合的方法求解NLS方程，并給出該問題全離散格式的相容性分析。

1 預備知識

給定m+1個節點xj和其函數值gj（j=0，1，…，m），設插值多項式p（x）是g（x）的近似值，且滿足p（xj）=gj。將p（x）改寫成Lagrange插值多項式，即

g（x）≈p（x）=∑mj=0Lj（x）gj。（4）

式（4）中：Lj（x）為Lagrange插值基函數，且Lj（x）=∏mi=0，i≠j（x-xi）/∏mi=0，i≠j（xj-xi）。

為了保證良好的數值穩定性，重心Lagrange插值逼近未知函數時使用第2類Chebyshev節點，即

xj=cosjmπ，" j=0，1，…，m。（5）

p（x）在節點xj處的v階導數可以表示為

p（v）（xi）=dvp（xi）dxv=∑mj=0L（v）j（xi）pj=∑mj=0D（v）i，jpj，" v=1，2，…，m。（6）

由文獻［22］可知，二階微分矩陣D（2）的計算公式為

D（2）i，j=L″j（xi）=-2wj/witi-tj∑k≠iwk/witi-tk+1ti-tj，" j≠i，D（2）i，i=-∑mj=0，j≠iD（2）i，j。（7）

2 NLS方程的龍格庫塔配點格式

對一維NLS方程的半離散格式進行推導，設xi（i=0，1，…，m）為Chebyshev空間節點，時間節點為tk=kτ，k=0，1，2，…，n，τ=Tn。

在空間方向上采用重心Lagrange插值配點法得到半離散格式，即

iut（xi，t）=-ρ∑mk=0L″i（x）u（xi，t）-v（xi）u（xi，t）-βu（xi，t）2u（xi，t）。（8）

式（8）的矩陣形式為

i（uh（t））t=-ρ（D（2）Im）uh（t）-Vuh（t）-βuh（t）2uh（t）。（9）

式（9）中：符號表示矩陣的Kronecker積；uh（t）=［u0（t），u1（t），…，um（t）］′；V=diag［v（x0），v（x1），…，v（xm）］。

類似可得二維NLS方程的空間半離散格式，即

i（uh（t））t=-ρ（（C（2）Im）uh（t）+（InD（2））uh（t））-Vuh（t）-βuh（t）2uh（t）。（10）

定義算子A為

A=i［ρ（C（2）Im）+ρ（InD（2））+V+diag（βuh（t）2）。（11）

式（10）可改寫為

（uh（t））t=Auh（t）。（12）

其次，令ukh=uh（tk），在時間方向上用4階龍格庫塔方法進行離散，得到的全離散格式為

uk+1h=ukh+τ6（k1+2k2+2k3+k4）。（13）

式（13）中：k1=Aukh；k2=Aukh+k1τ2；k3=Aukh+k2τ2；k4=A（ukh+k3τ）。

3 相容性分析

設p（x，y）是u（x，y）的拉格朗日插值函數，滿足p（xi，yj）=u（xi，yj），i=0，1，…，m；j=0，1，…，n。 定義誤差函數為

r（x，y）=u（x，y）-p（x，y）。（14）

引理1［21］ 假設u∈Cm+1（a，b），則有

‖r（x，y）‖∞≤C1‖u（m+1）‖∞elx2mm+C2‖u（n+1）‖∞ely2nn，‖rx，x（x，y）‖∞≤C**1‖u（m+1）‖∞elx2（m-2）m-2+C2‖u（n+1）‖∞ely2nn，‖ry，y（x，y）‖∞≤C1‖u（m+1）‖∞elx2mm+C**2‖u（n+1）‖∞ely2（n-2）（n-2）。

式中：lx=b-a2；ly=d-c2；e是自然常數。

基于引理1，可得定理1。

定理1［23］若u（x，y，t）∈Cm+1（Ω）×C2（0，T］，Ω=［a，b］×［c，d］，則有

‖u（x，y，t）-u（xi，yj，t）‖∞≤C**1‖u（m+1）‖∞elx2（m-2）m-2+C**2‖u（n+1）‖∞ely2（n-2）n-2。

定理2為時間方向上基于4階龍格庫塔方法的全離散格式相容性分析。

定理2 若u（x，y，t）∈Cm+1（Ω）×C2（0，T］，Ω=［a，b］×［c，d］，則有

‖u（x，y，t）-u（xi，yj，tk+1）‖∞≤Cτ5+‖u（m+1）‖∞elx2（m-2）m-2+‖u（n+1）‖∞ely2（n-2）n-2。

證明：設u（xi，yj，t）是u（x，y，t）空間方向基于重心插值配點法離散的數值解，則

iut（xi，yj，t）=-ρΔu（xi，yj，t）-v（xi，yj）u（xi，yj，t）-βu（xi，yj，t）2u（xi，yj，t）+γi，j。（15）

令Th（xi，yj，t）=-ρΔu（xi，yj，t）-v（xi，yj）u（xi，yj，t）-βu（xi，yj，t）2u（xi，yj，t），（16）

則式（15）可簡化為

iut（xi，yj，t）=Thu（xi，yj，t）+γi，j。（17）

式（17）中：γi，j是空間截斷誤差。

由定理1可知

γi，j≤C**1‖u（m+1）‖∞elx2（m-2）m-2+C**2‖u（n+1）‖∞ely2（n-2）n-2。（18）

設uk+1h=u（xi，yj，tk+1）是u（x，y，t）時間方向基于4階龍格庫塔方法離散的數值解，則由泰勒展開公式有

uk+1h=ukh+τ6（k1+2k2+2k3+k4）+γi，j+O（τ5）。（19）

式（19）中：k1=Thukh；k2=Thukh+k1τ2；k3=Thukh+k2τ2；k4=Th（ukh+k3τ）。

證明完畢。

4 數值算例

定義L∞誤差（Err∞）和L2誤差（Err2）分別為

Err∞=max1≤i，j≤MUi，j-ui，j，Err2=h∑Mi，j=1（Ui，j-ui，j）2。

式中：ui，j表示點xi，j處的精確解；Ui，j表示點xi，j處的數值解。

算例1 ρ=12，β=-1，式（1）的真解為u（x，t）=exp（-3it2）·sin x，初始條件為u（x，0）=sin x，邊界條件為u（0，t）=u（2π，t）=0，勢函數為v（x）=-cos2x。

選取區域Ω=［0，2π］×［0，1］，時間步長τ=10-4。重心Lagrange插值配點格式和2階有限差分格式的L∞誤差和L2誤差，如表1所示。表1中：M為節點數。由表1可知：重心Lagrange插值配點格式在空間上用較少的點就可以達到更高的精度。

重心插值配點格式的誤差及時間收斂階，如表2所示。表2中：Rate1，Rate2分別為L∞誤差收斂階，L2誤差收斂階。由表2可知：龍格庫塔配點格式的時間收斂階是4階。

令M=48，重心Lagrange插值配點格式的數值解、精確解（M=48），如圖1所示。由圖1可知：數值解圖像與精確解圖像逼近。

不同數值格式的空間收斂階，如圖2所示。由圖2可知：有限差分格式的收斂階是2階，而重心Lagrange插值配點格式的收斂階滿足指數收斂的性質，明顯優于經典的有限差分方法。

重心Lagrange插值配點格式對守恒量的保持情況（M=48），如圖3所示。圖3中：M（t），E（t）分別為質量、能量函數。由圖3可知：龍格庫塔配點格式滿足質量守恒和能量守恒，與理論相符。

算例2 考慮一維NLS方程，選取M=138，τ=1/40 000，對應的計算區間為［-2π，10π］×［0，5］。ρ=1，β=2，選取初始條件為u（x，0）=sech xexp（2ix），邊界條件為u（a，t）=u（b，t）=0，勢函數v（x）=0。

重心Lagrange插值配點格式的孤立波波形，如圖4所示。由圖4可知：龍格庫塔配點格式的孤立波波形隨時間的改變而不斷演化，并且在演化的過程中波形并沒有出現絲毫震蕩，模擬結果充分體現出該格式的穩定性。

算例3 考慮如下二維NLS方程，即

iut+12Δu-v（x，y）u-u2u=0，u（x，y，0）=sin xsin y，u（0，y，t）=u（2π，y，t）=0，u（x，0，t）=u（x，2π，t）=0。

勢函數v（x，y）=1-sin2xsin2y，真解u（x，y，t）=exp（-2it）·sin xsin y。

選取區域Ω×［0，1］，Ω=［0，2π］2，τ=10-3。重心插值配點格式和二階有限差分格式的L∞誤差和L2誤差（τ=10-3），如表3所示。由表3可知：重心Lagrange插值配點格式在空間上用較少的點就可達到更高的精度，計算精度明顯優于經典的有限差分方法。

重心插值配點格式的誤差及時間收斂階（M=16），如表4所示。由表4可知：龍格庫塔配點格式的時間收斂階是4階。

令M=N=40，重心Lagrange插值配點格式的數值解、精確解（M=40），如圖5所示。由圖5可知：數值解圖像和精確解圖像基本吻合。

不同數值格式的空間收斂階，如圖6所示。

重心Lagrange插值配點格式對守恒量的保持情況（M=40），如圖7所示。

5 結束語

將龍格庫塔與重心Lagrange 插值Chebyshev 配點法相結合，時間精度可達到4階，空間精度滿足指數收斂，給出全離散格式的相容性分析，通過數值算例驗證所提格式的高精度和有效性。通過與經典的差分格式進行比較，表明提出的格式可以用較少的點達到較高的精度。該方法可推廣到其他非線性微分方程，從而為解決同類問題提供一種高精度數值求解方法。

參考文獻：

［1］ SCHRDINGER E.The present status of quantum mechanics［J］.Die Naturwissenschaften，1935，23：1-26.DOI：10.48550/arXiv.2104.09945.

［2］ HU Hanqing，HU Hanzhang.Maximum norm error estimates of fourth-order compact difference scheme for the nonlinear Schr dinger equation involving a quintic term［J］.Journal of Inequalities and Applications，2018，2018（1）：1-15.DOI：10.1186/s13660-018-1775-y.

［3］ GUO Feng，DAI Weizhong.A new absorbing layer approach for solving the nonlinear Schr dinger equation［J］.Applied Numerical Mathematics，2023，189：88-106.DOI：10.1016/j.apnum.2023.04.003.

［4］ FENG Xiaobing，LI Buyang，MA Shu.High-order mass-and energy-conserving SAV-Gauss collocation finite element methods for the nonlinear Schr dinger equation［J］.SIAM Journal on Numerical Analysis，2021，59：1566-1591.DOI：10.1137/20M1344998.

［5］ WANG Junjun，LI Meng，GUO Lijuan.Superconvergence analysis for nonlinear Schr dinger equation with two-grid finite element method［J］.Applied Mathematics Letters，2021，122：107553.DOI：10.1016/j.aml.2021.107553.

［6］ FU Yayun，XU Zhuangzhi.Explicit high-order conservative exponential time differencing Runge-Kutta schemes for the two-dimensional nonlinear Schr dinger equation［J］.Computers and Mathematics with Applications，2022，119：141-148.DOI：10.1016/j.camwa.2022.05.021.

［7］ HU Hanzhang，CHEN Yanping.A conservative difference scheme for two dimensional nonlinear Schr dinger equation with wave operator［J］.Numerical Methods for Partial Differential Equations，2016，32（3）：862-876.DOI：10.1002/num.22033.

［8］ ZHAI Shuying，WANG Dongling，WENG Zhifeng，et al.Error analysis and numerical simulations of strang splitting method for space fractional nonlinear Schr dinger equation［J］.Journal of Scientific Computing，2019，81：965-989.DOI：10.1007/s10915-019-01050-w.

［9］ DENG Beichuan，SHEN Jie，ZHUANG Qingqu.Second-order SAV schemes for the nonlinear Schr dinger equation and their error analysis［J］.Journal of Scientific Computing，2021（69）：88.DOI：10.1007/s10915-021-01576-y.

［10］ LI Jin，QU Jinzheng.Barycentric Lagrange interpolation collocation method for solving the Sine-Gordon equation［J］.Wave Motion，2023，120：103159.DOI：10.1016/j.wavemoti.2023.103159.

［11］ HU Yudie，PENG Ao，CHEN Liquan，et al.Analysis of the barycentric interpolation collocation scheme for the Burgers equation［J］.Science Asia，2021，47：758-765.DOI：10.2306/ scienceasia1513-1874.2021.081.

［12］ 胡玉蝶，彭澳，陳麗權，等.有限差分-配點法求解二維Burgers方程［J］.純粹數學與應用數學，2023，39（1）：100-112.DOI：10.3969/j.issn.1008-5513.2023.01.008.

［13］ 羅詩棟，凌永輝.Rosenau-Burgers方程的一種高精度有限差分格式［J］.閩南師范大學學報（自然科學版），2022，35（4）：5-12.DOI：10.16007/j.cnki.issn2095-7122.2022.04.013.

［14］ ORUC .Two meshless methods based on local radial basis function and barycentric rational interpolation for solving 2D viscoelastic wave equation［J］.Computational and Applied Mathematics，2020，79：3272-3288.DOI：10.1016/j.camwa.2020.01.025.

［15］ DENG Yangfang，WENG Zhifeng.Barycentric interpolation collocation method based on Crank-Nicolson scheme for the Allen-Cahn equation［J］.AIMS Mathematics，2021，6：3857-3873.DOI：10.3934/math.2021229.

［16］ DENG Yangfang，WENG Zhifeng.Operator splitting scheme based on barycentric Lagrange interpolation collocation method for the Allen-Cahn equation［J］.Journal of Applied Mathematics and Computing，2022，68（5）：3347-3365.DOI：10.1007/s12190-021-01666-y.

［17］ 黃蓉，鄧楊芳，翁智峰.SAV/重心插值配點法求解Allen-Cahn方程［J］.應用數學和力學，2023，44（5）：573-582.DOI：10.21656/1000-0887.430149.

［18］ 黃蓉，翁智峰.時間分數階Allen-Cahn方程的重心插值配點法［J］.華僑大學學報（自然科學版），2022，43（4）：553-560.DOI：10.11830/ISSN.1000-5013.202101060.

［19］ HUANG Rong，WENG Zhifeng.A numerical method based on barycentric interpolation collocation for nonlinear convection-diffusion optimal control problems［J］.Networks and Heterogeneous Media，2023，18（2）：562-580.DOI：10.3934/nhm.2023024.

［20］ 黃蓉，姚夢麗，翁智峰.對流擴散方程最優控制問題的重心插值配點格式［J］.華僑大學學報（自然科學版），2023，44（3）：407-416.DOI：10.11830/ISSN.1000-5013.202203023.

［21］ YI Shichao，YAO Linquan.A steady barycentric Lagrange interpolation method for the 2D higher-order time-fractional telegraph equation with nonlocal boundary condition with error analysis［J］.Numerical Methods for Partial Differential Equations，2019（35）：1694-1716.DOI：10.1002/num.22371.

［22］ BERRUT J P，TREFETHEN L N.Barycentric Lagrange interpolation［J］.SIAM Review，2004，46：501-507.DOI：10.1137/S0036144502417715.

［23］ SUN Haoran，HUANG Siyu，ZHOU Mingyang，et al.A numerical investigation of nonlinear Schr dinger equation using barycentric interpolation collocation method［J］.AIMS Mathematics，2023，8（1）：361-381.DOI：10.3934/math.2023017.

（責任編輯：" 陳志賢" 英文審校： 黃心中）