上海交通大学硕士学位论文

目录

第一章引言

-2-

1.1耦合长短波方程的背景

-2-

1.1.1Schrodinger方程的介绍

-2-

1.1.2Kdv方程的介绍

-8-

1.1.3Schrodinger-KdV方程组及耦合长短波方程的由来

-10-

1.2Hamilton系统、辛算法及多辛算法

-11-

1.3原有的数值方法

-17-

1.3.1时间分裂方法

-17-

1.3.2Crank-Nicolson方法

-19-

第二章多辛格式

-22-

2.1EulerBox格式

-24-

2.2Preissman格式

-30-

2.3Fourier拟谱格式

-35-

第三章数值实验

-40-

3.1E_LS、CNI与Box_ls的误差及阶数

-40-

3.2TSS与LS1的误差及阶数

-49-

第四章结论

-55-

参考文献

-56-

致谢

-59-

第一章引言

1.1耦合长短波方程的背景

1.1.1Schrodinger方程的介绍

薛定谔方程(Schrodingerequation)是由奥地利物理学家薛定谔提出的量子力学中的一个基本方程，也是量子力学的一个基本假定，其正确性只能靠实验来检验。是将物质波的概念和波动方程相结合建立的二阶偏微分方程，可描述微观粒子的运动，每个微观系统都有一个相应的薛定谔方程式，通过解方程可得到波函数的具体形式以及对应的能量，从而了解微观系统的性质。

一．关于波的简单介绍

波的形式是多种多样的，平面波是指在波的传递过程中在一段距离内波幅变动不大的波。描述平面波的状态的波函数为：

(1.1)

为波函数。x与t分别表示x方向的距离变量与时间变量。函数计算结果的物理量为波幅，即长度单位为米。

A为波幅，单位为米。指平面波相对稳定的最大波幅绝对值。

cos(kx--wt)为波幅的变动系数，数值在-1，0，+1之间变动。

k为单位长度弧度数，即，单位为弧度／米。

为单位时间弧度数，即，单位为弧度／秒。

由于=波速度，=即波在传播过程中不同距离的时间差。波函数(1.1)式可以写成：

(1.2)

如果上述波函数(1.2)式去掉括号中的第一项，得到：

(1.3)

即为一个以时间为变量的振幅函数，这样的波并不向外传播，而是在原地上下振动，余弦cos的正负角的函数值是相等的，因此括号中的负号可以忽略。

如果波函数(1.2)式中去掉括号中的第二项同时保留第一项，得到：

(1.4)

即为一个以坐标为变量的函数，是描述波到达各种位置的振动位移情况。

波函数(1.1)式也可以写成以e为底的复数形式：

(1.5)

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

(1.6)

化简为一维形式的波动方程：

(1.7)

函数是直接描述变量之间的关系，而物理方程则是在复杂的情况下表述各种物理量的关系，它可以包含各种函数，但对函数中的变量之间关系的表述不是那么直接，函数表达式与方程之间，在一定条件下可以转换。下面给出具体的推导过程。

对宏观平面波函数(1.5)式进行转换，得：

(1.8)

对时间二次微商得：

整理得：

对坐标二次微商得：

整理得：

将两个二次微商的结果连接得：

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

这就是前面的一维形式的波动方程(1.7)式。可见波函数与波动方程在一定条件下可以相互转换。应当注意到由波函数转换为波动方程的过程中，波幅常量A消失了，说明这个波动方程中含有的波函数与波幅常量没有直接联系。如果由波动方程反解出波函数，这个常量A需经过归一化处理重新找回来。

通常，参照一般的波动描述方式，似乎可以建立一个类似的电磁波的波函数表达式：

然而，我们以宏观的平面波方式来描述量子波时，必须考虑到量子波与平面波的重要区别：平面波的波函数描述的是一个波群(波束)，而量子波由于其量子性，表现为个体的量子波。即使一个有限的空间在一段有限的时间里，只有一个电子或者一个光子在运动，我们也得承认这个电子或者光子的波动性，应当也有波函数描述它。因此用波函数来描述量子波时，显然应当与描述光束或者电子束的方式有所区别，或者说不能用描述一般平面波的波函数的方法来描述量子波。

