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

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

中圖分類號： O 241.82文獻(xiàn)標(biāo)志碼： 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

在經(jīng)典力學(xué)中，質(zhì)點的狀態(tài)采用質(zhì)點的坐標(biāo)和速度進(jìn)行描述，質(zhì)點的運動方程就是牛頓運動方程。而在量子力學(xué)中，微觀粒子的狀態(tài)則采用波函數(shù)進(jìn)行描述，且決定粒子狀態(tài)變化的方程不再是牛頓運動方程，而是Schr dinger方程 ［1］。1926年，奧地利物理學(xué)家薛定諤提出了Schr dinger方程，該方程是量子力學(xué)領(lǐng)域的基本方程，其在量子力學(xué)的重要性毫不遜色于牛頓運動定律在經(jīng)典力學(xué)中的重要性.非線性Schr dinger（NLS）方程在量子力學(xué)、非線性光學(xué)、電磁學(xué)等離子理論和固體物理等領(lǐng)域中有著廣泛的應(yīng)用。帶初邊值條件的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為虛數(shù)單位；ρ，β為實常數(shù)；Ω∈Rd（d=1，2）是有界區(qū)域；Δ為Laplace算子，復(fù)值函數(shù)u（x，t）為波函數(shù)；實值函數(shù)v（x）為外場的勢能；u0（x）為足夠光滑的函數(shù)。

帶初邊值條件的NLS方程的計算結(jié)果很好地反映量子力學(xué)效應(yīng)和微觀系統(tǒng)性質(zhì)，且能很好地描述微觀粒子的狀態(tài)隨時間的變化情況。NLS方程滿足質(zhì)量守恒，即

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方程的數(shù)值研究引起了廣大學(xué)者的關(guān)注。Hu等［2］研究了帶5次項的NLS方程4階緊致差分格式，并證明最大范數(shù)下的無條件穩(wěn)定性和收斂性。Guo等［3］ 利用兩級有限差分格式結(jié)合吸收邊界方法求NLS方程。Feng 等［4］構(gòu)造NLS方程的高階守恒SAV-Gauss配置有限元格式。Wang等［5］提出NLS方程的兩層網(wǎng)格有限元格式，并對其進(jìn)行超收斂性分析。Fu等［6］研究二維NLS方程的顯式高階保指數(shù)差分龍格庫塔格式。Hu等［7］提出帶波動算子的二維NLS方程的守恒型差分格式。Zhai等［8］利用算子分裂法求解空間分?jǐn)?shù)階NLS方程，并進(jìn)行嚴(yán)格的誤差分析和數(shù)值模擬。Deng等［9］ 提出NLS方程的二階SAV格式，并給出嚴(yán)格的誤差分析。

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

1 預(yù)備知識

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

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

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

為了保證良好的數(shù)值穩(wěn)定性，重心Lagrange插值逼近未知函數(shù)時使用第2類Chebyshev節(jié)點，即

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

p（x）在節(jié)點xj處的v階導(dǎo)數(shù)可以表示為

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

由文獻(xiàn)［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方程的半離散格式進(jìn)行推導(dǎo)，設(shè)xi（i=0，1，…，m）為Chebyshev空間節(jié)點，時間節(jié)點為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階龍格庫塔方法進(jìn)行離散，得到的全離散格式為

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 相容性分析

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

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

引理1［21］ 假設(shè)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是自然常數(shù)。

基于引理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。

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

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）

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

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 數(shù)值算例

定義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處的數(shù)值解。

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

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

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

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

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

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

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

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

算例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。

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

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

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

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

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

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

5 結(jié)束語

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

參考文獻(xiàn)：

［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］ 胡玉蝶，彭澳，陳麗權(quán)，等.有限差分-配點法求解二維Burgers方程［J］.純粹數(shù)學(xué)與應(yīng)用數(shù)學(xué)，2023，39（1）：100-112.DOI：10.3969/j.issn.1008-5513.2023.01.008.

［13］ 羅詩棟，凌永輝.Rosenau-Burgers方程的一種高精度有限差分格式［J］.閩南師范大學(xué)學(xué)報（自然科學(xué)版），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］.應(yīng)用數(shù)學(xué)和力學(xué)，2023，44（5）：573-582.DOI：10.21656/1000-0887.430149.

［18］ 黃蓉，翁智峰.時間分?jǐn)?shù)階Allen-Cahn方程的重心插值配點法［J］.華僑大學(xué)學(xué)報（自然科學(xué)版），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］ 黃蓉，姚夢麗，翁智峰.對流擴(kuò)散方程最優(yōu)控制問題的重心插值配點格式［J］.華僑大學(xué)學(xué)報（自然科學(xué)版），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.

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