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

1、1目 录 微分方程数值微分方程数值实验实验.21. 孤立波简介 .22. 有限差分格式 .42.1 蛙跳格式 .52.2 跳点格式 .72.3 分裂格式 .83. 实例 .94. 附录：MATLAB源码 .122微分方程数值实验微分方程数值实验KDV 方程是历史上最著名的方程之一，其起源是水波中的孤立波。1. 孤立波简介我们自己先设想一下：在一条窄河道中，当迅速前进的船突然停下时，在船头形成的一个孤立的水波迅速离开船头,继续向前运动，如下图(1)所示。图1.孤立波的运动图孤立波是由是苏格兰一位优秀的造船工程师 Russell(John Scott Russell 1808-1882)在 183

2、4 年研究船舶在运动中所受到的阻力时发现的。在一次学术报告中,他是这样描述的：“我把注意力集中在船舶给予流体的运动上，立刻就观察到一个非同寻常而又非常绚丽的现象，它是如此之重要，以致我将首先详细描述它所表现出来的外貌。当我正在观察一只高速运动的船舶，让它突然停止时，在船舶周围所形成的小波浪中，一个紊乱的扰动现象吸引了我的注意。在船身长度的中部附近，许多水聚集在一起，形成一个廓线很清楚的水堆，最后还出现一尖峰，并以相当高的速度开始向前运动，到船头后,继续保持它的形状不变，在静止流体的表面上，完全孤立地向前运动，成为一孤立行进波,直到河道的转弯处才开始消失掉。 ”纵观力学和物理学的发展史不难看到，

3、每当开始引入一种新思想或新概念的时候，往往会受到怀疑和非难，并常会引起激烈的争论，孤立波的命运也是如此。3当时科学界的权威们对 Russell 的这些结果，一开始时就表示了怀疑和反对。甚至连英国流体力学家斯托克斯(George Gabriel Stokes,1819-1903)爵士也对此提出质疑，怀疑在静止水面上能存在不变形的行波。他们的怀疑的问题主要有：一个完整的波动为什么会全部在水面上，而不是一部分在水面上，一部分在水面下；波在传播的过程中，为什么波幅不会衰减；波的运动速度也与他们的研究结果不符。这一争论延续到 19 世纪 70 年代才初步得到解决。H.E.巴津(Bazin,H.E.)和英

4、国科学家瑞利(John William Strut Rayleigh，1842-1919)对孤立波进行了一系列的研究，证明了 Russell 的工作是正确的。Russell 与艾里，斯托克斯的争论，最终于 1895 年由数学家 D.J.科尔特弗(Korteweg ，D.J.)和他的学生 G.德.弗里斯(Vires ,G.de)所解决。他们在小振幅与长波的假定下，从流体动力学导出了关于孤立波的方程(后人称它为 KdV 方程)。这一方程的行波解，在波长趋于无限的情况下，正是 Russell 所发现的孤立波。KdV 方程的提出，从理论上阐明了孤立波的存在，给这场争论划上了句号。从 Russell 的