分析光子的运动：光子在振动且以光速运动，以位置作为变量，光子的相位随之而变，换以时间作为变量，光子的相位也是随之而变。光子的波幅、一个波动周期内的速度也是在随位置和时间变动的，但是波动量的大小与相位的变动是相关的。因此首先选取相位作为一个波动周期的因变量。

将原来描述光束的函数式(1.8)改为描述单个光子的波函数：

(1.9)

函数的值就是相位，单位为弧度。

将上述描述光子相位的函数用三角函数表示：

(1.10)

得到的是无量纲数，数值在-1，0，1之间变动，这就是我们描述概率幅的函数。加上适当的常数项A，使得波函数有适当的物理量，将模平方并按一定的物理量积分后，得到的结果是有量纲数，即概率密度：

(1.11)

根据(1.8)式对电磁波的函数表达式，按照式(1.6)的通用波动方程，可以建立一个电磁波的波动方程：

(1.12)

方程中E表示什么呢?根据前面对一般波函数与波动方程的转换分析我们知道，这个E无法代表电场强度，因为由波函数转换为波动方程的过程中，电场强度常量E已经消失了，现在这个E函数计算得到的值是无量纲数，在一定条件下按归一化处理后可以有一定的物理量纲，就是概率幅。由此可见，即使人们将波动视为量子的出现概率的似波性，在进行归一化处理之前，波动方程并不能完整地描述量子的波动状态，因为波动方程计算出来的概率幅数值在-1,0,1之间变动，必须考虑行程、波数等因素，进行适当的转换后才能得到实际波动的概率幅数值。因此归一化处理，实际上是使波动方程描述的物理现象回归实在的波动。笔者将在后面的实例中进一步说明，波动方程加上波动的边界条件，再加上归一化处理，是可以描述量子实在的。

二．薛定谔方程

建立描述电子运动的波函数

既然电子的运动可以形成波，那么应当如同描述光子的波动一样，有一个波动方程描述电子的波动。参照前面描述光子的波动方程(1.10)式建立描述电子的波函数：

上式是量子单位制形式。将上面的波函数按下式描述为经典形式：

由于动量，即动能／速度。能量E=，此处表示电子库仑势能。普朗克常数h=1，

(1.13)

将电子的波函数改写为复数形式并进行微商

将(1.13)式改写为复数形式：

(1.14)

这里为了波函数的完整性加上了A，在后面进行微商处理时A自然消失，而在归一化处理时A又回到公式中，那时它的物理意义与我们进行归一化处理方式有关。

将上述函数对时间微商得：

(1.15)

将上述函数对坐标微商得：

(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)

如果不考虑时间因素，得到三维定态的薛定谔方程：

(1.23)

其中：

简化为一维定态的薛定谔方程：

(1.24)

定态是指粒子的概率密度分布不随时间变化的状态，定态薛定谔方程示，任何时候粒子概率幅的空间分布，不随时间变化。

1.1.2Kdv方程的介绍

1834年英国科学家ScottRussell偶然观察到了一种奇妙的水波。1844年，他在《英国科学促进协会第14届会议报告》上发表的《论波动》一文中，对此现象作了生动的描述：“我观察过一次船的运动，这条船被两匹马拉着沿狭窄的运动迅速前进着，突然，船停了下来，而被船所推动的大堆水却并不停止，它们积聚在船头周围激烈地扰动着，然后水浪突然呈现在一个滚圆而平滑、轮廓分明的巨大孤立波峰，它以巨大的速度向前滚动着，急速地离开了船头，在行进中它的形状和速度并没有明显的改变，我骑在马上紧跟着观察，它以每小时约八、九英里的速度滚滚向前，并保持着长约30英尺，高约1-1.5英尺的原始形状，渐渐地它的高度下降了。当我跟踪1-2英里后，它终于消失在逶迤的河道之中”。这就是Russell观察到的奇特现象，进而他认为这种孤立的波动是流体运动的一个稳定解，并成它为“孤立波”。Russell当时未能成功地证明并使物理学家们信服他的论断，从而埋汰数学家未能从已知的流体运动方程预言出这一现象，之后有关孤立波的问题在当时许多物理学家中引起了广泛的争论。直到60年后的1895年，Kerteweg，DeVries研究了浅水波的运动，在长波近似和小振幅的假定下，建立了单向运动的浅水波运动方程

下面就简单推导一个Kdv方程：

在竖直平面考察平面内流体的运动，考虑下面定解问题

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

