基于徑向基函數的局部近似特解法求解二維薛定諤方程

洪永興　陳　文　林　繼

(河海大學 力學與材料學院， 南京　211100)

?

基于徑向基函數的局部近似特解法求解二維薛定諤方程

洪永興陳文林繼

(河海大學 力學與材料學院， 南京211100)

摘要：基于徑向基函數的局部近似特解法具有形式簡單、易編程、精度高、收斂速度快等優點，是一種純無網格配點方法．它通過將計算域劃分為若干子區域并利用各個區域內的節點構造局部低階矩陣，然后再將該矩陣拓展為全局形式，從而構造一個全局稀疏矩陣，以便于快速計算．本文采用局部近似特解法數值模擬二維薛定諤方程，首先采用隱式歐拉差分格式對時間項進行離散，并利用基于Multiquadric(MQ)函數的局部近似特解法對空間項進行離散．數值實驗表明，局部近似特解法求解精度高、收斂速度快且計算耗時少，具有較好的工程應用前景．

關鍵詞：隱式歐拉差分；二維薛定諤方程；Multiquadric函數；局部近似特解法；全局近似特解法

薛定諤方程是重要的波動方程之一，廣泛存在于光學、生物學、海洋學等領域[1-3]．國內外許多學者對薛定諤方程進行了深入的研究[4-7]．但是因為薛定諤方程形式的復雜性，解析求解變得相當困難，所以尋找一個簡單、快速、有效的數值方法具有重要的理論意義和應用價值．

徑向基函數配點法作為一種重要的無網格方法，被廣泛用于研究科學和工程問題[8-12]．該方法具有精度高，編程容易，思路簡單等優點．但是，采用全局的徑向基函數配點方法在求解過程中可能產生病態稠密矩陣，影響計算結果，并且計算耗時長阻礙了其在大規模科學和工程問題中的應用．為了克服該缺點，學者們提出了局部化的徑向基函數方法，此類方法能快速高效地處理大規模科學和工程問題．目前，基于徑向基函數的局部方法有以下幾種：局部緊支徑向基函數方法[13]、局部徑向基函數配點法[14]、基于徑向基函數的局部微分求積法[15]和局部近似特別解法[16]等．其中局部近似特別解法最早由Chen和Yao提出，用于心血管內鈣元素動力學行為模擬[17]等大規模物理力學問題．隨后局部近似特解法得到學者們的關注，被用于波和熱傳導問題的數值模擬[18-20]．

本文采用局部近似特解法數值模擬二維薛定諤方程，首先采用隱式歐拉差分法對時間項進行離散，然后采用基于MQ函數的局部近似特解法對空間項進行離散，從而離散整個薛定諤方程．通過與全局近似特解法進行系統的比較發現，局部近似特解法求解精度高、收斂速度快且計算耗時少．

1薛定諤方程及時間項離散

1.1二維薛定諤方程

本文考慮如下簡單形式的二維薛定諤方程：

(1)

考慮如下初始條件和邊界條件：

(2)

(3)

1.2時間項離散

首先采用隱式歐拉差分格式對式(1)時間項進行離散，可得

(4)

其中dt為時間步長dt=tn+1-tn，tn=n·dt，以及un=(x，y，tn)，Δun=?2u(x，y，tn)/?x2+?2u(x，y，tn)/?y2．

將包含un+1和un的項分別整理到等號的左邊和右邊，可得

(5)

2近似特解法與局部近似特解法

2.1近似特解法

不失一般性，考慮如下的偏微分方程：

(6)

(7)

其中未知數x=(x1，x2)，u(x)為待求的函數，Ω為問題的求解區域，?Ω為區域邊界滿足Dirichlet邊界條件，f(x)和g(x)為已知函數．那么未知函數u(x)可以表示為徑向基函數的線性組合，如下：

(8)

(9)

(10)

其中，φ=ΔΦ，ni為內部點數，nb=n-ni為邊界點數，n為總點數．由此可以計算出式(8)的未知系數{αj}．

近似特解法和Kansa方法的主要區別在于：近似特解法是對徑向基函數φ積分求得Φ的表達式，而Kansa方法是對Φ進行求導得到φ的具體表達式[21]．本文采用MQ函數作為徑向基函數，分別采用特解法和Kansa方法推導得到如下Φ和φ：

