《海洋數(shù)值模擬》實(shí)驗(yàn)2指導(dǎo)書解析

微分方程數(shù)值實(shí)驗(yàn)錯(cuò)誤！未定義書簽。微分方程數(shù)值實(shí)驗(yàn)TOC\o"1-5"\h\z\o"CurrentDocument"孤立波簡(jiǎn)介 2\o"CurrentDocument"有限差分格式 4\o"CurrentDocument"2.1蛙跳格式 5\o"CurrentDocument"2.2跳點(diǎn)格式 7\o"CurrentDocument"2.3分裂格式 8實(shí)例 9附錄：Matlab源碼 12微分方程數(shù)值實(shí)驗(yàn)KDV方程是歷史上最著名的方程之一，其起源是水波中的孤立波。孤立波簡(jiǎn)介我們自己先設(shè)想一下：在一條窄河道中，當(dāng)迅速前進(jìn)的船突然停下時(shí)，在船頭形成的一個(gè)孤立的水波迅速離開船頭,繼續(xù)向前運(yùn)動(dòng)，如下圖(1)所示。孤立波是由是蘇格蘭一位優(yōu)秀的造船工程師 Russell(JohnScottRussell1808-1882)在1834年研究船舶在運(yùn)動(dòng)中所受到的阻力時(shí)發(fā)現(xiàn)的。在一次學(xué)術(shù)報(bào)告中，他是這樣描述的：“我把注意力集中在船舶給予流體的運(yùn)動(dòng)上，立刻就觀察到一個(gè)非同尋常而又非常絢麗的現(xiàn)象，它是如此之重要，以致我將首先詳細(xì)描述它所表現(xiàn)出來(lái)的外貌。當(dāng)我正在觀察一只高速運(yùn)動(dòng)的船舶，讓它突然停止時(shí)，在船舶周圍所形成的小波浪中，一個(gè)紊亂的擾動(dòng)現(xiàn)象吸引了我的注意。在船身長(zhǎng)度的中部附近，許多水聚集在一起，形成一個(gè)廓線很清楚的水堆，最后還出現(xiàn)一尖峰，并以相當(dāng)高的速度開始向前運(yùn)動(dòng)，到船頭后，繼續(xù)保持它的形狀不變，在靜止流體的表面上，完全孤立地向前運(yùn)動(dòng)，成為一孤立行進(jìn)波，直到河道的轉(zhuǎn)彎處才開始消失掉?！笨v觀力學(xué)和物理學(xué)的發(fā)展史不難看到，每當(dāng)開始引入一種新思想或新概念的時(shí)候，往往會(huì)受到懷疑和非難，并常會(huì)引起激烈的爭(zhēng)論，孤立波的命運(yùn)也是如此。當(dāng)時(shí)科學(xué)界的權(quán)威們對(duì)Russell的這些結(jié)果，一開始時(shí)就表示了懷疑和反對(duì)。甚至連英國(guó)流體力學(xué)家斯托克斯(GeorgeGabrielStokes,1819-1903)爵士也對(duì)此提出質(zhì)疑，懷疑在靜止水面上能存在不變形的行波。他們的懷疑的問(wèn)題主要有：一個(gè)完整的波動(dòng)為什么會(huì)全部在水面上，而不是一部分在水面上，一部分在水面下;波在傳播的過(guò)程中，為什么波幅不會(huì)衰減；波的運(yùn)動(dòng)速度也與他們的研究結(jié)果不符。這一爭(zhēng)論延續(xù)到19世紀(jì)70年代才初步得到解決。H.E.巴津(Bazin,H.E.)和英國(guó)科學(xué)家瑞利(JohnWilliamStrutRayleigh，1842-1919)對(duì)孤立波進(jìn)行了一系列的研究，證明了Russell的工作是正確的。Russell與艾里，斯托克斯的爭(zhēng)論，最終于1895年由數(shù)學(xué)家D.J.科爾特弗(Korteweg，D.J.)和他的學(xué)生G.德.弗里斯(Vires,G.de)所解決。他們?cè)谛≌穹c長(zhǎng)波的假定下，從流體動(dòng)力學(xué)導(dǎo)出了關(guān)于孤立波的方程(后人稱它為KdV方程)。這一方程的行波解，在波長(zhǎng)趨于無(wú)限的情況下，正是Russell所發(fā)現(xiàn)的孤立波。KdV方程的提出，從理論上闡明了孤立波的存在，給這場(chǎng)爭(zhēng)論劃上了句號(hào)。從Russell的發(fā)現(xiàn)到KdV方程的提出，大約經(jīng)歷了60年時(shí)間，孤立波才為學(xué)術(shù)界普遍接受。數(shù)學(xué)家Korteweg,D.J和GdeVires從流體動(dòng)力學(xué)中導(dǎo)出了關(guān)于孤立波的方程(后人稱它為KdV方程)du du d3u八TOC\o"1-5"\h\z——+以u(píng)——+p =0，dt dx dx3第二項(xiàng)空間采用中心差分蛙跳大部分時(shí)間(1)在數(shù)學(xué)物理方程中，可解得其精確解一行波解［9］,如下\o"CurrentDocument"u(x,t)=3Csech2"x-Bt+D］， (2)其中C和D是常數(shù)；A=上,B=aACo2\：P下面從物理學(xué)和數(shù)學(xué)角度理解孤立波的特點(diǎn)，物理學(xué)：孤立波是物質(zhì)非線性效應(yīng)的一種特殊產(chǎn)物。KdV方程中含有非線性項(xiàng)duau—。dx數(shù)學(xué)：某些非線性偏微分方程的一類穩(wěn)定的，能量有限的不彌散解。孤立波可謂是孤立，如下在傳播過(guò)程中保持波形和速度不變，像是“透明地”穿過(guò)對(duì)方。碰撞時(shí)兩個(gè)孤立波重疊在一起,其高度低于碰撞前孤立波高度較高的一個(gè)(這表明在非線性過(guò)程中，不存在線性疊加原理)兩個(gè)孤立波碰撞時(shí)互相穿透且維持原來(lái)的波形和速度；孤立波的波幅愈高，其傳播速度愈快，如下圖(2)。圖2.兩孤立波的傳播的交會(huì)情形，(a)t=t;(b)t=t>t;(a)t=t>t;(a)t=t>t。1 2 1 3 2 3 4碰撞后孤立波的軌道與碰撞前有些偏離(即發(fā)生了相移)。在海洋中也存在孤立波，如下涌浪向緩慢傾斜的海岸傳播時(shí)，隨著涌浪接近海岸，波峰變尖，到水深很淺的地方，波峰成了隆起的狀態(tài)，被很長(zhǎng)而又很平的波谷隔開，一個(gè)個(gè)波峰好像成了孤立波，因而可引入孤立波理論【⑶。在分層流體中也發(fā)現(xiàn)了孤立波的現(xiàn)象，這種內(nèi)孤立波在兩層流體之間沿某一方向傳播。海洋中有大振幅內(nèi)潮時(shí)，其可能會(huì)演變出的強(qiáng)孤立內(nèi)波會(huì)干擾海洋工程作業(yè)，對(duì)建成的石油鉆井平臺(tái)和海底油氣管道構(gòu)成嚴(yán)重的威脅；內(nèi)潮導(dǎo)致的水體擾動(dòng)還會(huì)對(duì)海上艦船的行駛產(chǎn)生不利影響［14］。此外，大氣中也存在孤立波。孤立子理論應(yīng)用如此廣泛，可見(jiàn)研究孤立波的差分格式及其數(shù)值模擬的誤差對(duì)生產(chǎn)具有非常重要的意義。有限差分格式有限差分方法是解微分方程和積分微分方程數(shù)值解的方法?；舅枷胧前堰B續(xù)的定解區(qū)域用有限個(gè)離散點(diǎn)構(gòu)成的網(wǎng)格來(lái)代替，這些離散點(diǎn)稱作網(wǎng)格的節(jié)點(diǎn)；把連續(xù)定解區(qū)域上的連續(xù)變量的函數(shù)用在網(wǎng)格上定義的離散變量函數(shù)來(lái)近似；把原方程和定解條件中的微商用差商來(lái)近似，積分用積分和來(lái)近似，于是原微分方程和定解條件就近似地代之以代數(shù)方程組，即有限差分方程組，解此方程組就可以得到原問(wèn)題在離散點(diǎn)上的近似解。然后再利用插值方法便可以從離散解得到定解問(wèn)題在整個(gè)區(qū)域上的近似解?？偟膩?lái)說(shuō)，有限差分法的基本概念是用差商代替微商。有限差分法的主要內(nèi)容包括:如何根據(jù)問(wèn)題的特點(diǎn)將定解區(qū)域作網(wǎng)格剖分；如何把原微分方程離散化為差分方程組以及如何解此代數(shù)方程組。此外為了保證計(jì)算過(guò)程的可行和計(jì)算結(jié)果的正確，還需從理論上分析差分方程組的性態(tài)，包括