5、发现到 KdV 方程的提出，大约经历了 60 年时间,孤立波才为学术界普遍接受。数学家 Korteweg ,D.J 和 G.de Vires 从流体动力学中导出了关于孤立波的方程(后人称它为 KdV 方程) ， 033xuxuutu第二项空间采用中心差分第二项空间采用中心差分 蛙跳大部分时间蛙跳大部分时间 (1)在数学物理方程中，可解得其精确解行波解9，如下， (2)DBtAxhCtxu2sec3),(其中和是常数；,。CDCA21ACB下面从物理学和数学角度理解孤立波的特点,物理学：孤立波是物质非线性效应的一种特殊产物。KdV方程中含有非线性项。xuu数学：某些非线性偏微分方程的一类稳定的，

6、能量有限的不弥散解。孤立波可谓是孤立，如下1. 在传播过程中保持波形和速度不变,像是“透明地”穿过对方。2. 碰撞时两个孤立波重叠在一起,其高度低于碰撞前孤立波高度较高的一个(这表明在非线性过程中，不存在线性叠加原理)两个孤立波碰撞时互相穿透且维持原来的波形和速度；孤立波的波幅愈高，其传播速度愈快，如下图(2)。4 图 2.两孤立波的传播的交会情形，(a)；(b)； (a)； (a)1tt 12ttt23ttt。43ttt3. 碰撞后孤立波的轨道与碰撞前有些偏离(即发生了相移)。在海洋中也存在孤立波，如下1. 涌浪向缓慢倾斜的海岸传播时，随着涌浪接近海岸，波峰变尖，到水深很浅的地方，波峰成了隆

7、起的状态，被很长而又很平的波谷隔开，一个个波峰好像成了孤立波，因而可引入孤立波理论13。2. 在分层流体中也发现了孤立波的现象，这种内孤立波在两层流体之间沿某一方向传播。海洋中有大振幅内潮时，其可能会演变出的强孤立内波会干扰海洋工程作业，对建成的石油钻井平台和海底油气管道构成严重的威胁；内潮导致的水体扰动还会对海上舰船的行驶产生不利影响14。此外，大气中也存在孤立波。孤立子理论应用如此广泛，可见研究孤立波的差分格式及其数值模拟的误差对生产具有非常重要的意义。2. 有限差分格式有限差分方法是解微分方程和积分微分方程数值解的方法。基本思想是把连续的定解区域用有限个离散点构成的网格来代替，这些离散点

8、称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解。总的来说，有限差分法的基本概念是用差商代替微商。有限差分法的主要内容包括:如何根据问题的特点将定解区域作网格剖分；5如何把原微分方程离散化为差分方程组以及如何解此代数方程组。此外为了保证计算过程的可行和计算结果的正确，还需从理论上分析差分方程组的性态，包括解的唯一性，存在性和

9、差分格式的相容性，收敛性和稳定性。对于一个微分方程建立的各种差分格式，为了有实用意义，一个基本要求是它们能够任意逼近微分方程，这就是相容性要求。另外，一个差分格式是否有用，最终要看差分方程的精确解能否任意逼近微分方程的解，这就是收敛性的概念。此外，还有一个重要的概念必须考虑，即差分格式的稳定性。因为差分格式的计算过程是逐层推进的，在计算第 n1 层的近似值时要用到第 n 层的近似值，直到与初始值有关。前面各层若有舍入误差，必然影响到后面各层的值，如果误差的影响越来越大，以致差分格式的精确解的面貌完全被掩盖，这种格式是不稳定的，相反如果误差的传播是可以控制的，就认为格式是稳定的。只有在这种情形，

10、差分格式在实际计算中的近似解才可能任意逼近差分方程的精确解。关于差分格式的构造一般有以下 3 种方法。最常用的方法是数值微分法，比如用差商代替微商等。另一方法叫积分插值法，因为在实际问题中得出的微分方程常常反映物理上的某种守恒原理，一般可以通过积分形式来表示。此外还可以用待定系数法构造一些精度较高的差分格式。下面给出 KdV 方程的三种有限差分格式。我们假设在网格点上，对于有近似值，其中njtxu,:nmnuu，和，xmxmNm, 2 , 1 , 0tntn0n其中是空间步长，是时间步长。对空间步长和时间步长，有xtx，xtr这里 是网比。为了描述有限差分格式，我们使用了中心差分算子15，r0

11、 x定义如下。 (3)nmnmnmxuuu1102.1 蛙跳格式对于时间微分项采用中心差分格式(蛙跳格式)，有tu； tuutunmnmnm211(4)时间微分，第m个网格点不变；时间步长对于非线性项中也采用中心差分格式，有xu6， (5)xuuxunmnmnm211根据中心差分算子，(5)可化为0 x； xuxunmxnm20(6)对于3阶微分项，先对2阶微分采用中心差分，再对二阶微分进行一阶33xu中心差分，有， (7)xxuxuxuxxunmnmnmnm21221222233将和泰勒展开，得nmu1nmu13332221! 31! 21xxuxxuxxuuunmnmnmnmnm， (8)

12、444! 41xxunm 3332221! 31! 21xxuxxuxxuuunmnmnmnmnm， (9)444! 41xxunm将(8)与(9)相加得到二阶导数的差分格式， (10)211222xuuuxunmnmnmnm根据中心差分算子，(5)可化为0 x； (11)31103322xuuuxunmnmnmxnm 所以，将(4)，(6)，(11)代入到方程(1)中，得到方程(1)的近似差分方程，如下， (12)022223110011xuuuxuutuunmnmnmxnmxnmnmnm其中，并化简得nmnmnmnmuuuu11317， (13)最023110011xuuuxuutuunm

13、nmnmxnmxnmnmnm终差分格式，由matlab编写将网比代入上式，得xtr nmnmnmnmnmnmnmuuuuuruu11111131， (14)nmnmnmnmuuuuxr2112222其稳定性条件15是。 (15)14max2xur2.2 跳点格式对于时间微分项采用向前差分，有tu； tuutunmnmnm1(16)至于对空间变量的微分，跳点格式是根据的和的不同，有不同的差分nm,格式，其叙述如下。当时，对于非线性项中空间微分，在时刻采用中心差偶数 nmxutn分格式， 即(5)式。对于3阶微分项，在时刻先对2阶微分采用中心差分，再对二阶微分33xutn进行一阶中心差分，即(11

14、)式。 将(5)式，(11)式，(16)式代入到方程(1)中，得， (17)0222311001xuuuxuutuunmnmnmxnmxnmnmnm其中，上式可化简，得nmnmnmuuu1121 (18)0222311001xuuuxftuunmnmnmxnmxnmnm其中，化简后，得22uf 。 (19)nmnmnmnmnmxnmnmuuuuxrfruu2112201222121当时，对空间的微分项都是采用的隐式格式，即对空间的微分奇数 nm8都是取在时刻，表述如下。tn1对于非线性项中的，仍然是在时刻采用中心差分格式，得xutn1； (20)xuuxunmnmnm211111对于3阶微分项

15、中的，在时刻先对2阶微分采用中心差分，再进33xutn1行一阶中心差分，得。 (21)311111013322xuuuxunmnmnmxnm 将(16)式，(20)式，(21)式代入到方程(1)中得， (22)022231111101011xuuuxuutuunmnmnmxnmxnmnmnm其中，上式可化简，得1111121nmnmnmuuu， (23)02223111110101xuuuxftuunmnmnmxnmxnmnm其中，化简后，得22uf 。 (24)121111122101222121nmnmnmnmnmxnmnmuuuuxrfruu跳点格式的稳定性条件15是 。 (25)122

16、xur2.3 分裂格式将非线性项和频散项分开后独立计算。对流方程可用迎风格0 xuutu式或 Lax 等耗散格式，频散方程可取一般格式如中心差分格式，033xutuCN 格式等。例如对流方程取 Lax 格式，频散方程取时间前差空间中心差分格式，则得，nmnmnmnmnmffruuu11112121，nmnmnmnmnmuuuxruu1121221精度为，稳定性条件用(15)式，当系数较大时，频散项起重)(2xt33xu9要的作用，如该项用显式格式，要求，而用全隐式增加计算量，为此3xt可采用半隐式格式，得， (26) 022321111111201xuuuuuuxftuunmnmnmnmnmn

17、mnmxnmnm化简得 nmxnmnmfruu01121， (27) nmnmnmnmnmnmuuuuuuxr21111111222这是一个三层隐式格式，稳定性条件15为 。(28)1ur3. 实例对于方程(1)，系数和参数的选取：，;0 . 141084. 43 . 0C0 . 6D时间步长，空间步长。根据系数之间的关系，得0001. 0t02. 0 x ，448.1221CA。3.7345ACB根据行波解的特性，此时孤立波向前传播的速度可以求出 。3000. 0ABv用Matlab对三种差分格式进行编程计算(程序见附录)，至于初始值和边界值，是根据精确解(即行波解)给出的，即 (29)。D

18、BtAxhCtxuuDBtAxhCtxuuDBtAxhCtxuuDBtAxhCtxuuDBtAxhCtxuuDAxhCtxuunMnMnMnMnMnMnnnnnnmmmmmm21211121102001211200sec3,sec3,sec3,sec3,sec3,sec3,计算后作出时，精确解以及三种差分格式的解的图，如下1 . 0t10 (a) (b) (c) (d)图3.各种格式的解得图：(a)精确解；(b)蛙跳格式；(c)跳点格式；(d)分裂格式；上面的(c)图明显不同于另外三幅图，是因为模拟的结果中有负值存在，如表1中所示，而精确解(2)是一平方表达式，不可能小于零，蛙跳格式和分裂格式

19、在其他时刻某些位置也有负值(如附表2)，只是在时所取的点处没有负1 . 0t值。表1给出了时的振幅，表2给出了的位置不同时刻的振幅。1 . 0t2x表1. 时不同位置的振幅1 . 0tx精确解蛙跳跳点分裂01.0481e-0051.0481e-0051.0481e-0051.0481e-0050.10.00152230.0025050.00212250.00212980.20.196560.200540.198890.199420.30.325620.320950.323080.322590.40.00276450.00285950.00287410.00286930.51.9046e-005

20、2.0717e-0052.0292e-0052.0246e-0050.61.3102e-0071.4048e-0071.3921e-0071.3733e-0070.79.0125e-0109.4558e-0109.6225e-0109.4272e-0100.86.1996e-0126.5007e-0126.6283e-0126.4889e-012110.94.2646e-0144.4718e-0144.5597e-0144.4638e-01412.9336e-0163.0761e-0163.1366e-0163.0706e-0161.12.018e-0182.116e-0182.1576e-0

21、182.1122e-0181.21.3882e-0201.4556e-0201.4888e-0201.453e-0201.39.549e-0231.0013e-0228.893e-0239.9949e-0231.46.5687e-0256.8878e-0253.7904e-0226.8754e-0251.54.5185e-0274.7381e-0275.7569e-0174.7295e-0271.63.1083e-0293.2593e-0291.2507e-0133.2534e-0291.72.1381e-0312.242e-031-1.0754e-0092.238e-0311.81.4708

22、e-0331.5423e-033-2.592e-0071.5395e-0331.91.0118e-0351.0609e-0352.6574e-0061.059e-03526.9597e-0386.9597e-0386.9597e-0386.9597e-038表2.的位置不同时刻的振幅2xt精确解蛙跳跳点分裂01.39e-0161.39e-0161.39e-0161.39e-0160.011.4978e-0161.505e-0161.5049e-0161.5046e-0160.021.614e-0161.6294e-0161.6293e-0161.6287e-0160.031.7392e-016

23、1.7641e-0161.764e-0161.7631e-0160.041.874e-0161.9099e-0161.9098e-0161.9085e-0160.052.0194e-0162.0678e-0162.0678e-0162.0659e-0160.062.176e-0162.2388e-0162.2389e-0162.2363e-0160.072.3447e-0162.4239e-0162.4247e-0162.4208e-0160.082.5265e-0162.6243e-0162.6282e-0162.6205e-0160.092.7225e-0162.8412e-0162.85

24、69e-0162.8366e-0160.12.9336e-0163.0761e-0163.1366e-0163.0706e-0160.113.1611e-0163.3304e-0163.5613e-0163.3239e-0160.123.4062e-0163.6058e-0164.4844e-0163.598e-0160.133.6704e-0163.9037e-0167.2436e-0163.8947e-0160.143.955e-0164.2247e-0161.6914e-0154.2134e-0160.154.2617e-0164.5738e-0165.2713e-0154.5196e-

25、0160.164.5922e-0165.3788e-0161.8691e-0144.3975e-0160.174.9484e-0161.5447e-0156.8746e-0143.4968e-0170.185.3321e-0161.5021e-0142.528e-013-2.7561e-0150.195.7456e-0161.536e-0139.1822e-013-1.3001e 0140.26.1912e-0161.2862e-0123.2994e-0126.8679e-015时，三种差分格式的误差曲线沿空间的变化趋势大体相同，如图 4 所示。1 . 0t12图 4. 时，三种差分格式的误差

26、曲线沿空间的变化趋势1 . 0t由图 4 得，在的范围内误差很小，几乎紧贴零线，表明差分格式是收敛的；1x在内，误差的变化波动较大。在内三种差分格式的误差变化如 1 , 0 x 1 , 0 x下图 5。图 5. 内三种差分格式的误差变化图 1 , 0 x4. 附录：附录：Matlab 源码源码蛙跳格式程序：蛙跳格式程序：13运行首先工作取决运行首先工作取决DebugDebug runrun工作空间出现所有参数工作空间出现所有参数右击查看参数信息右击查看参数信息%蛙跳格式clear all，clc 清空%系数参数赋值 a=1.0;b=4.84e-4;C=0.3;D=-6;%取步长tb=0.000

27、1; %时间步长xb=0.02; %空间步长倍数可改变调试 A=0.5*(a*C/b)0.5;B=a*A*C; k1=-1/3*a*tb/xb;k2=-tb*b/(xb)3; mmax=201;nmax=2001; 精确解空间网格点数 for i=1:mmax %利用精确解给出初值，n=0和n=1时 x(i)=(i-1)*xb; u(i,1)=3*C*(sech(A*x(i)+D)2; u(i,2)=3*C*(sech(A*x(i)-B*tb+D)2; uju(i,1:2)=0;end for j=1:nmax %利用精确解给出边值，m=0,m=1,m=199,m=200 t(j)=(j-1)

28、*tb; u(1,j)=3*C*(sech(A*x(1)-B*t(j)+D)2; u(2,j)=3*C*(sech(A*x(2)-B*t(j)+D)2; u(mmax-1,j)=3*C*(sech(A*x(mmax-1)-B*t(j)+D)2; u(mmax,j)=3*C*(sech(A*x(mmax)-B*t(j)+D)2; uju(1:2,j)=0; uju(mmax-1:mmax,j)=0;end for m=1:mmax %计算精确解 x(m)=(m-1)*xb; for n=1:nmax t(n)=(n-1)*tb; uj(m,n)=3*C*(sech(A*x(m)-B*t(n)+D

29、)2;14 endend for n=3:nmax %蛙跳格式计算 for m=3:mmax-2 u(m,n)=u(m,n-2)+k1*(u(m+1,n-1)+u(m,n-1)+u(m-1,n-1)*(u(m+1,n-1)-u(m-1,n-1)+k2*(u(m+2,n-1)-2*u(m+1,n-1)+2*u(m-1,n-1)-u(m-2,n-1); endend跳点格式程序跳点格式程序：%跳点格式clear all，clc%系数参数赋值a=1.0;b=4.84e-4;C=0.3;D=-6; tb=0.0001; %时间步长xb=0.02; %空间步长 A=0.5*sqrt(a*C/b);B=a

30、*A*C; k1=-1/4*a*tb/xb;k2=-tb*b/(2*(xb)3); mmax=201;nmax=2001; for i=1:mmax %利用精确解给出初值条件,t=0时 x(i)=(i-1)*xb; u(i,1)=3*C*(sech(A*x(i)+D)2; uju(i,1)=0;end for j=1:nmax %利用精确解给出边值条件 t(j)=(j-1)*tb; u(1,j)=3*C*(sech(A*x(1)-B*t(j)+D)2; u(2,j)=3*C*(sech(A*x(2)-B*t(j)+D)2; u(mmax-1,j)=3*C*(sech(A*x(mmax-1)-B

31、*t(j)+D)2; u(mmax,j)=3*C*(sech(A*x(mmax)-B*t(j)+D)2; uju(1:2,j)=0; uju(mmax-1:mmax,j)=0;end for m=3:mmax-2 %迭代中未知点的初始赋值 for n=2:nmax u(m,n)=0.5; uu(m,n)=0.6;15 endend for n=2:nmax %跳点格式计算 for m=3:mmax-2 %计算显式值 if (mod(m+n),2)=0) u(m,n)=u(m,n-1)+k1*(u(m+1,n-1)2-(u(m-1,n-1)2)+k2*(u(m+2,n-1)-2*u(m+1,n-1)+2*u(m-1,n-1)-u(m-2,n-1); uu(m,n)=u(m,n); end end j=0; while (1) %迭代计算隐式值 j=j+1; a=0; for m=3:mmax-2 if (mod(m+n),2)=1) uu(m,n)=u(m,n-1)+k1*(u(m+1,n)2-(u(m-1,n)2)+k2*(u(m+2,n)-2*u(m+1,n)+2*u(m-1,n)-u(m-2,n); end end for m=3:mmax-2 %计算总的误差 i