上海交通大學碩士學位論文

目錄

第一章引言

-2-

1.1耦合長短波方程的背景

-2-

1.1.1Schrodinger方程的介紹

-2-

1.1.2Kdv方程的介紹

-8-

1.1.3Schrodinger-KdV方程組及耦合長短波方程的由來

-10-

1.2Hamilton系統(tǒng)、辛算法及多辛算法

-11-

1.3原有的數(shù)值方法

-17-

1.3.1時間分裂方法

-17-

1.3.2Crank-Nicolson方法

-19-

第二章多辛格式

-22-

2.1EulerBox格式

-24-

2.2Preissman格式

-30-

2.3Fourier擬譜格式

-35-

第三章數(shù)值實驗

-40-

3.1E_LS、CNI與Box_ls的誤差及階數(shù)

-40-

3.2TSS與LS1的誤差及階數(shù)

-49-

第四章結論

-55-

參考文獻

-56-

致謝

-59-

第一章引言

1.1耦合長短波方程的背景

1.1.1Schrodinger方程的介紹

薛定諤方程(Schrodingerequation)是由奧地利物理學家薛定諤提出的量子力學中的一個基本方程，也是量子力學的一個基本假定，其正確性只能靠實驗來檢驗。是將物質波的概念和波動方程相結合建立的二階偏微分方程，可描述微觀粒子的運動，每個微觀系統(tǒng)都有一個相應的薛定諤方程式，通過解方程可得到波函數(shù)的具體形式以及對應的能量，從而了解微觀系統(tǒng)的性質。

一．關于波的簡單介紹

波的形式是多種多樣的，平面波是指在波的傳遞過程中在一段距離內(nèi)波幅變動不大的波。描述平面波的狀態(tài)的波函數(shù)為：

(1.1)

為波函數(shù)。x與t分別表示x方向的距離變量與時間變量。函數(shù)計算結果的物理量為波幅，即長度單位為米。

A為波幅，單位為米。指平面波相對穩(wěn)定的最大波幅絕對值。

cos(kx--wt)為波幅的變動系數(shù)，數(shù)值在-1，0，+1之間變動。

k為單位長度弧度數(shù)，即，單位為弧度／米。

為單位時間弧度數(shù)，即，單位為弧度／秒。

由于=波速度，=即波在傳播過程中不同距離的時間差。波函數(shù)(1.1)式可以寫成：

(1.2)

如果上述波函數(shù)(1.2)式去掉括號中的第一項，得到：

(1.3)

即為一個以時間為變量的振幅函數(shù)，這樣的波并不向外傳播，而是在原地上下振動，余弦cos的正負角的函數(shù)值是相等的，因此括號中的負號可以忽略。

如果波函數(shù)(1.2)式中去掉括號中的第二項同時保留第一項，得到：

(1.4)

即為一個以坐標為變量的函數(shù)，是描述波到達各種位置的振動位移情況。

波函數(shù)(1.1)式也可以寫成以e為底的復數(shù)形式：

(1.5)

通用的波動方程是如下形式：

(1.6)

化簡為一維形式的波動方程：

(1.7)

函數(shù)是直接描述變量之間的關系，而物理方程則是在復雜的情況下表述各種物理量的關系，它可以包含各種函數(shù)，但對函數(shù)中的變量之間關系的表述不是那么直接，函數(shù)表達式與方程之間，在一定條件下可以轉換。下面給出具體的推導過程。

對宏觀平面波函數(shù)(1.5)式進行轉換，得：

(1.8)

對時間二次微商得：

整理得：

對坐標二次微商得：

整理得：

將兩個二次微商的結果連接得：

引入波速u，將上式整理得到：

這就是前面的一維形式的波動方程(1.7)式?？梢姴ê瘮?shù)與波動方程在一定條件下可以相互轉換。應當注意到由波函數(shù)轉換為波動方程的過程中，波幅常量A消失了，說明這個波動方程中含有的波函數(shù)與波幅常量沒有直接聯(lián)系。如果由波動方程反解出波函數(shù)，這個常量A需經(jīng)過歸一化處理重新找回來。

通常，參照一般的波動描述方式，似乎可以建立一個類似的電磁波的波函數(shù)表達式：

然而，我們以宏觀的平面波方式來描述量子波時，必須考慮到量子波與平面波的重要區(qū)別：平面波的波函數(shù)描述的是一個波群(波束)，而量子波由于其量子性，表現(xiàn)為個體的量子波。即使一個有限的空間在一段有限的時間里，只有一個電子或者一個光子在運動，我們也得承認這個電子或者光子的波動性，應當也有波函數(shù)描述它。因此用波函數(shù)來描述量子波時，顯然應當與描述光束或者電子束的方式有所區(qū)別，或者說不能用描述一般平面波的波函數(shù)的方法來描述量子波。

分析光子的運動：光子在振動且以光速運動，以位置作為變量，光子的相位隨之而變，換以時間作為變量，光子的相位也是隨之而變。光子的波幅、一個波動周期內(nèi)的速度也是在隨位置和時間變動的，但是波動量的大小與相位的變動是相關的。因此首先選取相位作為一個波動周期的因變量。

將原來描述光束的函數(shù)式(1.8)改為描述單個光子的波函數(shù)：

(1.9)

函數(shù)的值就是相位，單位為弧度。

將上述描述光子相位的函數(shù)用三角函數(shù)表示：

(1.10)

得到的是無量綱數(shù)，數(shù)值在-1，0，1之間變動，這就是我們描述概率幅的函數(shù)。加上適當?shù)某?shù)項A，使得波函數(shù)有適當?shù)奈锢砹?，將模平方并按一定的物理量積分后，得到的結果是有量綱數(shù)，即概率密度：

(1.11)

根據(jù)(1.8)式對電磁波的函數(shù)表達式，按照式(1.6)的通用波動方程，可以建立一個電磁波的波動方程：

(1.12)