解的唯一性，存在性和差分格式的相容性，收斂性和穩(wěn)定性。對(duì)于一個(gè)微分方程建立的各種差分格式，為了有實(shí)用意義，一個(gè)基本要求是它們能夠任意逼近微分方程，這就是相容性要求。另外，一個(gè)差分格式是否有用，最終要看差分方程的精確解能否任意逼近微分方程的解，這就是收斂性的概念。此外，還有一個(gè)重要的概念必須考慮，即差分格式的穩(wěn)定性。因?yàn)椴罘指袷降挠?jì)算過(guò)程是逐層推進(jìn)的，在計(jì)算第n+1層的近似值時(shí)要用到第n層的近似值，直到與初始值有關(guān)。前面各層若有舍入誤差，必然影響到后面各層的值，如果誤差的影響越來(lái)越大，以致差分格式的精確解的面貌完全被掩蓋，這種格式是不穩(wěn)定的，相反如果誤差的傳播是可以控制的，就認(rèn)為格式是穩(wěn)定的。只有在這種情形，差分格式在實(shí)際計(jì)算中的近似解才可能任意逼近差分方程的精確解。關(guān)于差分格式的構(gòu)造一般有以下3種方法。最常用的方法是數(shù)值微分法，比如用差商代替微商等。另一方法叫積分插值法，因?yàn)樵趯?shí)際問(wèn)題中得出的微分方程常常反映物理上的某種守恒原理，一般可以通過(guò)積分形式來(lái)表示。此外還可以用待定系數(shù)法構(gòu)造一些精度較高的差分格式。下面給出KdV方程的三種有限差分格式。我們假設(shè)在網(wǎng)格點(diǎn)上，對(duì)于uX,tn)<■有近似值un:=｛加｝，其中x=mAx，m=0,1,2,,N和t=nAt，n>0，m n其中A是空間步長(zhǎng)，At是時(shí)間步長(zhǎng)。對(duì)空間步長(zhǎng)A和時(shí)間步長(zhǎng)，有_AtAx(3)(4)(5)這里r是網(wǎng)比。為了描述有限差分格式，我們使用了中心差分算子5o[15i定義如下(3)(4)(5)50Un=Un-Un。2.1蛙跳格式"’心對(duì)于時(shí)間微分項(xiàng)竺采用中心差分格式(蛙跳格式)，有所‘8")nUn+1—Un-13)?m2Atm；m時(shí)間微分，第m個(gè)網(wǎng)格點(diǎn)不變；時(shí)間步長(zhǎng)對(duì)于非線性項(xiàng)中坐也采用中心差分格式，有8x8u'n Un—un8x) 2Axm