Kansa方法：

特解法：

其中c為形狀參數．需要指出的是，采用全局的近似特解法(Globalmethodofapproximateparticularsolution，GMAPS)將可能產生病態稠密矩陣，會增加計算耗時并影響計算精度，難以處理大規模科學和工程問題[8]．

2.2局部近似特解法

局部近似特解法(Localizedmethodofapproximateparticularsolution，LMAPS)可以從一定程度上減少病態稠密矩陣對結果的影響．它將全局空間劃分成若干子區域，對每個區域利用徑向基函數進行插值．一般構造影響區域的方式是以某一插值點為中心尋找該插值點鄰近的ns個點作為它影響區域，稱ns為區域值．下面簡單介紹局部近似特解法的基本步驟，由式(8)可得

(11)

由式(9)～(11)可得

(12)

其中uhI=[f(x1) f(x2) … f(xni)]T，uhB=[g(xni+1) g(xni+2) … g(xn)]T，AI=[φij]ni×n，AB=[Φij]nb×n．那么

(13)

接著在每個子區域Ωc內，采用近似特解法構造低階矩陣：

(14)

其中Ωc，Ac分別表示以xc為中心的影響區域和該區域的滿秩矩陣．

由式(14)和(13)可得

(15)

(16)

其中

此外局部近似特解法中的形狀參數不是常數，而是與距離相關的變量

(17)

其中dr是區域中心點到各點的最大距離，c為常數．

3數值結果和討論

本節采用局部近似特解法求解二維薛定諤方程．為了驗證該方法的有效性，本文求解兩種區域問題：規則和不規則區域．誤差分析采用標準誤差(Rootmeansquareerror， RMS)，相對誤差(Relativeerror， RE)和最大相對誤差(Maximumrelativeerror， MRE)，定義如下：

MRE=max(RE(xk))k=1，2，…，M

式中，M為測試點總數．其中RMSr表示實部標準誤差，RMSi表示虛部標準誤差，REr表示實部相對誤差，REi表示虛部相對誤差．

考慮如下二維薛定諤方程：

(18)

初始條件u(x，y，0)和邊界條件u(x，y，t)分別為：

(19)

(20)

用于比較的解析解為：

(21)

3.1算例1：規則區域

首先考慮如圖 1所示的規則區域，其中邊界點數nb=108，內部點數ni=415，總點數n=523．

圖1　規則區域及配點示意圖

圖2給出了局部近似特解法和全局近似特解法的誤差，其中dt=0.01，ns=13，c=15(LMAPS)，c=0.1(GMAPS)．數值結果表明兩者均能較好地模擬二維薛定諤方程，并且LMAPS方法的數值精度在某些時刻較GMAPS高．

由文獻可知，局部近似特解法采用與節點距離有關的形狀參數可以有效地降低形狀參數c對精度的影響[20]．圖3顯示了誤差與形狀參數c之間的關系．從圖3可以看出，局部近似特解法形狀參數c的有效取值范圍比全局近似特解法的大得多，該結果驗證了上述結論．

圖2　LMAPS和GMAPS誤差比較

圖3　誤差與形狀參數c的關系，dt=0.01，T=40，ns=13

圖4給出了局部近似特解法的誤差隨著ns的變化曲線．從圖中可以看出，適當增加影響區域的插值點數可以提高精度．這是因為增加影響區域的插值點數可以提高控制方程的條件數[22]．雖然此時同時增加計算量，但是由于區域插值點數總小于總點數，該方法的計算量總比全局近似特解法的小．當ns增大到總點數時，兩者的計算量相同．

圖4　誤差隨ns的變化曲線dt=0.01，c=13

其次，分別減小空間步長dx，dy和時間步長dt，分析局部近似特解法和全局近似特解法的穩定性．圖5表明局部近似特解法的收斂速度快，且可以通過采用較少的配點，保證同等的精度要求，因此可以進一步提高計算效率．圖6表明隨著時間步長dt的減小，局部近似特解法可收斂到一個更高的精度．圖5和圖6同時說明，局部近似特解法具有更高的穩定性．

圖5　誤差隨dx，dy的變化曲線dt=0.01，c和ns取dx，dy對應的最優值