方程中E表示什么呢?根據(jù)前面對一般波函數(shù)與波動方程的轉換分析我們知道，這個E無法代表電場強度，因為由波函數(shù)轉換為波動方程的過程中，電場強度常量E已經(jīng)消失了，現(xiàn)在這個E函數(shù)計算得到的值是無量綱數(shù)，在一定條件下按歸一化處理后可以有一定的物理量綱，就是概率幅。由此可見，即使人們將波動視為量子的出現(xiàn)概率的似波性，在進行歸一化處理之前，波動方程并不能完整地描述量子的波動狀態(tài)，因為波動方程計算出來的概率幅數(shù)值在-1,0,1之間變動，必須考慮行程、波數(shù)等因素，進行適當?shù)霓D換后才能得到實際波動的概率幅數(shù)值。因此歸一化處理，實際上是使波動方程描述的物理現(xiàn)象回歸實在的波動。筆者將在后面的實例中進一步說明，波動方程加上波動的邊界條件，再加上歸一化處理，是可以描述量子實在的。

二．薛定諤方程

建立描述電子運動的波函數(shù)

既然電子的運動可以形成波，那么應當如同描述光子的波動一樣，有一個波動方程描述電子的波動。參照前面描述光子的波動方程(1.10)式建立描述電子的波函數(shù)：

上式是量子單位制形式。將上面的波函數(shù)按下式描述為經(jīng)典形式：

由于動量，即動能／速度。能量E=，此處表示電子庫侖勢能。普朗克常數(shù)h=1，

(1.13)

將電子的波函數(shù)改寫為復數(shù)形式并進行微商

將(1.13)式改寫為復數(shù)形式：

(1.14)

這里為了波函數(shù)的完整性加上了A，在后面進行微商處理時A自然消失，而在歸一化處理時A又回到公式中，那時它的物理意義與我們進行歸一化處理方式有關。

將上述函數(shù)對時間微商得：

(1.15)

將上述函數(shù)對坐標微商得：

(1.16)

為了將能量E引入，將(1.16)式再對x求一次偏導，得到：

(1.17)

轉換中用上的公式，是由德布羅意波的動量與能量關系決定的。

得到薛定諤方程

由(1.15)式得到：

(1.18)

由(1.17)式得到：

(1.19)

將(1.18)式與(1.19)式連接起來，就得到自由粒子的一維含時不計算勢能的薛定諤方程：

(1.20)

如果考慮到量子的庫侖勢能，得到一維含時的薛定諤方程：

(1.21)

將上式擴展到三維，得到三維含時的薛定諤方程：

(1.22)

如果不考慮時間因素，得到三維定態(tài)的薛定諤方程：

(1.23)

其中：

簡化為一維定態(tài)的薛定諤方程：

(1.24)

定態(tài)是指粒子的概率密度分布不隨時間變化的狀態(tài)，定態(tài)薛定諤方程示，任何時候粒子概率幅的空間分布，不隨時間變化。

1.1.2Kdv方程的介紹

1834年英國科學家ScottRussell偶然觀察到了一種奇妙的水波。1844年，他在《英國科學促進協(xié)會第14屆會議報告》上發(fā)表的《論波動》一文中，對此現(xiàn)象作了生動的描述：“我觀察過一次船的運動，這條船被兩匹馬拉著沿狹窄的運動迅速前進著，突然，船停了下來，而被船所推動的大堆水卻并不停止，它們積聚在船頭周圍激烈地擾動著，然后水浪突然呈現(xiàn)在一個滾圓而平滑、輪廓分明的巨大孤立波峰，它以巨大的速度向前滾動著，急速地離開了船頭，在行進中它的形狀和速度并沒有明顯的改變，我騎在馬上緊跟著觀察，它以每小時約八、九英里的速度滾滾向前，并保持著長約30英尺，高約1-1.5英尺的原始形狀，漸漸地它的高度下降了。當我跟蹤1-2英里后，它終于消失在逶迤的河道之中”。這就是Russell觀察到的奇特現(xiàn)象，進而他認為這種孤立的波動是流體運動的一個穩(wěn)定解，并成它為“孤立波”。Russell當時未能成功地證明并使物理學家們信服他的論斷，從而埋汰數(shù)學家未能從已知的流體運動方程預言出這一現(xiàn)象，之后有關孤立波的問題在當時許多物理學家中引起了廣泛的爭論。直到60年后的1895年，Kerteweg，DeVries研究了淺水波的運動，在長波近似和小振幅的假定下，建立了單向運動的淺水波運動方程

下面就簡單推導一個Kdv方程：

在豎直平面考察平面內(nèi)流體的運動，考慮下面定解問題

其中是水的深度，h(t,x)是水波的函數(shù)。

引進小參數(shù)a是平面波振幅，是波長。作下列變換：

則上述方程組變?yōu)椋?/p>

把去掉還原成：

用構造

所以成立

構造的滿足因為

將代入到得

將代入到得

下面就建立關于h的偏微分方程：

在上述兩式中忽略一次以及以上的項得：

取可得，展開

在上述兩式中忽略二次以及以上的項得

將代入到剛得到的兩項中分別可得：

取可得

這就是經(jīng)典的Kdv方程。

1.1.3Schrodinger-KdV方程組及耦合長短波方程的由來

各種客觀環(huán)境已經(jīng)在研究長短波之間的相互作用現(xiàn)象，尤其是在流體力學、等離子物理和化學物理。短波通常用Schrodinger方程來描述，長波一般用帶有色散的某種類型的波方程來描述。本文中，我們用多辛積分[1,2,3]格式來研究耦合長短波方程(LSIE)。

(1.25)

滿足

(1.26)

復值函數(shù)表示短波的包絡層，實值函數(shù)表示長波的振幅。是實常數(shù)。在方程(1.25)被用來在重力和微小范圍模式下建立表面波。這個方程是Schrodinger-KdV方程組的一個特殊例子。

(1.27)

關于初始邊界問題(1.25)--(1.26)，可以得到下面命題

命題1.1帶有初始邊界問題(1.25)--(1.26)的解滿足下面的守恒律

是復函數(shù)φ的共軛函數(shù)。

方程組(1.25)具有獨特的數(shù)學特點因為研究發(fā)現(xiàn)它具有完整的可積結構。特別的，它有一個逆散射變換和顯示的N孤子解[5,6]。各種數(shù)值技術尤其是有限元方法、譜方法和時間分裂步方法都已經(jīng)用來解各種形式的長短波相互作用的方程。

1.2Hamilton系統(tǒng)、辛算法及多辛算法

力學研究中一個非常重要的是對稱的幾何觀點，它不論是從基本的原理公式出發(fā)，還是到具體的應用，都特別強調了幾何方法和力學研究的密不可分的關系。其中Hamilton提出的力學定理可以用微分流形來研究剛體體系及太陽系等復雜系統(tǒng)的力學性質；可以用相應的Hamilton函數(shù)的對稱性來研究能量、線性動量與角動量等Hamilton系統(tǒng)的守恒性質。