根據(jù)中心差分算子50,(5)可化為了TOC\o"1-5"\h\z即丫 50un 心I—— -xm ； (6)"辦) 2Axm對(duì)于3階微分項(xiàng)來(lái)，先對(duì)2階微分采用中心差分，再對(duì)二階微分進(jìn)行一階中OX3心差分，有傳2u丫 仔2u丫—— ———"辦2) "8X2)2Ax將Un和Un泰勒展開，得m+1 m-1Unm+1_ (OUnm+1_ (Ou、=um*［云)mnAx+—一 Ax22!Ox21(83u+ 3!Ox3刊八〃—— Ax31(O1(O4u'+ 4!IOx4nAx4+…，(8)Unm-1=uUnm-1=unmOx)mnA1(O2U、nAx+ 2!Ox2Ax21(O3"丫—3!"Ox3)mAx31(O4u丫+ 4!"Ox4)將(8)與⑼相加得到二階導(dǎo)數(shù)的差分格式nun —2un+un鋁一nun —2un+un鋁一m+1 m m-1，Ax2"Ox2)m根據(jù)中心差分算子50，(5)可化為(10)」(c)50un—2un+un，鋁——x m+ m m^—2Ax3(11)所以，將(4)，(6)，(11)代入到方程(1)中，其中unmun+1—un-1m2At"'1(3un+un得到方程(1)的近似差分方程，如下50un 50u一2un+un+以u(píng)n-―xm+p—xm+ m m1=0，m2Ax 2Ax3)，并化簡(jiǎn)得+unm-1紡1一％1+以u(píng)nAt m50un50^un —2un+un——x——m+p―x m+ m m——=0，Ax Ax3(12)(13)最終差分格式，由matlab編寫將網(wǎng)比r=A代入上式，得Ax