圖6　誤差隨dt變化趨勢ns=13，c=13

下面簡要分析局部近似特解法和全局近似特解法的計算耗時．將計算分為兩個階段，零時刻的賦值階段和時間方向的離散階段．首先第一階段，后者的矩陣需要賦值數為ni·ni+nb，前者的為ns·ni+nb，由于ns?ni，因此前者的計算耗時比后者的少；其次是第二階段，由于前者采用稀疏矩陣，在時間方向上每一步的計算耗時均比后者少．因此，局部近似特解法的計算耗時比全局近似特解法的少，圖7驗證了上述結論．

圖7　計算耗時ns=13

3.2算例2：不規則區域

為進一步驗證該方法的有效性，考慮如圖8的不規則區域，其中邊界點nb=136，內部點ni=625，總點數n=761．圖9給出了對應不同時刻的相對誤差曲面圖(c=22，ns=13，dt=0.001)，圖10～12顯示了三個隨機配點處模擬值與真實值的吻合程度，圖13顯示了布點數增加到10萬時的計算耗時，數值結果表明局部近似特解法具有較高的精度和時間穩定性，且計算耗時少，可以有效地求解此類二維薛定諤方程．

圖8　不規則區域及配點示意圖

圖9　不同時間T對應的相對誤差曲面圖

圖10　隨機配點(0.175，0.720)

圖11　隨機配點(0.130，0.275)

圖12　隨機配點(0.945，0.570)

圖13　計算耗時ns=13

4結語

本文采用基于MQ函數的局部近似特解法，模擬二維薛定諤方程．上文的算例表明，該方法能很好地解決這類問題，并且克服了全局方法對形狀參數的敏感性和稠密矩陣帶來的求解困難等問題．同時利用稀疏矩陣，大幅度提高了求解時域問題的計算效率．總之，該方法很好地繼承了全局近似特解法的優點，同時改善了它的部分缺點，并在一定程度上提高了它的穩定性．鑒于該方法無網格的優點，可以將它推廣到三維和大規模工程問題上[23]．對于三維問題，其矩陣將更加稠密．因此，需要對本文所提方法進行改進，如引進更有效的形狀參數c，引入多尺度徑向基函數等．

參考文獻：

[1]宋祥，李畫眉．變系數耦合非線性薛定諤方程的矢量孤子解：暗-亮孤子解[J]．浙江師范大學學報：自然科學版，2013(3)：294-298．

[2]宋詩艷，王晶，孟俊敏，等．深海內波非線性薛定諤方程的研究[J]．物理學報，2010(2)：1123-1129．

[3]李宏，孟祥東，張斯淇，等．用薛定諤方程和路徑積分方法研究中子雙縫衍射[J]．南開大學學報：自然科學版，2014(3)：68-73．

[4]劉登云，王劍波．用因式分解法求解薛定愕方程[J]．大學物理，1990(7)：13-15．

[5]宮建平．有限差分法求解薛定諤方程[J]．晉中學院學報，2014，31(3)：1-6．

[6]Tang D J． Generalized Matrix Numerov Solutions to the Schr?dinger equation[D]． Singapore：National University of Singapore，2014．

[7]Filiz A，Ekici M，Sonmezoglu A． F-Expansion Method and New Exact Solutions of the Schr?dinger-KdV Equation [J]．The Scientific World Journal，2014：1-14．

[8]王嬋媛，石記松，陳文．基于徑向基函數的Hermite配點法的薄膜自振分析[A]．中國土木工程學會教育工作委員會．第六屆全國土木工程研究生學術論壇論文集[C]．中國土木工程學會教育工作委員會，2008：1．

[9]Kansa E J． Multiquadrics-A Scattered Data Approximation Scheme with Applications to Computational Fluid-dynamics-I Surface Approximations and Partial Derivative Estimates[J]．Computers & Mathematics with Applications，1990，19(8-9)：147-161．

[10] Duan Y，Hon Y C，Zhao W． Stability Estimate on Meshless Unsymmetric Collocation Method for Solving Boundary Value Problems [J]．Engineering Analysis with Boundary Elements，2013，37：666-672．