Hamilton系統(tǒng)是一種重要的力學系統(tǒng)，廣泛的出現(xiàn)在物理、力學、工程、純數(shù)學與應用數(shù)學領域。通?？梢哉J為，一切耗散可忽略不計的真實物理過程都可以表示成Hamilton方程的形式，而它們的共同基礎都是辛幾何。辛幾何歷史可以追溯到十九世紀英國物理學家和數(shù)學家Hamilton，他為了研究Newton力學，引入廣義坐標和廣義動量來表示系統(tǒng)的能量，即Hamilton函數(shù)。

考慮Newton運動方程，設表示質點位置，它滿足拉格朗日方程其中表示動能，表示勢能。引入共軛動量，其滿足

下面定義則

令，那么上述方程可表示如下：

其中，為n階單位矩陣，稱為有限維Hamilton系統(tǒng)，H(z)為Hamilton函數(shù)。

Hamilton系統(tǒng)有兩個重要特性：守恒性與辛結構。我們知道Hamilton系統(tǒng)的解是一個單參數(shù)的保測變換,即辛變換。因此在研究Hamilton系統(tǒng)的計算方法時候，都要使得離散后的方程保持原有系統(tǒng)的辛結構。用傳統(tǒng)的數(shù)值算法模擬Hamilton系統(tǒng)，會破壞系統(tǒng)的辛結構，進而使數(shù)值模擬不能保持長時間的穩(wěn)定。因此構造一種能保持Hamilton系統(tǒng)的辛結構的算法，具有重要意義，稱能保持Hamilton系統(tǒng)辛結構的算法為辛算法。

無限維Hamilton系統(tǒng)一般情況下可以表示為：

(1.28)

其中，z=z(x,t)為狀態(tài)向量，H為Hamilton函數(shù)：

(1.29)

為H的Frechet導數(shù)。

等式(1.28)也可通過Lagrange泛函導出?？紤]如下非線性Klein-Gordon方程，其中，為某一光滑非線性函數(shù)。

對以下的Lagrange泛函作變分

(1.30)

其中Lagrange密度為：

就可以得到Euler-Lagrange方程：

對Lagrange密度L作Legendre變換，即令得到一階方程組

(1.31)

令

則(1.31)等價于Hamilton系統(tǒng)(1.28)。相應的辛二形式為：

從以上過程可知，Legendre變換僅對時間方向進行。若對空間方向也進行Legendre變換，即可得到相應一階方程

(1.32)

Hamilton函數(shù)為：

相應空間方向的辛二形式為：

大量的孤立波方程，可表示為無限維Hamilton系統(tǒng)。例如：

非線性Schrodinger方程：

其Hamilton函數(shù)為：

其中

KdV方程：

其Hamilton形式為：

將求解有限維Hamilton系統(tǒng)的辛算法推廣到無限維最直接有效的方法是：先對空間方向進行離散，這樣離散以后的系統(tǒng)為有限維Hamilton系統(tǒng)，然后應用辛算法。常用的半離散方法為差分法[20,21,22,23,24]，擬譜方法和配置點法[25]等。

在上例導出非線性Klein-Gordon方程的無限維Hamilton形式時，采用的是不完全Legendre變換。若對方程(1.29)采用完全的Legendre變換，即同時令，可得到相應一階方程組：

(1.33)

記(1.33)等價于：

(1.34)

其中

Bridges和Reich稱形如(1.34)的方程為具有多辛結構的Hamilton系統(tǒng)。

求解有限維Hamilton系統(tǒng)的辛算法與其它方法相比具有許多優(yōu)點，其中有一點就是辛算法可以對守恒量長時間數(shù)值模擬。但是當用辛算法對無限維Hamilton進行離散時具有局限性，具體表現(xiàn)在守恒量是全局性的。為了克服此局限性，Marsden等[26,27,28,29]和Bridges[30,31,32]分別從Lagrange系統(tǒng)和Hamilton系統(tǒng)出發(fā)，提出了多辛Hamilton系統(tǒng)和多辛算法的概念。Marsden等是從變分原理開始，由邊界項得到辛形式，限制到具體方程上得到該方程的多辛守恒律；Bridges是將有限維的Hamilton系統(tǒng)推廣到無限維Hamilton系統(tǒng)，使得偏微分方程在時間方向和空間方向上都有各自不同的辛結構。

多辛結構的良好特征是它可以導出一個包含微分二形式的守恒律，即多辛守恒律。與多辛守恒律相對應的是多辛Hamilton系統(tǒng)具有能量守恒律和動量守恒律。由于這些守恒律不依賴邊界，因而都是局部守恒律。

Bridges和Reich考慮了一般情形的多辛Hamilton系統(tǒng)。設是有限維相空間，S為Q上的函數(shù)，M，K為兩個反對稱矩陣，則稱

(1.35)

為多辛Hamilton系統(tǒng)。其中，為反對稱矩陣；為狀態(tài)向量及為光滑函數(shù)；表示函數(shù)關于z的梯度，為底空間的兩個獨立變量。

具有m個空間變量的一般多辛Hamilton系統(tǒng)可表示為：

(1.36)

洪佳林等對(1.36)作了詳細的討論[33,34,35]，為簡單起見，我們僅討論具有一個空間變量，且系數(shù)矩陣為常系數(shù)的情形，即(1.35)中。

Bridges和Reich證明了多辛Hamilton系統(tǒng)(1.35)具有如下局部守恒律，即多辛守恒律、局部能量守恒律和局部動量守恒律。

Hamilton系統(tǒng)(1.35)滿足多辛守恒律：

(1.37)

其中分別對應于t和x方向的兩個不同的辛結構。

多辛Hamilton系統(tǒng)(1.35)具有局部能量守恒律：

(1.38)

和局部動量守恒律：

(1.39)

利用矩陣分解：

Moore[36]得到了多辛守恒律、局部能量和動量守恒律的另一種表示，即

與無限維Hamilton系統(tǒng)比較，多辛Hamilton系統(tǒng)具有的守恒律都是局部的，但是它們在任意的時空區(qū)域內(nèi)成立，而且不依賴于邊界條件。若對局部守恒律在空間方向積分，并利用適當?shù)倪吔鐥l件，就可導出相應的整體守恒律即辛守恒律、能量守恒律和動量守恒律。而這些適當?shù)倪吔鐥l件顯然也是整體守恒律的必要條件。因此，局部守恒律蘊涵著整體守恒律。另外多辛Hamilton系統(tǒng)的相空間是有限維的，因而可應用有限維的理論來分析方程的性質。