1( )( )Un+1=Un-1一—raU+Un1( )( )Un+1=Un-1一—raU+Un+Un/Un—Un/邛( )云物+2一2U：+1+2U：—1一U：—2，(14)其穩(wěn)定性條件［15是—48、一rmaxa|u|+ <1。I 或2)2.2跳點(diǎn)格式(15)對(duì)于時(shí)間微分項(xiàng)號(hào)采用向前差分’有Un+1—UnmAtm；(16)至于對(duì)空間變量的微分，跳點(diǎn)格式是根據(jù)m,n的和的不同，有不同的差分格式，其敘述如下。當(dāng)m+n=偶數(shù)時(shí)，對(duì)于非線性項(xiàng)中空間微分坐，在nAt時(shí)刻采用中心差分dx格式，艮即(5)式。對(duì)于3階微分項(xiàng)竺，在nAt時(shí)刻先對(duì)2階微分采用中心差分，再對(duì)二階微分dx3進(jìn)行一階中心差分，即(11)式。將(5)式，(11)式，(16)式代入到方程(1)中，得Un+1—Un n80Un 80\n—2un+Un―m m-+au-—xm+p—xm+1 m m—1—=0，At m2Ax 2Ax3其中Un=1(Un+Un)，上式可化簡(jiǎn)，得m2m+1 m—1(17)Un+1—Un 80fn 80―m m-+a-——+8—x-At 2Ax力4 )n —2Un+Unm+1 m m—1—=U2Ax3(18)其中f=竺，化簡(jiǎn)后，得2, 1廠Un+1=, 1廠Un+1=Un——ra80f(19)一2云U+2一:+1+2un—1—Un—2。當(dāng)m+n=奇數(shù)時(shí)，對(duì)空間的微分項(xiàng)都是采用的隱式格式，即對(duì)空間的微分都是取在(n+1i\t時(shí)刻，表述如下。對(duì)于非線性項(xiàng)中的笑，仍然是在G+山時(shí)刻采用中心差分格式，得空、n空、n+1dx/mUn+1—Un+1

鋁一m+1 m一1；2Ax(20)對(duì)于3階微分項(xiàng)中的竺，在(n+1)At時(shí)刻先對(duì)2階微分采用中心差分，再進(jìn)dx3行一階中心差分，得n+】 60U行一階中心差分，得n+】 60Un+1—2un+1+Un+1)2Ax3其中u+1m、8x3/m將(16)式，(20)式，(21)式代入到方程(1)中得Un+1—Un n+160Un+1 60Un+1—2"n+1+Un+1―m m?+au ―xm—+p―xm+ m m——=0，At m2Ax 2Ax3n+1+Un+1)，上式可化簡(jiǎn)，得2m+1 m—1un+1— un60fn+1 60 ^Un+1 —2un+1 +Un+1—m m +以一—x m +p―x m+1 m m—1— =0，At 2Ax 2Ax3其中f=竺，化簡(jiǎn)后，得21 1rP( )Un+1=Un——尸0/60fn+1— Un+1—2Un+1+2Un+1—Un+1，。mm2xm2Ax2 m+2 m+1m—1 m-2跳點(diǎn)格式的穩(wěn)定性條件［15是rau—^<1。Ax22.3分裂格式(21)(22)(23)(24)(25)將非線性項(xiàng)和頻散項(xiàng)分開后獨(dú)立計(jì)算。對(duì)流方程告章u齊=。可用迎風(fēng)格式或Lax等耗散格式，頻散方程*+p°巴=0可取一般格式如中心差分格式，Ot Ox3CN格式等。例如對(duì)流方程取Lax格式，頻散方程取時(shí)間前差空間中心差分格式，則得Un=—^Un+Un )一-!-ra^fn—fn，，m2 m+1 m—1 2m+1m—11rpum+=um—2Ax2'「U—2un+unm+1 m m—1V )精度為O(At+Ax2)，穩(wěn)定性條件用(15)式，當(dāng)系數(shù)P較大時(shí)，頻散項(xiàng)孕起重要Ox3用半隱式格式，得的作用，如該項(xiàng)用顯式格式，要求At?Ax3,而用全隱式增加計(jì)算量，為此可采用半隱式格式，得Un+1—un 60fn Un—Un+Un—1)+^Un+1+Un—1)—Un―m m+a-—x——m+P―mI2 m+ m+ m— m— m2=0， (26)At 2Ax 2Ax3化簡(jiǎn)得1Un+1Un+1=um——ra5mm2—2Ax2m+20fnxm^.n +un—1 (n+1+un—1)—un 1，m+1 m+1 m—1 m—1 m—2(27)這是一個(gè)三層隱式格式，穩(wěn)定性條件［15為(28)raU<1。(28)3.實(shí)例對(duì)于方程(1)，系數(shù)和參數(shù)的選取:-=1.0，6=4.84x10-4，C=0.3，D=—6.0;時(shí)間步長(zhǎng)M=0.0001，空間步長(zhǎng)軟=0.02。根據(jù)系數(shù)之間的關(guān)系aCA= -=12.448\：6B=aAC=3.7345。根據(jù)行波解的特性，此時(shí)孤立波向前傳播的速度可以求出v=-=0.3000A用Matlab對(duì)三種差分格式進(jìn)行編程計(jì)算(程序見(jiàn)附錄)，至于初始值和邊界值，是根據(jù)精確解(即行波解)給出的，即TOC\o"1-5"\h\zu0=u(x,t)=3Csech2^Ax+D)m m0、 mu1=u(x,t)=3Csech2(Ax—Bt+D)m,m1、 /m1、(29)un=u(x,t)=3Csech2(Ax—Bt+D)(29)v0 /0n、 /0 幾、un=uG,t)=3Csech2(Ax—Bt+D)un=utx ,t)=3Csech2(Ax —Bt+D)M—1,M—1、n , M—1 n、un=uvx,t)=3Csech2(Ax—Bt+D)。計(jì)算后作出t=0.1時(shí)，精確解以及三種差分格式的解的圖，如下(c) (d)圖3.各種格式的解得圖：(a)精確解；(b)蛙跳格式；(c)跳點(diǎn)格式；(d)分裂格式；上面的(c)圖明顯不同于另外三幅圖，是因?yàn)槟M的結(jié)果中有負(fù)值存在，如表1中所示，而精確解(2)是一平方表達(dá)式，不可能小于零，蛙跳格式和分裂格式在其他時(shí)刻某些位置也有負(fù)值(如附表2)，只是在t=0.1時(shí)所取的點(diǎn)處沒(méi)有負(fù)值。表1給出了t=0.1時(shí)的振幅，表2給出了x=2的位置不同時(shí)刻的振幅。表1.t=0.1時(shí)不同位置的振幅x精確解蛙跳跳點(diǎn)分裂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-0052.0717e-0052.0292e-0052.0246e-0050.61.3102e-0071.4048e-0071.3921e-0071.3733e-0070.79.0125e-0109.4558e-0109.6225e-0109.4272e-010