引进小参数a是平面波振幅，是波长。作下列变换：

则上述方程组变为：

把去掉还原成：

用构造

所以成立

构造的满足因为

将代入到得

将代入到得

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

在上述两式中忽略一次以及以上的项得：

取可得，展开

在上述两式中忽略二次以及以上的项得

将代入到刚得到的两项中分别可得：

取可得

这就是经典的Kdv方程。

1.1.3Schrodinger-KdV方程组及耦合长短波方程的由来

各种客观环境已经在研究长短波之间的相互作用现象，尤其是在流体力学、等离子物理和化学物理。短波通常用Schrodinger方程来描述，长波一般用带有色散的某种类型的波方程来描述。本文中，我们用多辛积分[1,2,3]格式来研究耦合长短波方程(LSIE)。

(1.25)

满足

(1.26)

复值函数表示短波的包络层，实值函数表示长波的振幅。是实常数。在方程(1.25)被用来在重力和微小范围模式下建立表面波。这个方程是Schrodinger-KdV方程组的一个特殊例子。

(1.27)

关于初始边界问题(1.25)--(1.26)，可以得到下面命题

命题1.1带有初始边界问题(1.25)--(1.26)的解满足下面的守恒律

是复函数φ的共轭函数。

方程组(1.25)具有独特的数学特点因为研究发现它具有完整的可积结构。特别的，它有一个逆散射变换和显示的N孤子解[5,6]。各种数值技术尤其是有限元方法、谱方法和时间分裂步方法都已经用来解各种形式的长短波相互作用的方程。

1.2Hamilton系统、辛算法及多辛算法

力学研究中一个非常重要的是对称的几何观点，它不论是从基本的原理公式出发，还是到具体的应用，都特别强调了几何方法和力学研究的密不可分的关系。其中Hamilton提出的力学定理可以用微分流形来研究刚体体系及太阳系等复杂系统的力学性质；可以用相应的Hamilton函数的对称性来研究能量、线性动量与角动量等Hamilton系统的守恒性质。

Hamilton系统是一种重要的力学系统，广泛的出现在物理、力学、工程、纯数学与应用数学领域。通常可以认为，一切耗散可忽略不计的真实物理过程都可以表示成Hamilton方程的形式，而它们的共同基础都是辛几何。辛几何历史可以追溯到十九世纪英国物理学家和数学家Hamilton，他为了研究Newton力学，引入广义坐标和广义动量来表示系统的能量，即Hamilton函数。

考虑Newton运动方程，设表示质点位置，它满足拉格朗日方程其中表示动能，表示势能。引入共轭动量，其满足

下面定义则

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

其中，为n阶单位矩阵，称为有限维Hamilton系统，H(z)为Hamilton函数。

Hamilton系统有两个重要特性：守恒性与辛结构。我们知道Hamilton系统的解是一个单参数的保测变换,即辛变换。因此在研究Hamilton系统的计算方法时候，都要使得离散后的方程保持原有系统的辛结构。用传统的数值算法模拟Hamilton系统，会破坏系统的辛结构，进而使数值模拟不能保持长时间的稳定。因此构造一种能保持Hamilton系统的辛结构的算法，具有重要意义，称能保持Hamilton系统辛结构的算法为辛算法。

无限维Hamilton系统一般情况下可以表示为：

(1.28)

其中，z=z(x,t)为状态向量，H为Hamilton函数：

(1.29)

为H的Frechet导数。

等式(1.28)也可通过Lagrange泛函导出。考虑如下非线性Klein-Gordon方程，其中，为某一光滑非线性函数。

对以下的Lagrange泛函作变分

(1.30)

其中Lagrange密度为：

就可以得到Euler-Lagrange方程：

对Lagrange密度L作Legendre变换，即令得到一阶方程组

(1.31)

令

则(1.31)等价于Hamilton系统(1.28)。相应的辛二形式为：

从以上过程可知，Legendre变换仅对时间方向进行。若对空间方向也进行Legendre变换，即可得到相应一阶方程

(1.32)

Hamilton函数为：

相应空间方向的辛二形式为：

大量的孤立波方程，可表示为无限维Hamilton系统。例如：

非线性Schrodinger方程：

其Hamilton函数为：

其中

KdV方程：

其Hamilton形式为：