1.3原有的數(shù)值方法

1.3.1時間分裂方法

時間分裂方法[31-35]是近些年研究非線性偏微分方程數(shù)值解時使用較多的方法，目前常見的有分裂差分方法[36-40]、分裂譜方法[41-46]等等，時間分裂方法的主要優(yōu)點是把方程的非線性項與線性項分開計算，把較難處理的非線性項單獨處理。

時間分裂方法的基本原理[37]如下：

對方程：

(1.40)

方程(1.40)改寫為如下形式：

(1.41)

其中分別為方程的線性項和非線性項。

在時間步上，由(1.41)有：

(1.42)

即有：

(1.43)

設是的近似，由(1.42)可得

(1.44)

則可以得到如下的分裂格式：

(1.45)

(1.46)

其中為函數(shù)在時間層對u的逼近。

因此方程(1.40)可以采用下面的時間分裂步來進行計算：

分裂步1：(1.47)

分裂步2：(1.48)

方程(1.47)是一個非線性的常微分方程，可以精確的求解，方程(1.48)可以采用不同的數(shù)值方法來計算，常見的有分裂譜方法，分裂差分方法，如分裂C-N格式等等。

分裂格式(1.47)-(1.48)在時間方向上為一階的。根據(jù)Strang[31]的分裂思想，可構造如下在時間方向上二階的分裂格式：

(1.49)

(1.50)

(1.51)

本格式相當于對(1.47)分別在時間步和進行了兩次的計算。

對本問題利用時間分裂方法，其步驟如下：

Step1.方程組(1.25)的第一個方程可分解為：

(1.52)

Step2.對方程組(1.52)的第一個方程做傅里葉變換得：

得

進一步解得

Step3.對方程組(1.25)的第二個方程做傅里葉變換得：

然后對上述方程做差分得進而可求出：

Step4.直接解方程組(1.52)的第二個方程可得：

Step5.最后求出：

這一小節(jié)先是對時間分裂步方法的基本原理的一個簡單介紹，然后簡單推導得出時間分裂步方法在本篇文章當中對耦合長短波方程的實現(xiàn)步驟，當然以上五個步驟也是在用matlab進行數(shù)值試驗時的基本編程思路。在本文的數(shù)值試驗部分，取特定的參數(shù)及變量值就可以按照上述步驟編寫程序得出相應的結論。

1.3.2Crank-Nicolson方法

對方程組(1.25)進行以下的差分離散：

記如下：

其中下標N是時間節(jié)點數(shù)。

將上述方程組的第一個等式按照定義的合并空間指標j之后得到下面這個等式：

其中A是單位矩陣，B和運算是如下格式：

上式兩邊同時乘以：

化簡得：

令上式可變?yōu)椋?/p>

同理方程組第二式可化為：

其中

綜上所述最后推出迭代格式為：

本文研究的方程是耦合的長短波方程，它是Schrodinger-KdV方程組的一個特例。本文所做的工作主要是建立耦合的長短波方程的幾種多辛格式，并通過數(shù)值試驗與已有的數(shù)值算法比較。

在本文的第一章引言中，會介紹一下耦合的長短波方程的背景，然后敘述Hamilton系統(tǒng)、辛算法和多辛算法，進而會介紹兩種現(xiàn)在認可的數(shù)值算法：時間分裂方法和Crank-Nicolson方法，并對相應的格式進行簡單的推導。

在本文的第二章中，會詳細敘述新建立的三種多辛格式，它們是EulerBox格式、多辛Preissman格式和Fourier擬譜格式，給出相應的性質并格式的誤差的階數(shù)做了詳細的說明。三種多辛數(shù)值算法用來解決周期性帶有初始問題的長短波相互作用方程。

第三章是本文的數(shù)值試驗，針對上述5種數(shù)值格式，取三組時間步長dt和相應的節(jié)點數(shù)N算出各自的誤差，并根據(jù)公式求出相應的誤差的階數(shù)。

第四章是本文的結論說明，通過理論分析和數(shù)值實驗來分析新構造的三種格式，給出本篇文章的結論。

第二章多辛格式

根據(jù)多辛的定義[9]，如果偏微分方程可表示為：

(2.1)

則稱方程（2.1）為多辛Hamilton系統(tǒng)。其中：M,K是反對稱矩陣；z(x，t)是狀態(tài)變量的函數(shù)。S：是光滑函數(shù)；是Hamiltonian函數(shù)S=S(z)的梯度。

多辛格式滿足多辛守恒律

(2.2)

是外微分格式

(2.3)

定義了一個對稱的時間空間結構。多辛結構由Nothers定理[8]自然地可以得到局部守恒律。事實上，對于Hamiltonian系統(tǒng)中的偏微分方程(2.1)，當S與x，t無關，可以得到局部能量守恒律和局部動量守恒律。

(2.4)

(2.5)

對于周期邊界或者在邊界條件為0的F(z)和G(z)，可以推出整體能量和動量守恒律

(2.6)

其中，

能保持離散多辛守恒律的數(shù)值方法為多辛算法[9,10,11]。為了得到長短波相互作用方程的多辛格式，令方程組(1.25)可以用下列一階方程組表示：

(2.7)

反對稱矩陣M,K表示如下。

方程(2.1)等號右邊的S(z)為：

直接計算就可以得出方程組(3.1)滿足多辛守恒律：

(2.8)

(2.4)和(2.5)中定義的密度函數(shù)可以下面給出

(2.9)

加上適當?shù)倪吔鐥l件可以得到整體守恒律。例如在周期或者消失在無限邊界條件的情況下可對E和I在空間區(qū)域上積分，就可以得到整體能量守恒律(1.6)和動量守恒律(1.7)。

(2.10)

對于長短波相互作用方程這是兩個很重要的整體守恒律。

2.1EulerBox格式

根據(jù)[13]通過引入兩個矩陣M和K的分裂矩陣可以獲得方程組(2.1)的Euler-box格式。例如，分裂矩陣M,K如下：

(2.11)

相應的格式變?yōu)椋?/p>

(2.12)

格式(2.12)滿足離散多辛守恒律：

