童 慧，孔令華，王 蘭

(江西師范大學(xué)數(shù)學(xué)與信息科學(xué)學(xué)院，江西南昌 330022)

0 引言

自由粒子的運(yùn)動(dòng)在數(shù)學(xué)上可以用無量綱化的Dirac方程來描述.本文考慮Dirac方程的第1類初邊值問題［1］:

容易驗(yàn)證以上問題滿足電荷守恒，即

上式是量子物理中的重要守恒律，數(shù)值方法對(duì)此守恒律的保持具有重要的意義和作用.對(duì)這一問題的數(shù)值方法的研究是近年來的熱點(diǎn)問題［2-9］.

保持動(dòng)力系統(tǒng)辛或多辛結(jié)構(gòu)的數(shù)值格式在長(zhǎng)時(shí)間數(shù)值模擬方面具有一定的優(yōu)勢(shì)，其重要構(gòu)造方法之一就是Runge-Kutta方法，但是這種方法得到的多辛格式往往是完全隱式的.對(duì)于非線性問題需要求解非線性代數(shù)方程組［10］，應(yīng)用起來非常消耗資源，效率比較低.而分裂步多辛算法很好地克服了這一弊端［11-15］.其基本思想是把原來的多辛哈密爾頓系統(tǒng)分裂成若干個(gè)辛或多辛哈密爾頓子系統(tǒng)，然后再對(duì)這些子系統(tǒng)構(gòu)造辛或多辛算法.這一算法已經(jīng)初步顯示了它的優(yōu)越性和高效的計(jì)算能力.為進(jìn)一步提高計(jì)算效率和計(jì)算精度，在空間方向經(jīng)常用具有高精度的能夠保持辛結(jié)構(gòu)的緊致格式離散.

多辛哈密爾頓系統(tǒng)的一般形式為

其中z∈Rm，M，K∈Rm×m是反對(duì)稱矩陣，S(z)是光滑實(shí)函數(shù)，又稱哈密爾頓能量泛函.它滿足多辛守恒律:

其中微分形式ω =1/2d z∧M d z，κ=1/2d z∧K d z.

1 Dirac方程的多辛結(jié)構(gòu)及其分裂子系統(tǒng)

Dirac方程(1)具有天然的多辛結(jié)構(gòu)，即不需要引進(jìn)正則動(dòng)量就可以得到其多辛哈密爾頓系統(tǒng)的形式.為此令 ψj=pj+i qj，j=1，2;pj，qj為實(shí)函數(shù)，z=(q1，p1，q2，p2)T.把實(shí)部和虛部分開，得到相應(yīng)的實(shí)函數(shù)方程

從而得到辛結(jié)構(gòu)矩陣

哈密爾頓函數(shù)為

相應(yīng)的多辛守恒律為

由S(z)的表達(dá)式可知，這是1個(gè)不可分的哈密爾頓系統(tǒng)，只能構(gòu)造完全隱式的多辛格式.為了有效的求解方程(1)，可以利用分裂方法，將方程(1)分裂成線性子問題和非線性問題:

(2)式中的線性子系統(tǒng)具有與原系統(tǒng)一樣的多辛結(jié)構(gòu)，非線性子系統(tǒng)的辛結(jié)構(gòu)為

易知非線性子系統(tǒng)具有點(diǎn)點(diǎn)的概率守恒:?x，t，

這一點(diǎn)點(diǎn)守恒律將使得非線性子問題可以精確求解.用分離變量法計(jì)算可得其精確解為

2 辛Euler格式及分裂格式

利用平行線 xj=-L+jh(j=0，1，2，…，J)，tn=nτ(n=0，1，2，…，N)，剖分時(shí)空區(qū)域，其中h，τ分別是x與t方向的步長(zhǎng)，用記號(hào)ψnj表示 ψ(x，t)在點(diǎn)(xj，tn)的近似.

2.1 空間的離散

在空間方向用高階緊致格式去離散.考慮1階導(dǎo)數(shù)fj'的x離散.運(yùn)用文獻(xiàn)［11］中的離散公式，有如下形式:

其中α，a1，b1都是待定參數(shù)，它們的數(shù)值由階數(shù)的大小決定.

對(duì)上面的方程中各項(xiàng)在x=xj處作泰勒展開，可以得到方程組

選取α為自由參數(shù)，可解得a1=2(α+2)/3，b1=(4α-1)/3.此時(shí)的截?cái)嗾`差主項(xiàng)為 4(3α-1)f(5)h4/5!.若取 α =1/3，則 a1=14/9，b1=1/9，截?cái)嗾`差主項(xiàng)為4f(7)h6/7!.

把fx的離散格式寫成線性算子的形式:

其中 R fj=(fj-1+3fj+fj+1)/3，I fj=(-fj-2-28fj-1+28fj+1+fj+2)/(36h).方程(3)對(duì)應(yīng)的矩陣形式為

這里

注1 A是循環(huán)對(duì)稱矩陣，B是循環(huán)反對(duì)稱矩陣，從而A－1B是循環(huán)反對(duì)稱矩陣，因而可以利用循環(huán)矩陣的有關(guān)性質(zhì)計(jì)算它的特征值.

2.2 時(shí)間的離散

時(shí)間方向用辛Euler格式離散，即

或

則非線性Dirac方程(2)的離散格式為

格式(4)中的前兩式是對(duì)非線性子問題的精確求解，從而滿足離散的辛守恒，而后兩式是多辛格式.因而整個(gè)格式是辛格式，而且具有多辛格式的特征.

2.3 格式穩(wěn)定性分析

這里討論Dirac方程的緊致分裂多辛格式(4)的穩(wěn)定性.非線性子問題是精確求解，而且解是穩(wěn)定，因而對(duì)格式穩(wěn)定性的研究只需要關(guān)注線性子問題的穩(wěn)定性.為此令r=τ/h，D=hA-1B，則離散格式(4)的增長(zhǎng)矩陣為

當(dāng)增長(zhǎng)矩陣的譜半徑小于等于1時(shí)，格式穩(wěn)定.考慮矩陣A，B的階為N=2k.因?yàn)锳，B是循環(huán)矩陣，利用循環(huán)矩陣的相關(guān)性質(zhì)，可得它們的特征值分別為

而λN-j-2與λj共軛，這里的特征值都為實(shí)數(shù)，

而λN-j-2與λj共軛，且它們都是純虛數(shù).又因?yàn)锳，B都是循環(huán)矩陣，則存在同一個(gè)正交矩陣Q，使得

即為h A-1B=D的特征值.按照上述對(duì)角線上特征值的排列順序，將G的特征值依次記為 λ1(D)，λ2(D)，…，λN(D)，則矩陣G的特征方程為

當(dāng)矩陣A，B的階為奇數(shù)時(shí)，λ1(A)=5/3，λj(A)=1+(2/3)cos［2(j-1)π/N］，j=2，…，k+1.而 λN+2-j與 λj共軛..而 λN+2-j與 λj共軛.接下的討論與上述類似.

3 數(shù)值實(shí)驗(yàn)

本部分考察Dirac方程的緊致分裂多辛格式的數(shù)值行為

其中

取Λ =0.75，根據(jù)初始條件，方程的孤立波解為

為了驗(yàn)證該數(shù)值方法的理論精度為O(τ+h6)，并且驗(yàn)證此格式比傳統(tǒng)多辛格式在計(jì)算效率方面更具有優(yōu)勢(shì).考慮x∈［-50，50］，取不同的h和τ計(jì)算，計(jì)算的最終時(shí)刻是T=1.表1和表2記錄了數(shù)值模擬的結(jié)果，其中格式(a)為本文緊致分裂多辛格式，格式(b)為文獻(xiàn)［9］中Lie-Trotter分裂的格式，格式(c)為文獻(xiàn)［4］中2階隱多辛龍格庫(kù)塔格式.

表1 當(dāng)τ=1×10-6時(shí)ψ1的誤差、收斂階、各個(gè)計(jì)算消耗的時(shí)間