将求解有限维Hamilton系统的辛算法推广到无限维最直接有效的方法是：先对空间方向进行离散，这样离散以后的系统为有限维Hamilton系统，然后应用辛算法。常用的半离散方法为差分法[20,21,22,23,24]，拟谱方法和配置点法[25]等。

在上例导出非线性Klein-Gordon方程的无限维Hamilton形式时，采用的是不完全Legendre变换。若对方程(1.29)采用完全的Legendre变换，即同时令，可得到相应一阶方程组：

(1.33)

记(1.33)等价于：

(1.34)

其中

Bridges和Reich称形如(1.34)的方程为具有多辛结构的Hamilton系统。

求解有限维Hamilton系统的辛算法与其它方法相比具有许多优点，其中有一点就是辛算法可以对守恒量长时间数值模拟。但是当用辛算法对无限维Hamilton进行离散时具有局限性，具体表现在守恒量是全局性的。为了克服此局限性，Marsden等[26,27,28,29]和Bridges[30,31,32]分别从Lagrange系统和Hamilton系统出发，提出了多辛Hamilton系统和多辛算法的概念。Marsden等是从变分原理开始，由边界项得到辛形式，限制到具体方程上得到该方程的多辛守恒律；Bridges是将有限维的Hamilton系统推广到无限维Hamilton系统，使得偏微分方程在时间方向和空间方向上都有各自不同的辛结构。

多辛结构的良好特征是它可以导出一个包含微分二形式的守恒律，即多辛守恒律。与多辛守恒律相对应的是多辛Hamilton系统具有能量守恒律和动量守恒律。由于这些守恒律不依赖边界，因而都是局部守恒律。

Bridges和Reich考虑了一般情形的多辛Hamilton系统。设是有限维相空间，S为Q上的函数，M，K为两个反对称矩阵，则称

(1.35)

为多辛Hamilton系统。其中，为反对称矩阵；为状态向量及为光滑函数；表示函数关于z的梯度，为底空间的两个独立变量。

具有m个空间变量的一般多辛Hamilton系统可表示为：

(1.36)

洪佳林等对(1.36)作了详细的讨论[33,34,35]，为简单起见，我们仅讨论具有一个空间变量，且系数矩阵为常系数的情形，即(1.35)中。

Bridges和Reich证明了多辛Hamilton系统(1.35)具有如下局部守恒律，即多辛守恒律、局部能量守恒律和局部动量守恒律。

Hamilton系统(1.35)满足多辛守恒律：

(1.37)

其中分别对应于t和x方向的两个不同的辛结构。

多辛Hamilton系统(1.35)具有局部能量守恒律：

(1.38)

和局部动量守恒律：

(1.39)

利用矩阵分解：

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

与无限维Hamilton系统比较，多辛Hamilton系统具有的守恒律都是局部的，但是它们在任意的时空区域内成立，而且不依赖于边界条件。若对局部守恒律在空间方向积分，并利用适当的边界条件，就可导出相应的整体守恒律即辛守恒律、能量守恒律和动量守恒律。而这些适当的边界条件显然也是整体守恒律的必要条件。因此，局部守恒律蕴涵着整体守恒律。另外多辛Hamilton系统的相空间是有限维的，因而可应用有限维的理论来分析方程的性质。

1.3原有的数值方法

1.3.1时间分裂方法

时间分裂方法[31-35]是近些年研究非线性偏微分方程数值解时使用较多的方法，目前常见的有分裂差分方法[36-40]、分裂谱方法[41-46]等等，时间分裂方法的主要优点是把方程的非线性项与线性项分开计算，把较难处理的非线性项单独处理。

时间分裂方法的基本原理[37]如下：

对方程：

(1.40)

方程(1.40)改写为如下形式：

(1.41)

其中分别为方程的线性项和非线性项。

在时间步上，由(1.41)有：

(1.42)

即有：

(1.43)

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

(1.44)

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

(1.45)

(1.46)

其中为函数在时间层对u的逼近。

因此方程(1.40)可以采用下面的时间分裂步来进行计算：

分裂步1：(1.47)

分裂步2：(1.48)

方程(1.47)是一个非线性的常微分方程，可以精确的求解，方程(1.48)可以采用不同的数值方法来计算，常见的有分裂谱方法，分裂差分方法，如分裂C-N格式等等。

分裂格式(1.47)-(1.48)在时间方向上为一阶的。根据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.最后求出：