[11] Hardy R L． Multiquadric Equations of Topography and Other Irregular Surfaces[J]．Journal of Geophysical Research，1971，8 (76)：1905 -1915．

[12] Pang G F，Chen W，Fu Z J． Space-fractional Advection-dispersion Equations by the Kansa Method[J]．Journal of Computational Physics，2015，293：80-96．

[13] Yao G M． Local Radial Basis Function Methods for Solving Partial Differential Equations[D]． Hattiesburg： University of Southern Mississippi，2010．

[14] Islam S U，Vertnik R，Sarlera B． Local Radial Basis Function Collocation Method Along with Explicit Time Stepping for Hyperbolic Partial Differential Equations[J]．Applied Numerical Mathematics， 2013， 67：136-151．

[15] Khoshfetrat A，Abedini M J． Numerical Modeling of Long Waves in Shallow Water Using LRBF-DQ and Hybrid DQ/LRBF-DQ[J]．Ocean Modeling，2013，65：1-10．

[16] Yao G M，Kolibal J，Chen C S． A Localized Approach for the Method of Approximate Particular Solusions[J]．Computers and Mathematics with Applications，2011，61：2376-2387．

[17] Omohundro S M． Efficient Algorithms with Neural Network Behaviour [J]．Journal of Complex Systems， 1987，1：273-347．

[18] Su J B，Zhu F，Geng Y，et al． Numerical Study of Wave Overtopping Based on Local Method of Approximate Particular Solution Method [J]．Advanced in Mechanical Engineering，2014：1-9．

[19] Lin C Y，Gu M H，Young D L，et al． Localized Method of Approximate Particular Solutions with Cole-Hopf Transformation for Multi-dimensional Burgers Equations[J]．Engineering Analysis with Boundary Elements，2014，40：78-92．

[20] Zhang X Y，Tian H Y，Chen W． Local Method of Approximate Particular Solutions for Two-dimensional Unsteady Burgers' Equations[J]．Computers and Mathematics with Applications，2014，66：2425-2432．

[21] Fu Z J，Chen W，Ling L． Method of Approximate Particular Solutions for Constant-and Variable-order Fractional Diffusion Models[J]．Engineering Analysis with Boundary Elements，2015，57：37-46．

[22] Lin J，Chen W，Sze K Y． A New Radial Basis Function for Helmholtz Problems[J]．Engineering Analysis with Boundary Elements，2012，36(12)：1923-1930．

[23] 劉從建，陳文，王海濤，等．自適應快速多極正則化無網格法求解大規模三維位勢問題[J]．應用數學和力學，2013，3(3)：259-271．

[責任編輯王康平]

收稿日期：2015-09-08

基金項目：國家杰出青年基金資助項目(11125208)；國家自然科學基金(11302069,11372097)；江蘇省自然科學基金項目(BK20150795)；中央高校基本科研業務費專項資金資助(2015B11914)

通信作者：陳文(1967-)，男，教授，博士，主要研究方向為計算力學算法和軟件．E-mail：chenwen@hhu.edu.cn

DOI：10.13393/j.cnki.issn.1672-948X.2016.01.011

中圖分類號：O413.1

文獻標識碼：A

文章編號：1672-948X(2016)01-0051-06

Numerical Simulation of 2D Schr?dinger Equation by Localized Method of Approximate Particular Solution Based on Radial Basis Functions

Hong YongxingChen WenLin Ji

(College of Mechanics & Materials, Hohai Univ., Nanjing 211100, China)

AbstractThe localized method of approximate particular solution (LMAPS) is a truly meshless collocation method with merits of high accuracy, rapid convergence and easy-to-program. In the LMAPS, the interest domain is divided into several subdomains; and in each subdomain, the localized low-rank matrix is formed based on the local nodes. And then a sparse system is constructed by reformulating the localized formulation into globalized form, which can be solved efficiently. In this paper, the LMAPS is applied to simulate 2D Schr?dinger equation where the implicit-Euler scheme is used for time discretization; and the LMAPS is utilized for space discretization. Numerical experiments reveal the high accuracy, fast convergence and lower computational time of the LMAPS. And it has potential to engineering applications.

KeywordsImplicit-Euler;2D Schr?dinger equation; Multiquadric function; localized method of approximate particular solution;global method of approximate particular solution