(2.13)

其中很明顯矩陣的分裂是不唯一的。但是在本文中，分別是M,K的上三角矩陣。這種特殊的選擇使得Euler-box格式變?yōu)橐韵滦问剑?/p>

(2.14)

離散方程組(2.14)也具有離散多辛守恒律(2.13)的形式。其中只需取為消去所有引入的變量p，q，r和s，與長短波相互作用方程組(1.25)相比可以得到下面的多辛積分：

(2.15)

其中，注意這個格式是線性隱式的，不需要每一步都解線性代數(shù)方程?？梢宰C明此格式的局部截斷誤差是。

命題2.1格式(2.15)的局部階段誤差的階是。

證明：(1)對于等式可得

截斷誤差為：

在進行泰勒展開得：

由(2.7)中的第一個等式可得

所以

進一步得知截斷誤差的階是

(2)對于等式可得

截斷誤差T在進行泰勒展開得：

由(2.7)中的第二個等式可得：

進一步得知的截斷誤差階是

(3)對于等式可得

在進行泰勒展開得：

將代入得

將等式代入上式中得：

由(2.7)中的第四、五、六式可得：

由于

所以得：

進一步得知截斷誤差的階是。證畢

命題2.2多辛Euler-box格式(2.15)滿足下面的離散守恒律：

證明：方程組(2.15)的第三個方程是利用的定義展開上式得，

兩邊同時乘以，

得對空間方向j求和，得到

由周期邊界條件可知

進一步可知

再由可以得到所需要的結論。證畢。

2.2Preissman格式

多辛Preissman格式有具有許多非常好的性質。Bridges證明了多辛Preissman格式應用于線性Hamilton系統(tǒng)時，滿足離散局部能量和局部動量守恒律;對于非線性Hamilton系統(tǒng)，Moore和Reich證明了其僅滿足半離散局部能量和動量守恒律。

由多種方法可導出多辛Preissman格式，其中最典型的方法為時間和空間方向分別利用隱式中點辛格式。多辛Preissman格式可表示為

其中

.

和分別表示空間步長和時間步長。Bridges和Reich證明了Preissman格式能精確滿足多辛守恒律。

Preissman格式滿足多辛守恒律

其中

對于線性Hamilton系統(tǒng)，多辛Preissman能精確保持局部能量和動量守恒律，如果Hamilton函數(shù)，其中A為對稱矩陣，則Preissman格式全離散局部能量守恒律

和全離散局部動量守恒律

其中:

然而，一般而言，對非線性Hamilton系統(tǒng)，離散局部能量和動量守恒律不能精確滿足。Moore和Reich提出了一種半離散能量和動量守恒律。

半離散格式

滿足半離散能量守恒律

其中

而半離散格式

滿足半離散動量守恒律

其中

由此可見，對于全離散多辛格式，它能精確保持Hamilton系統(tǒng)的離散多辛結構，但并不意味著系統(tǒng)的其它守恒量，如局部能量和動量守恒律，整體能量和動量，以及決定相空間結構的其它整體不變量守恒。而大量的數(shù)值計算顯示，多辛格式能使局部守恒律在長時間計算中保持很好。

針對多辛Preissman格式，給出能量密度E，能量流F，動量密度I和動量流G相應的離散形式。分別如下

在空間和時間中都用辛隱式中點離散可以得到Preissman格式：

(2.16)

格式(2.16)滿足下面的離散多辛守恒律：

(2.17)

其中

把多辛Preissman格式(2.16)應用到方程(1.27)的多辛形式中，得到

(2.18)

(2.18)是一個二階的多辛格式。所以它滿足相應的多辛離散守恒律(2.17)，其中只需?。?/p>

(2.19)

方程(2.18)的第一式乘以，并將(2.18)的第三式代入其中，同時(2.18)的第二式乘以，并將(2.18)的第四式代入其中得：

方程(2.18)的第五、第六和第七式消去s，r得

所以消除所有的引入的變量u，v，p，q，r和s，只用變量φ和ψ來表示多辛格式。于是我們可以得到一個新的隱式多辛格式：

(2.20)

下面由(2.20)通過消去D和M得到所需計算的迭代格式：

(1)由定義可知，所以

進一步可得：

(2)由定義可知，所以

進一步可得：

(3)由可得

假設為周期邊界，下面引入矩陣A，B以及向量：

綜合上述(1)(2)(3)的推導，(2.20)的第一式可化為：

.

.

由于(2.20)中的第一式是非線性方程，只能用迭代的來算，兩邊同時乘以得：

令得：

同理由(2.20)的第二式可得：

即(2.20)轉化為：

注意，方程(2.20)是隱式的,想要得到解需要解非線性方程組。因此,這個簡單的迭代算法使用。通過泰勒展式可知多辛格式(2.20)的截斷誤差階是。

命題2.3多辛的Preissman格式(2.18)滿足下面的守恒性質：

(2.21)

(2.22)

證明：等式(2.21)是很顯然的。很容易證明等式(2.23)是正確的。

(2.23)

用方程組(2.18)的第一個第二個等式分別乘以和，然后合并成為下面的等式：

(2.24)

由等式(2.23)和方程組(2.18)的第三個和第四個方程可以得到：

把上式代入等式(2.24)，對j，k求和，可以得到：

所以，等式(2.22)也成立。證畢。

2.3Fourier擬譜格式

譜方法是求解微分方程的高精度算法，是近二十年來發(fā)展最快的數(shù)值方法之一，它已成功地應用于許多領域的數(shù)值計算，例如流體力學、理論物理、量子力學、非線性光學等。由于用無限可微的正交函數(shù)作為基函數(shù)，譜方法具有通常有限元和差分方法無與倫比的高精度，即如果精確解是解析的，則譜方法的精度是指數(shù)階的。

對于滿足周期邊界條件的多辛Hamilton系統(tǒng)，Bridges和Reieh[43]提出了基于Fourier變換的多辛Fourier離散。多辛Fourier變換導出了在Fourier空間中的多辛概念和半離散系統(tǒng)。Islas和Schober[65,66]對多辛Fourier離散進行進一步研究，證明了對于線性Hamilton系統(tǒng)，能保證局部能量和動量守恒。陳景波和秦孟兆[44]討論了在物理空間上進行離散的多辛Fourier擬譜方法，并應用于求解非線性Schrodinger方程。