这一小节先是对时间分裂步方法的基本原理的一个简单介绍，然后简单推导得出时间分裂步方法在本篇文章当中对耦合长短波方程的实现步骤，当然以上五个步骤也是在用matlab进行数值试验时的基本编程思路。在本文的数值试验部分，取特定的参数及变量值就可以按照上述步骤编写程序得出相应的结论。

1.3.2Crank-Nicolson方法

对方程组(1.25)进行以下的差分离散：

记如下：

其中下标N是时间节点数。

将上述方程组的第一个等式按照定义的合并空间指标j之后得到下面这个等式：

其中A是单位矩阵，B和运算是如下格式：

上式两边同时乘以：

化简得：

令上式可变为：

同理方程组第二式可化为：

其中

综上所述最后推出迭代格式为：

本文研究的方程是耦合的长短波方程，它是Schrodinger-KdV方程组的一个特例。本文所做的工作主要是建立耦合的长短波方程的几种多辛格式，并通过数值试验与已有的数值算法比较。

在本文的第一章引言中，会介绍一下耦合的长短波方程的背景，然后叙述Hamilton系统、辛算法和多辛算法，进而会介绍两种现在认可的数值算法：时间分裂方法和Crank-Nicolson方法，并对相应的格式进行简单的推导。

在本文的第二章中，会详细叙述新建立的三种多辛格式，它们是EulerBox格式、多辛Preissman格式和Fourier拟谱格式，给出相应的性质并格式的误差的阶数做了详细的说明。三种多辛数值算法用来解决周期性带有初始问题的长短波相互作用方程。

第三章是本文的数值试验，针对上述5种数值格式，取三组时间步长dt和相应的节点数N算出各自的误差，并根据公式求出相应的误差的阶数。

第四章是本文的结论说明，通过理论分析和数值实验来分析新构造的三种格式，给出本篇文章的结论。

第二章多辛格式

根据多辛的定义[9]，如果偏微分方程可表示为：

(2.1)

则称方程（2.1）为多辛Hamilton系统。其中：M,K是反对称矩阵；z(x，t)是状态变量的函数。S：是光滑函数；是Hamiltonian函数S=S(z)的梯度。

多辛格式满足多辛守恒律

(2.2)

是外微分格式

(2.3)

定义了一个对称的时间空间结构。多辛结构由Nothers定理[8]自然地可以得到局部守恒律。事实上，对于Hamiltonian系统中的偏微分方程(2.1)，当S与x，t无关，可以得到局部能量守恒律和局部动量守恒律。

(2.4)

(2.5)

对于周期边界或者在边界条件为0的F(z)和G(z)，可以推出整体能量和动量守恒律

(2.6)

其中，

能保持离散多辛守恒律的数值方法为多辛算法[9,10,11]。为了得到长短波相互作用方程的多辛格式，令方程组(1.25)可以用下列一阶方程组表示：

(2.7)

反对称矩阵M,K表示如下。

方程(2.1)等号右边的S(z)为：

直接计算就可以得出方程组(3.1)满足多辛守恒律：

(2.8)

(2.4)和(2.5)中定义的密度函数可以下面给出

(2.9)

加上适当的边界条件可以得到整体守恒律。例如在周期或者消失在无限边界条件的情况下可对E和I在空间区域上积分，就可以得到整体能量守恒律(1.6)和动量守恒律(1.7)。

(2.10)

对于长短波相互作用方程这是两个很重要的整体守恒律。

2.1EulerBox格式

根据[13]通过引入两个矩阵M和K的分裂矩阵可以获得方程组(2.1)的Euler-box格式。例如，分裂矩阵M,K如下：

(2.11)

相应的格式变为：

(2.12)

格式(2.12)满足离散多辛守恒律：

(2.13)

其中很明显矩阵的分裂是不唯一的。但是在本文中，分别是M,K的上三角矩阵。这种特殊的选择使得Euler-box格式变为以下形式：

(2.14)

离散方程组(2.14)也具有离散多辛守恒律(2.13)的形式。其中只需取为消去所有引入的变量p，q，r和s，与长短波相互作用方程组(1.25)相比可以得到下面的多辛积分：

(2.15)

其中，注意这个格式是线性隐式的，不需要每一步都解线性代数方程。可以证明此格式的局部截断误差是。