0.86.1996e-0126.5007e-0126.6283e-0126.4889e-0120.94.2646e-0144.4718e-0144.5597e-0144.4638e-01412.9336e-0163.0761e-0163.1366e-0163.0706e-0161.12.018e-0182.116e-0182.1576e-0182.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.4708e-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.x二2的位置不同時(shí)刻的振幅t精確解蛙跳跳點(diǎn)分裂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-0161.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.8569e-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-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.3001e0140.26.1912e-0161.2862e-0123.2994e-0126.8679e-015f=0.1時(shí)，三種差分格式的誤差曲線沿空間的變化趨勢(shì)大體相同，如圖4所示。

由圖4得，在x>1的范圍內(nèi)誤差很小，幾乎緊貼零線，表明差分格式是收斂的;在xet),1］內(nèi)，誤差的變化波動(dòng)較大。在e!d,1］內(nèi)三種差分格式的誤差變化如下圖5。*0.00018| 圖5.xeh1］內(nèi)三種差分格式的誤差變化圖4.附錄：Matlab源碼蛙跳格式程序:運(yùn)行首先工作取決Debugrun工作空間出現(xiàn)所有參數(shù)右擊查看參數(shù)信息%蛙跳格式clearall，clc清空%系數(shù)參數(shù)賦值a=1.0;b=4.84e-4;C=0.3;D=-6;%取步長(zhǎng)tb=0.0001; %時(shí)間步長(zhǎng)xb=0.02; %空間步長(zhǎng)倍數(shù)可改變調(diào)試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;精確解空間網(wǎng)格點(diǎn)數(shù)fori=1:mmax %利用精確解給出初值，n=0和n=1時(shí)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;endforj=1:nmax %利用精確解給出邊值，m=0,m=1,m=199,m=200t(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(mmaxT,j)=3*C*(sech(A*x(mmaxT)-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;endform=1:mmax %計(jì)算精確解x(m)=(m-1)*xb;forn=1:nmaxt(n)=(n-1)*tb;uj(m,n)=3*C*(sech(A*x(m)-B*t(n)+D))“2;endendforn=3:nmax %蛙跳格式計(jì)算form=3:mmax-2u(m,n)=u(m,n-2)+k1*(u(m+1,nT)+u(m,nT)+u(mT,nT))*(u(m+1,nT)-u(mT,nT))+k2*(u(m+2,nT)-2*u(m+1,nT)+2*u(mT,nT)-u(m-2,nT));endend跳點(diǎn)格式程序：%跳點(diǎn)格式clearall,clc%系數(shù)參數(shù)賦值a=1.0;b=4.84e-4;C=0.3;D=-6;tb=0.0001; %時(shí)間步長(zhǎng)xb=0.02; %空間步長(zhǎng)A=0.5*sqrt(a*C/b);B=a*A*C;k1=-1/4*a*tb/xb;k2=-tb*b/(2*(xb)”3);mmax=201;nmax=2001;fori=1:mmax %利用精確解給出初值條件，t=0時(shí)x(i)=(i-1)*xb;u(i,1)=3*C*(sech(A*x(i)+D))“2;uju(i,1)=0;endforj=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(mmaxT,j)=3*C*(sech(A*x(mmaxT)-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;endform=3:mmax-2 %迭代中未知點(diǎn)的初始賦值forn=2:nmaxu(m,n)=0.5;uu(m,n)=0.6;endendforn=2:nmax %跳點(diǎn)格式計(jì)算form=3:mmax-2 %計(jì)算顯式值if(mod((m+n),2)==0)u(m,n)=u(m,nT)+k1*((u(m+1,nT))“2-(u(mT,nT))“2)+k2*(u(m+2,nT)-2*u(m+1,nT)+2*u(mT,nT)-u(m-2,nT));uu(m,n)=u(m,n);endendj=0；while(1