Fourier擬譜方法分二個基本步驟。第一步為通過解在配置點上三角多項式插值構造解的離散形式；第二步利用解在配置點上的值求出導數(shù)值，即構造微分矩陣。為簡單起見，設所討論的區(qū)間。對任何整數(shù)N>0，令為空間步長，

為配置點；

為插值空間，其中為

插值算子助定義如下：對任意

利用

及正交性

得：

對定義離散內(nèi)積和離散模

在配置點上的導數(shù)值可由函數(shù)值和微分矩陣求得

其中表示一階Fourier微分矩陣，其元素為

顯然，為反對稱矩陣。

傅里葉變換可以保持偏微分方程的多辛性質不變。離散傅立葉方程恢復了標準譜離散的性質，這導致哈密頓常微分方程系統(tǒng)標準辛[11,14]可積。

對方程(1.27)用傅立葉擬譜方法，是一階的傅立葉擬譜微分矩陣，我們得到

(2.25)

其中。的元素定義如下：

由定義可知是反對稱矩陣。

從計算的角度，對傅立葉擬譜方法的評價是通過使用FFT算法，而不是譜差分矩陣。

多辛擬譜離散具有半離散多辛守恒律：

(2.26)

其中只需?。?/p>

(2.27)

(2.27)是多辛守恒律(2.2)的譜離散形式。因此傅立葉擬譜離散保留了方程(2.1)的多辛結構。

反對稱矩陣保證了整體辛守恒。因為和，對空間方向j加和，方程(2.26)變?yōu)?/p>

(2.28)

由(2.28)可知在時間上也具有整體辛守恒。

把Euler中點格式應用到半離散格式(2.25)中，可以得到

(2.29)

格式(2.29)是多辛的，換句話說，它滿足N次離散多辛守恒律：

(2.30)

其中只需?。?/p>

消去引入的變量u，v，p，q，r和s，可以得到：

(2.31)

其中。截斷誤差階為，m是解的光滑度。

第三章數(shù)值實驗

在本章中，將上述的五種格式進行分類對比，Eulerbox格式(簡記box_ls)、CrankNicolson格式(簡記CNI)與多辛Preissman(簡記E_LS)進行比較，當dt/dx分別取0.002/0.2，0.001/0.1，0.0005/0.05比較上述格式的無窮范數(shù)誤差及階數(shù)；同時，時間分裂步方法(簡記TSS)與Fourier擬譜方法(簡記ls1)比較，當dt=0.0001時，N=128，256，512和N=128時，dt=0.05，0.02，0.01的u和v的誤差和階數(shù)。比較上述格式的無窮范數(shù)誤差及階數(shù)。其中dx=L/N，L為距離長度，N為離散節(jié)點。數(shù)值試驗時，L不變，dx是隨著N的增加而減小。

給出誤差的階數(shù)的定義：

U_Order≈(lg(eu1/eu2))/(lg(h1/h2)),

V_Order≈(lg(ev1/ev2))/(lg(h1/h2))

h1、h2是兩次不同數(shù)值試驗的空間步長，即上述定義的dx；eu1、eu2是相對應的u的誤差，同理ev1、ev2是相對應的v的誤差。

3.1E_LS、CNI與Box_ls的誤差及階數(shù)

1、多辛Preissman在第二章已經(jīng)推導出相應的格式：

這是非線性方程組，采取迭代方法。

考慮第n+1層的計算方法(前n層已經(jīng)算好)：

Step1.初始值的選取.

Step2.進行迭代，下面是第k+1步的結果.

Step3.迭代終止條件.

當時迭代結束.

2、CNI的有關推導在第二章已經(jīng)給出

其中矩陣形式在前面章節(jié)有具體的說明，這里不再重復。

下面分別取下面三組數(shù)值：

dt=0.002,L=40,N=201,dx=L/N=0.2

dt=0.001,L=40,N=401,dx=L/N=0.1

dt=0.0005,L=40,N=801,dx=L/N=0.05

利用上述e_ls、CNI和box_ls三種數(shù)值方法分別求出u和v的無窮范數(shù)誤差、二范數(shù)誤差以及相對應的階數(shù)，eu、ev分別代表u、v的誤差,Order是對應的誤差的階數(shù)。

經(jīng)過運算，得到以下結果：

dt/dx

SchemeⅠ(box_ls)

Order

SchemeⅡ(CNI)

Order

0.002/0.2

eu=1.2176e-002

ev=2.8485e-003

eu=4.3666e-003

ev=1.7120e-002

0.001/0.1

eu=3.1101e-003

ev=6.5453e-004

U_Order=1.9690

V_Order=2.1217

eu=1.0934e-003

ev=4.2602e-003

U_Order=1.9977

V_Order=2.0067

0.0005/0.05

eu=7.8168e-004

ev=1.6213e-004

U_Order=1.9923

V_Order=2.0133

eu=2.7314e-004

ev=1.0751e-003

U_Order=2.0011V_Order=1.9865

表1box_ls格式與CNI格式在dt/dx取表中三組值的無窮范數(shù)誤差及階數(shù)

Dt/dx

SchemeⅡ(CNI)

Order

SchemeⅢ(e_ls)

Order

0.002/0.2

eu=4.3666e-003

ev=1.7120e-002

eu=1.4137e-002

ev=3.9281e-002

0.001/0.1

eu=1.0934e-003

ev=4.2602e-003

U_Order=1.9977

V_Order=2.0067

eu=3.6487e-003

ev=9.5126e-003

U_Order=1.9540

V_Order=2.0459

0.0005/0.05

eu=2.7314e-004

ev=1.0751e-003

U_Order=2.0011V_Order=1.9865

eu=9.9859e-004

ev=2.3104e-003

U_Order=1.9694

V_Order=2.0417

表2CNI格式與e_ls格式在dt/dx取表中三組值的無窮范數(shù)誤差及階數(shù)

dt/dx

SchemeⅠ(box_ls)

Order

SchemeⅡ(CNI)

Order

0.002/0.2

eu=1.7294e-002

ev=3.1941e-003

eu=7.1059e-003

ev=1.6248e-002

0.001/0.1

eu=4.3933e-003

ev=7.8758e-004

U_Order=1.9769

V_Order=2.0199

eu=1.7749e-003

ev=4.0883e-003

U_Order=2.0013

V_Order=1.9907

0.0005/0.05

eu=1.1040e-003

ev=1.9646e-004

U_Order=1.9926

V_Order=2.0032