表2 當(dāng)h=0.1時(shí)ψ1的誤差、收斂階、各個(gè)計(jì)算消耗的時(shí)間

表1列出了當(dāng)時(shí)間步長(zhǎng)τ=1×10-6，空間網(wǎng)格點(diǎn)數(shù)N=100，200，400時(shí)在T=1時(shí)刻的數(shù)值解的誤差和階數(shù)，由此可驗(yàn)證該數(shù)值解在空間上以O(shè)(h6)逼近精確解，同時(shí)與文獻(xiàn)［4，9］的格式對(duì)比，可以看出緊致分裂多辛格式(a)在空間精度上比格式(b)和格式(c)更高，運(yùn)行時(shí)間也小于另外2種格式.

表2列出了當(dāng)時(shí)刻T=1，空間步長(zhǎng)h=0.1時(shí)采用不同的時(shí)間步長(zhǎng)的數(shù)值解的誤差和階數(shù)，由此可以看出該數(shù)值解在時(shí)間上以O(shè)(τ)逼近精確解，與理論結(jié)果相符，并與另外2種算法進(jìn)行比較，緊致分裂多辛格式的計(jì)算時(shí)間要少于另外2種格式的計(jì)算時(shí)間.

［1］ Alvarez A，Carreras B.Interaction dynamics for the solitary waves of a nonlinear Dirac model［J］.Phys Lett A，1981，86(6/7):327-332.

［2］ Alvarez A.Linear Crank-Nicholson scheme for nonlinear Dirac equations［J］.JComput Phys，1992，99(3):348-350.

［3］Alvarez A，Kuo P，L Vazquez.The numerical study of a nonlinear one-dimensional Dirac equation ［J］.Appl Math Comput，1983，13(1):1-15.

［4］Hong Jialin，LiChun.Multi-symplectic Runge-Kuttamethods for nonlinear Dirac equations ［J］.J Comput Phys，2006，211(2):448-472.

［5］Wang Han，Tang Huazhong.An efficient adaptivemesh redistribution method for a nonlinear Dirac equation［J］.J Comput Phys，2007，222(1):176-193.

［6］王健.Dirac方程的多辛格式［J］.上海交通大學(xué)學(xué)報(bào):自然科學(xué)版，2004，38(5):849-852.

［7］周祖珍，李淮江.Dirac方程的對(duì)稱性與守恒定律［J］.云南師范大學(xué)學(xué)報(bào):自然科學(xué)版，1992，12(3):34-36.

［8］邵嗣烘，湯華中.非線性Dirac方程的數(shù)值研究［J］.高等學(xué)校計(jì)算數(shù)學(xué)學(xué)報(bào)，2005，27(S1):123-126.

［9］丁建，徐君祥，張福保.非周期 Dirac方程的穩(wěn)態(tài)解［J］.中國(guó)科學(xué):數(shù)學(xué)，2011，41(6):517-534.

［10］Wang Yushun，Hong Jialin.Multi-symplectic algorithms for Hamiltonian partial differential equations［J］.Commun Appl Math Comput，2013，27(2):163-230.

［11］Hong Jialin，Kong Linghua.Novelmultisymplectic integrators for nonlinear fourth-order Schr?dinger equation with trapped term ［J］.Commun Comput Phys，2010，7(6):613-630.

［12］Kong Linghua，Cao Ying，Wang Lan.Split-step multisymplectic integrator for the fourth-order Schr?dinger equation with cubic nonlinear term ［J］.Chin J Comput Phys，2011，28(1):76-82.

［13］王蘭，符芳芳，童慧.Dirac方程的分裂步多辛格式［J］.江西師范大學(xué)學(xué)報(bào):自然科學(xué)版，2013，37(5):462-465.

［14］王蘭.多辛Preissmann格式及其應(yīng)用［J］.江西師范大學(xué)學(xué)報(bào):自然科學(xué)版，2009，33(1):42-46.

［15］黃紅，王蘭.薛定諤方程的局部1維多辛算法［J］.江西師范大學(xué)學(xué)報(bào):自然科學(xué)版，2011，35(5):455-458.