命题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系统时，满足离散局部能量和局部动量守恒律;对于非线性Hamilton系统，Moore和Reich证明了其仅满足半离散局部能量和动量守恒律。

由多种方法可导出多辛Preissman格式，其中最典型的方法为时间和空间方向分别利用隐式中点辛格式。多辛Preissman格式可表示为

其中

.

和分别表示空间步长和时间步长。Bridges和Reich证明了Preissman格式能精确满足多辛守恒律。

Preissman格式满足多辛守恒律

其中

对于线性Hamilton系统，多辛Preissman能精确保持局部能量和动量守恒律，如果Hamilton函数，其中A为对称矩阵，则Preissman格式全离散局部能量守恒律

和全离散局部动量守恒律

其中:

然而，一般而言，对非线性Hamilton系统，离散局部能量和动量守恒律不能精确满足。Moore和Reich提出了一种半离散能量和动量守恒律。

半离散格式

满足半离散能量守恒律

其中

而半离散格式

满足半离散动量守恒律

其中

由此可见，对于全离散多辛格式，它能精确保持Hamilton系统的离散多辛结构，但并不意味着系统的其它守恒量，如局部能量和动量守恒律，整体能量和动量，以及决定相空间结构的其它整体不变量守恒。而大量的数值计算显示，多辛格式能使局部守恒律在长时间计算中保持很好。

针对多辛Preissman格式，给出能量密度E，能量流F，动量密度I和动量流G相应的离散形式。分别如下

在空间和时间中都用辛隐式中点离散可以得到Preissman格式：

(2.16)

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

(2.17)

其中

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

(2.18)

(2.18)是一个二阶的多辛格式。所以它满足相应的多辛离散守恒律(2.17)，其中只需取：