eu=4.4421e-004

ev=1.0249e-003

U_Order=1.9984

V_Order=1.9960

表3box_ls格式與CNI格式在dt/dx取表中三組值的二范數(shù)誤差及階數(shù)

dt/dx

SchemeⅡ(CNI)

Order

SchemeⅢ(e_ls)

Order

0.002/0.2

eu=7.1059e-003

ev=1.6248e-002

eu=1.9566e-002

ev=4.1638e-002

0.001/0.1

eu=1.7749e-003

ev=4.0883e-003

U_Order=2.0013

V_Order=1.9907

eu=4.9466e-003

ev=9.9677e-003

U_Order=1.9838

V_Order=2.0626

0.0005/0.05

eu=4.4421e-004

ev=1.0249e-003

U_Order=1.9984

V_Order=1.9960

eu=1.3139e-003

ev=2.4099e-003

U_Order=1.9126

V_Order=2.0483

表4CNI格式與e_ls格式在dt/dx取表中三組值的二范數(shù)誤差及階數(shù)

3.2TSS與LS1的誤差及階數(shù)

1、時間分裂步方法(TSS)

用TSS進行數(shù)值試驗解耦合長短波方程，可分為以下幾個步驟，這也是用matlab編程求得結果的相應的算法。都表示對做Fourier變換，表示對做Fourier逆變換。針對此問題步驟如下：

Step1.方程組(1.25)的第一個方程可分解為：

Step2.對上述方程組的第一個方程做傅里葉變換得：

得，做逆變換進一步解得

Step3.對方程組(1.25)的第二個方程做傅里葉變換得：

然后對上述方程做差分得：

進而可求出，

Step4.直接解方程組的第二個方程，可得

Step5.最后求出

2、Fourier擬譜格式

在上一章針對耦合長短波的問題，利用Fourier擬譜方法推導出了下面的格式。

其中，是一階的傅立葉擬譜微分矩陣，的元素定義如下：

下面分別取下面三組數(shù)值：dt=0.0001時，N=128，256和512和N=128時，dt=0.05,0.02,0.01的u和v的誤差和階數(shù)。

利用上述TSS和ls1兩種數(shù)值方法分別求出u和v的無窮范數(shù)誤差、二范數(shù)誤差以及相對應的階數(shù)。與前面表示一致，eu、ev分別代表u、v的誤差。

N

SchemeⅣ(ls1)

SchemeⅤ(TSS)

128

eu=0.5926e-007

ev=1.1283e-005

eu=1.0862e-005

ev=5.0472e-005

256

eu=5.4559e-007

ev=2.5974e-006

eu=1.0798e-005

ev=4.6900e-005

512

eu=5.9126e-007

ev=2.6403e-006

eu=1.0812e-005

ev=4.6949e-005

表5當dt=0.001, N=128,256,512的ls1格式與TSS格式的無窮范數(shù)誤差

dt

SchemeⅣ(ls1)

SchemeⅤ(TSS)

0.001

eu=4.5926e-007

ev=1.1283e-005

eu=1.0862e-005

ev=5.0472e-005

0.0005

eu=9.2427e-007

ev=2.1668e-005

eu=5.4051e-005

ev=2.3803e-004

0.0002

eu=4.8666e-007

ev=1.3880e-005

eu=2.1660e-005

ev=9.7365e-005

表6當N=128, dt=0.02,0.01,0.005的ls1格式與TSS格式的二范數(shù)誤差

第四章結論

本文針對耦合的長短波方程，構造了Euler_box格式、多辛Pressiman格式和Fourier擬譜格式，從理論推導到數(shù)值實驗，我們可以得到以下結論：

1、從數(shù)值實驗的結果來看，所構造的三種格式計算精度都非常高，EulerBox格式是半顯式的，空間方向二階精度時間方向一階精度的數(shù)值方法；Preissman格式是隱式的，空間和時間方向都是二階精度的數(shù)值方法；Fourier擬譜方法是空間方向譜階精度時間方向二階精度。

2、數(shù)值實驗驗證我們的理論分析結果，新構造的三種數(shù)值算法都保持Hamilton系統(tǒng)的辛結構。

3、如果在長時間的積分區(qū)間內(nèi)計算，Hamilton系統(tǒng)能量誤差就可以控制在很小的范圍內(nèi)，充分顯示了多辛算法在守恒量的長時間數(shù)值模擬上的突出優(yōu)越性；

辛和多辛算法自創(chuàng)始以來已有一二十年，它的發(fā)展十分迅猛，特別是辛算法已發(fā)展得比較成熟。辛算法和多算方法還具有很大的研究價值，還存在著許多有待解決和完善的問題。本文雖然在這方面作了一些探討，但是仍然留有一些問題有待進一步研究。

參考文獻

[1]D.J.Benney,Ageneraltheoryforinteractionsbetweenshortandlongwaves,Stud.Appl.Math.,56(1977)81-94.

[2]M.Funakoshi,M.Oikawa,Theresonantinteractionbetweenalonginternalgravitywaveandasurfacegravitywavepacket,J.Phys.Soc.Japan,52(1983)1982-1995.

[3]N.Yajima,M.Oikawa,Formationandinteractionofsonic-Langmuirsolitons:inversescatteringmethod,Progr.Theoret.Phys.,56(1976)1719-1739.

[4]T.Kawahara,N.Sugimoto,T.Kakutani,Nonlinearinteractionbetweenshortandlongcapillary-gravitywaves,J.Phys.Soc.Japan,39(1975)1379-1386.

[5]Y.C.Ma,Thecompletesolutionofthelong-wave-short-waveresonanceequations,Stud.Appl.Math.,59(1978)201-221.

[6]Y.C.Ma,L.Redekopp,Somesolutionspertainingtotheresonantinteractionoflongandshortwaves,Phys.Fluids.,22(1979)1872-1876.

[7]M.Tsutsumi,S.Hatano,Well-posednessoftheCauchyproblemforthelongwave-shortwaveresonanceequations,NonlinearAnal.,22(1994),155-171.

[8]J.E.Marsden,G.P.Patrick,S.Shkoller,Multisymplecticgeometry,variationalintegra-torsandnonlinearPDEs,Comm.Math.Phys.,199(1999)351-395

[9]T.J.Bridges,Multi-symplecticstructuresandwavepropagation.MathematicalProceedingsoftheCambridgePhilosophicalSociety,121(1999)147-190

