《海洋数值模拟》实验2指导书

1目录微分方程数值实验.....................................................................错误!未定义书签。1.孤立波简介.............................................................22.有限差分格式..........................................................42.1蛙跳格式..............................................................52.2跳点格式..............................................................72.3分裂格式..............................................................83.实例...................................................................94.附录：MATLAB源码.......................................................122微分方程数值实验KDV方程是历史上最著名的方程之一，其起源是水波中的孤立波。1.孤立波简介我们自己先设想一下：在一条窄河道中，当迅速前进的船突然停下时，在船头形成的一个孤立的水波迅速离开船头,继续向前运动，如下图(1)所示。图1.孤立波的运动图孤立波是由是苏格兰一位优秀的造船工程师Russell(JohnScottRussell1808-1882)在1834年研究船舶在运动中所受到的阻力时发现的。在一次学术报告中,他是这样描述的：“我把注意力集中在船舶给予流体的运动上，立刻就观察到一个非同寻常而又非常绚丽的现象，它是如此之重要，以致我将首先详细描述它所表现出来的外貌。当我正在观察一只高速运动的船舶，让它突然停止时，在船舶周围所形成的小波浪中，一个紊乱的扰动现象吸引了我的注意。在船身长度的中部附近，许多水聚集在一起，形成一个廓线很清楚的水堆，最后还出现一尖峰，并以相当高的速度开始向前运动，到船头后,继续保持它的形状不变，在静止流体的表面上，完全孤立地向前运动，成为一孤立行进波,直到河道的转弯处才开始消失掉。”纵观力学和物理学的发展史不难看到，每当开始引入一种新思想或新概念的时候，往往会受到怀疑和非难，并常会引起激烈的争论，孤立波的命运也是如此。当时科学界的权威们对Russell的这些结果，一开始时就表示了怀疑和反对。甚至连英国流体力学家斯托克斯(GeorgeGabrielStokes,1819-1903)爵士也对此提3出质疑，怀疑在静止水面上能存在不变形的行波。他们的怀疑的问题主要有：一个完整的波动为什么会全部在水面上，而不是一部分在水面上，一部分在水面下；波在传播的过程中，为什么波幅不会衰减；波的运动速度也与他们的研究结果不符。这一争论延续到19世纪70年代才初步得到解决。H.E.巴津(Bazin,H.E.)和英国科学家瑞利(JohnWilliamStrutRayleigh，1842-1919)对孤立波进行了一系列的研究，证明了Russell的工作是正确的。Russell与艾里，斯托克斯的争论，最终于1895年由数学家D.J.科尔特弗(Korteweg，D.J.)和他的学生G.德.弗里斯(Vires,G.de)所解决。他们在小振幅与长波的假定下，从流体动力学导出了关于孤立波的方程(后人称它为KdV方程)。这一方程的行波解，在波长趋于无限的情况下，正是Russell所发现的孤立波。KdV方程的提出，从理论上阐明了孤立波的存在，给这场争论划上了句号。从Russell的发现到KdV方程的提出，大约经历了60年时间,孤立波才为学术界普遍接受。数学家Korteweg,D.J和G.deVires从流体动力学中导出了关于孤立波的方程(后人称它为KdV方程)033xuxuutu，第二项空间采用中心差分蛙跳大部分时间(1)在数学物理方程中，可解得其精确解—行波解[9]，如下DBtAxhCtxu2sec3),(，(2)其中C和D是常数；CA21,ACB。下面从物理学和数学角度理解孤立波的特点,物理学：孤立波是物质非线性效应的一种特殊产物。KdV方程中含有非线性项xuu。数学：某些非线性偏微分方程的一类稳定的，能量有限的不弥散解。孤立波可谓是孤立，如下1.在传播过程中保持波形和速度不变,像是“透明地”穿过对方。2.碰撞时两个孤立波重叠在一起,其高度低于碰撞前孤立波高度较高的一个(这表明在非线性过程中，不存在线性叠加原理)两个孤立波碰撞时互相穿透且维持原来的波形和速度；孤立波的波幅愈高，其传播速度愈快，如下图(2)。4图2.两孤立波的传播的交会情形，(a)1tt；(b)12ttt；(a)23ttt；(a)43ttt。3.碰撞后孤立波的轨道与碰撞前有些偏离(即发生了相移)。在海洋中也存在孤立波，如下1.涌浪向缓慢倾斜的海岸传播时，随着涌浪接近海岸，波峰变尖，到水深很浅的地方，波峰成了隆起的状态，被很长而又很平的波谷隔开，一个个波峰好像成了孤立波，因而可引入孤立波理论[13]。2.在分层流体中也发现了孤立波的现象，这种内孤立波在两层流体之间沿某一方向传播。海洋中有大振幅内潮时，其可能会演变出的强孤立内波会干扰海洋工程作业，对建成的石油钻井平台和海底油气管道构成严重的威胁；内潮导致的水体扰动还会对海上舰船的行驶产生不利影响[14]。此外，大气中也存在孤立波。孤立子理论应用如此广泛，可见研究孤立波的差分格式及其数值模拟的误差对生产具有非常重要的意义。2.有限差分格式有限差分方法是解微分方程和积分微分方程数值解的方法。基本思想是把连续的定解区域用有限个离散点构成的网格来代替，这些离散点称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解。总的来说，有限差分法的基本概念是用差商代替微商。有限差分法的主要内容包括:如何根据问题的特点将定解区域作网格剖分；如何把原微分方程离散化为差分方程组以及如何解此代数方程组。此外为了保证计算过程的可行和计算结果的正确，还需从理论上分析差分方程组的性态，包括5解的唯一性，存在性和差分格式的相容性，收敛性和稳定性。对于一个微分方程建立的各种差分格式，为了有实用意义，一个基本要求是它们能够任意逼近微分方程，这就是相容性要求。另外，一个差分格式是否有用，最终要看差分方程的精确解能否任意逼近微分方程的解，这就是收敛性的概念。此外，还有一个重要的概念必须考虑，即差分格式的稳定性。因为差分格式的计算过程是逐层推进的，在计算第n＋1层的近似值时要用到第n层的近似值，直到与初始值有关。前面各层若有舍入误差，必然影响到后面各层的值，如果误差的影响越来越大，以致差分格式的精确解的面貌完全被掩盖，这种格式是不稳定的，相反如果误差的传播是可以控制的，就认为格式是稳定的。只有在这种情形，差分格式在实际计算中的近似解才可能任意逼近差分方程的精确解。关于差分格式的构造一般有以下3种方法。最常用的方法是数值微分法，比如用差商代替微商等。另一方法叫积分插值法，因为在实际问题中得出的微分方程常常反映物理上的某种守恒原理，一般可以通过积分形式来表示。此外还可以用待定系数法构造一些精度较高的差分格式。下面给出KdV方程的三种有限差分格式。我们假设在网格点上，对于njtxu,有近似值}{:nmnuu，其中xmxm，Nm,,2,1,0和tntn，0n，其中x是空间步长，t是时间步长。对空间步长x和时间步长，有xtr，这里r是网比。为了描述有限差分格式，我们使用了中心差分算子0x[15]，定义如下nmnmnmxuuu110。(3)2.1蛙跳格式对于时间微分项tu采用中心差分格式(蛙跳格式)，有tuutunmnmnm211；(4)时间微分，第m个网格点不变；时间步长对于非线性项中xu也采用中心差分格式，有xuuxunmnmnm211，(5)6根据中心差分算子0x，(5)可化为xuxunmxnm20；(6)对于3阶微分项33xu，先对2阶微分采用中心差分，再对二阶微分进行一阶中心差分，有xxuxuxuxxunmnmnmnm21221222233，(7)将nmu1和nmu1泰勒展开，得3332221!31!21xxuxxuxxuuunmnmnmnmnm444!41xxunm，(8)3332221!31!21xxuxxuxxuuunmnmnmnmnm444!41xxunm，(9)将(8)与(9)相加得到二阶导数的差分格式211222xuuuxunmnmnmnm，(10)根据中心差分算子0x，(5)可化为31103322xuuuxunmnmnmxnm；(11)所以，将(4)，(6)，(11)代入到方程(1)中，得到方程(1)的近似差分方程，如下022223110011xuuuxuutuunmnmnmxnmxnmnmnm，(12)其中nmnmnmnmuuuu1131，并化简得023110011xuuuxuutuunmnmnmxnmxnmnmnm，(13)最终差分格式，由matlab编写将网比xtr代入上式，得7nmnmnmnmnmnmnmuuuuuruu11111131nmnmnmnmuuuuxr2112222，(14)其稳定性条件[15]是14max2xur。(15)2.2跳点格式对于时间微分项tu采用向前差分，有tuutunmnmnm1；(16)至于对空间变量的微分，跳点格式是根据nm,的和的不同，有不同的差分格式，其叙述如下。当偶数nm时，对于非线性项中空间微分xu，在tn时刻采用中心差分格式，即(5)式。对于3阶微分项33xu，在tn时刻先对2阶微分采用中心差分，再对二阶微分进行一阶中心差分，即(11)式。将(5)式，(11)式，(16)式代入到方程(1)中，得0222311001xuuuxuutuunmnmnmxnmxnmnmnm，(17)其中nmnmnmuuu1121，上式可化简，得0222311001xuuuxftuunmnmnmxnmxnmnm(18)其中22uf，化简后，得nmnmnmnmnmxnmnmuuuuxrfruu2112201222121。(19)当奇数nm时，对空间的微分项都是采用的隐式格式，即对空间的微分都是取在tn1时刻，表述如下。对于非线性项中的xu，仍然是在tn1时刻采用中心差分格式，得xuuxunmnmnm