(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拟谱格式

谱方法是求解微分方程的高精度算法，是近二十年来发展最快的数值方法之一，它已成功地应用于许多领域的数值计算，例如流体力学、理论物理、量子力学、非线性光学等。由于用无限可微的正交函数作为基函数，谱方法具有通常有限元和差分方法无与伦比的高精度，即如果精确解是解析的，则谱方法的精度是指数阶的。

对于满足周期边界条件的多辛Hamilton系统，Bridges和Reieh[43]提出了基于Fourier变换的多辛Fourier离散。多辛Fourier变换导出了在Fourier空间中的多辛概念和半离散系统。Islas和Schober[65,66]对多辛Fourier离散进行进一步研究，证明了对于线性Hamilton系统，能保证局部能量和动量守恒。陈景波和秦孟兆[44]讨论了在物理空间上进行离散的多辛Fourier拟谱方法，并应用于求解非线性Schrodinger方程。

Fourier拟谱方法分二个基本步骤。第一步为通过解在配置点上三角多项式插值构造解的离散形式；第二步利用解在配置点上的值求出导数值，即构造微分矩阵。为简单起见，设所讨论的区间。对任何整数N>0，令为空间步长，

为配置点；

为插值空间，其中为

插值算子助定义如下：对任意

利用

及正交性

得：

对定义离散内积和离散模

在配置点上的导数值可由函数值和微分矩阵求得

其中表示一阶Fourier微分矩阵，其元素为

显然，为反对称矩阵。

傅里叶变换可以保持偏微分方程的多辛性质不变。离散傅立叶方程恢复了标准谱离散的性质，这导致哈密顿常微分方程系统标准辛[11,14]可积。

对方程(1.27)用傅立叶拟谱方法，是一阶的傅立叶拟谱微分矩阵，我们得到

(2.25)

其中。的元素定义如下：

由定义可知是反对称矩阵。

从计算的角度，对傅立叶拟谱方法的评价是通过使用FFT算法，而不是谱差分矩阵。

多辛拟谱离散具有半离散多辛守恒律：

(2.26)

其中只需取：

(2.27)

(2.27)是多辛守恒律(2.2)的谱离散形式。因此傅立叶拟谱离散保留了方程(2.1)的多辛结构。

反对称矩阵保证了整体辛守恒。因为和，对空间方向j加和，方程(2.26)变为

(2.28)

由(2.28)可知在时间上也具有整体辛守恒。

把Euler中点格式应用到半离散格式(2.25)中，可以得到

(2.29)

格式(2.29)是多辛的，换句话说，它满足N次离散多辛守恒律：

(2.30)

其中只需取：

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

(2.31)

其中。截断误差阶为，m是解的光滑度。

第三章数值实验

在本章中，将上述的五种格式进行分类对比，Eulerbox格式(简记box_ls)、CrankNicolson格式(简记CNI)与多辛Preissman(简记E_LS)进行比较，当dt/dx分别取0.002/0.2，0.001/0.1，0.0005/0.05比较上述格式的无穷范数误差及阶数；同时，时间分裂步方法(简记TSS)与Fourier拟谱方法(简记ls1)比较，当dt=0.0001时，N=128，256，512和N=128时，dt=0.05，0.02，0.01的u和v的误差和阶数。比较上述格式的无穷范数误差及阶数。其中dx=L/N，L为距离长度，N为离散节点。数值试验时，L不变，dx是随着N的增加而减小。

给出误差的阶数的定义：

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

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

h1、h2是两次不同数值试验的空间步长，即上述定义的dx；eu1、eu2是相对应的u的误差，同理ev1、ev2是相对应的v的误差。

3.1E_LS、CNI与Box_ls的误差及阶数

1、多辛Preissman在第二章已经推导出相应的格式：

这是非线性方程组，采取迭代方法。

考虑第n+1层的计算方法(前n层已经算好)：

Step1.初始值的选取.

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

Step3.迭代终止条件.

当时迭代结束.

2、CNI的有关推导在第二章已经给出

其中矩阵形式在前面章节有具体的说明，这里不再重复。

下面分别取下面三组数值：

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三种数值方法分别求出u和v的无穷范数误差、二范数误差以及相对应的阶数，eu、ev分别代表u、v的误差,Order是对应的误差的阶数。

经过运算，得到以下结果：

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取表中三组值的无穷范数误差及阶数

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取表中三组值的无穷范数误差及阶数

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取表中三组值的二范数误差及阶数

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取表中三组值的二范数误差及阶数

3.2TSS与LS1的误差及阶数

1、时间分裂步方法(TSS)

用TSS进行数值试验解耦合长短波方程，可分为以下几个步骤，这也是用matlab编程求得结果的相应的算法。都表示对做Fourier变换，表示对做Fourier逆变换。针对此问题步骤如下：

Step1.方程组(1.25)的第一个方程可分解为：

Step2.对上述方程组的第一个方程做傅里叶变换得：

得，做逆变换进一步解得

Step3.对方程组(1.25)的第二个方程做傅里叶变换得：

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

进而可求出，

Step4.直接解方程组的第二个方程，可得

Step5.最后求出

2、Fourier拟谱格式

在上一章针对耦合长短波的问题，利用Fourier拟谱方法推导出了下面的格式。

其中，是一阶的傅立叶拟谱微分矩阵，的元素定义如下：

下面分别取下面三组数值：dt=0.0001时，N=128，256和512和N=128时，dt=0.05,0.02,0.01的u和v的误差和阶数。

利用上述TSS和ls1两种数值方法分别求出u和v的无穷范数误差、二范数误差以及相对应的阶数。与前面表示一致，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格式的无穷范数误差

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格式的二范数误差

第四章结论

本文针对耦合的长短波方程，构造了Euler_box格式、多辛Pressiman格式和Fourier拟谱格式，从理论推导到数值实验，我们可以得到以下结论：

1、从数值实验的结果来看，所构造的三种格式计算精度都非常高，EulerBox格式是半显式的，空间方向二阶精度时间方向一阶精度的数值方法；Preissman格式是隐式的，空间和时间方向都是二阶精度的数值方法；Fourier拟谱方法是空间方向谱阶精度时间方向二阶精度。

2、数值实验验证我们的理论分析结果，新构造的三种数值算法都保持Hamilton系统的辛结构。

3、如果在长时间的积分区间内计算，Hamilton系统能量误差就可以控制在很小的范围内，充分显示了多辛算法在守恒量的长时间数值模拟上的突出优越性；

辛和多辛算法自创始以来已有一二十年，它的发展十分迅猛，特别是辛算法已发展得比较成熟。辛算法和多算方法还具有很大的研究价值，还存在着许多有待解决和完善的问题。本文虽然在这方面作了一些探讨，但是仍然留有一些问题有待进一步研究。

参考文献

[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]秦孟兆,任意阶精度蛙跳格式稳定性分析,计算数学,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