[10]T.J.Bridges,S.Reich,Multi-symplecticintegrators:numericalschemesforHamiltonianPDEsthatconservesymplecticity,Phys.Lett.A,284(2001)184-193.

[11]T.J.Bridges,S.Reich,MultisymplecticspectraldiscretizationsfortheZakharov-Kuznetsovandshallowwaterequations,PhysicaD,152(2001)491-504.

[12]S.Reich,Multi-symplecticRunge-KuttacollocationmethodsforHamiltonianwaveequation,J.Comput.Phys.,157(2000)473-499.

[13]B.Moore,S.Reich,MultisymplecticintegrationmethodsforHamiltonianPDEs,FutureGener.Comput.Syst.,19(2003)395-402.

[14]J.B.Chen,M.Z.Qin,Multi-symplecticFourierpseudospectralmethodforthenonlinearSchrodingerequation,Electron.Trans.Numer.Anal.,12(2001)193-204.

[15]A.L.Islas,D.A.Karpeev,C.M.Schober,GeometricintegratorsfortheonlinearSchrodingerequation,J.Comput.Phys.,173(2001)116-148.

[16]A.L.Islas,C.M.Schober,MultisymplecticmethodsforgeneralizedSchrodingerequa-tions,FutureGener.Comput.Syst.,19(2003)403-413.

[17]A.L.Islas,C.M.Schober,Onthepreservationofphasespacestructureundermultisymplecticdiscretization,J.Comput.Phys.,197(2004)585-609.

[18]J.Wang,MultisymplecticFourierpseudospectralmethodforthenonlinearschrodingerequationswithwaveoperator,J.Comput.Math.,25(2007)31-48.

[19]J.Wang,AnoteonmultisymplecticFourierpseudospectralmethodforthenonlinearSchrodingerequation,Appl.Math.Comput.,191(2007)31-41

[20]J.Wang,MultisymplecticintegratoroftheZakharovsystem,Chin.Phys.Lett.,25(2008)3531-3534

[21]秦孟兆,任意階精度蛙跳格式穩(wěn)定性分析,計算數(shù)學,1(1992)1-9

[22]M.Z.Qin,W.J.Zhu,ConstructionofsymplecticschemesforwaveequationsviaHyperbolicFunctionssinh(x),cosh(x)andtanh(x),JComputersMath.Application,26(8)(1993)1-11.

[23]M.Z.Qin,AsymplecticdifferenceschemeforthePDEs,J.AdvancedMathematicss,3(1997)349-354.

[24]D.B.Dunean,SymplecticfinitedifferenceapproximationsofthenonlinearKlein-Gordonequation,SIAMJ.Number.Anal.34(1997)1742-1760.

[25]D.L.Wang,SemidisereteFourierspectralapproximationsofinfinitedimensionalHamiltoniansystemsandconservationlaw,J.ComputersMath.Application,21(4)(1991)63-75.

[26]J.E.Marsden,G.P.Patriek,S.Shkoller,multisymplecticgeometry,variationalintegrators,AndnonlinearPDEs,Comm.Math.Phys,199(1998)351-395.

[27]J.E.Marsden,S.Shkoller,MultisymplecticgeometrycovariantHamiltonianandwaterwaves.Math.Proc.Camb.Philos.Soc,125(1999),553

[28]J.E.Marsden,G.P.Patriek,S.Shkoller,Variationalmethods,multisymplecticgeometryandmechanics,J.Geom.Phys.,38(2001)253-284.

[29]J.E.Marsden,M.West,DiscretemechanicsandvariationalintegratorsActNumerical,10(2001)357-514.

[30]T.J.Bridges,Multisymplecticstructuresandwavepropagation.MathematicalProceedingsOftheCambridgePhilosophicalSociety121(1987)147-190.

[31]GilbertStrang,Ontheconstructionandcomparisonofdifferenceschemes,SIAMJ.Numer.Anal.,1968,5(3):506-517.

[32]RandallJ.LeBequeandJosephOliger,Numericalmethodsbasedonadditivesplittingsforhyperbolicpartialdifferentialequations,MathematicsofComputation.,1983,40(162):469-497.

[33]AntonArnoldandChristianRinghofer,AnoperatorsplittingmethodfortheWignerPoissonproblem,SIAMJ.Numer.Anal.,1996,33(4):1622-1643.

[34]TaoTang,Convergenceanalysisforoperator-splittingmethodsappliedtoconservationlawswithstiffsourceterms,SIAMJ.Numer.Anal.,1998,35(5):1939-1968.

[35]RandallJ.LeVeque,Intermediateboundaryconditionsfortime-splitmethodsappliedtohyperbolicpartialdifferentialequations,MathematicsofComputation.,1986,47(175):37-54.

[36]C.A.J.Fletcher,Time-splittingandthefiniteelementmethod,UniversityOfSydney,Sydney,NSW,2006,37-47.

[37]J.A.C.WeidemanandB.M.Herbst,Split-stepmethodsforthesolutionofthenonlinearSchrodingerequations,SIAMJ.Numer.Anal.,1986,23(3):485-507.

[38]QianshunChang,Yau-ShuWongandChi-KunLin,Numericalcomputationsforlong-waveshort-waveinteractionequations,JournalofComputationalPhysics.,2008,227:8489-8507.

[39]R.Zhang,K.ZhangandY.S.Zhou,Numericalstudyoftime-splittingandspace-timeadaptivewaveletschemeforSchrodingerequations,JournalofComputationalandAppliedMathematics.,2006,195:263-273.

[40]HanquanWang,Numericalstudiesonthesplit-stepfinitedifferencemethodfornonlinearSchrodingerequations,AppliedMathematicsandComputation,2005,170:17-35.

[41]G.M.Muslu,H.A.Erbay,Higher-ordersplit-stepFourierschemesforthegeneralizednonlinearSchrodingerequation,MathematicsandComputersinSimulation.,2005,67:581-595.

[42]WeizhuBao,ShiJinandPeterA.Markowich,Ontime-splittingspectralapproximationsfortheSchrodingerequation,JournalofComputationalPhysics.,2002,175:487-524.

[43]WeizhuBao,ShiJinandPeterA.Markowich,Numericalstudyoftime-splittingspectraldiscretizationsofnonlinearschrodingerequationinthesemiclassicalregimes,SIAMJ.sci.Compt.,2003,23:27-64.

[44]Hanquan