偏微分方程數(shù)值解

1、1偏微分方程數(shù)值解偏微分方程數(shù)值解(numerical solution of partial differential equations)主講：王曰朋主講：王曰朋2參考數(shù)目1. george j. haltiner, roger terry williams, numerical prediction and dynamic meteorology(2nd edition), the united states of america, 1979.2. curtis f.gerald and patrick o., applied numerical analysis, person edu

2、cation, inc., 2004.3. eugenia kalnay, atmospheric modeling, data assimilation and predictability, the press syndicate of the university of cambridge,2003.4. arieh iserles, a first course in the numerical analysis of differential equations, cambridge university press,1996.5. 李榮華，馮國(guó)忱. 微分方程數(shù)值解. 北京：人民教育

3、出版社，1980.6. 徐長(zhǎng)發(fā)，李紅. 實(shí)用偏微分方程數(shù)值解法. 華中科技大學(xué)出版社，2003.7. 沈桐立，田永祥等. 數(shù)值天氣預(yù)報(bào). 北京：氣象出版社，2007.3數(shù)值天氣預(yù)報(bào)pde數(shù)值解1. 挪威氣象學(xué)家v.bjerknes（1904）提出數(shù)值預(yù)報(bào)的思想：通過(guò)求解一組方程的初值問(wèn)題可以預(yù)報(bào)將來(lái)某個(gè)時(shí)刻的天氣思想；2. l.f.richardson(1922)：開(kāi)創(chuàng)了利用數(shù)值積分進(jìn)行預(yù)報(bào)天氣的先例，由于一些原因（如，計(jì)算穩(wěn)定性問(wèn)題“courant,1928”）并沒(méi)有取得預(yù)期的效果嘗試；3. charney, fjortoft, and von neumann(1950), 借助于princ

4、eton大學(xué)的的計(jì)算機(jī)(eniac)，利用一個(gè)簡(jiǎn)單的正壓渦度方程（c.g.rossby,1940）對(duì)500mb的天氣形式作了24小時(shí)預(yù)報(bào)-成功；4the electronic numerical integrator and computer (eniac).5常微分方程的數(shù)值解大氣科學(xué)中常微分方程和偏微分方程的關(guān)系 1. 大氣行星邊界層(近地面具有湍流運(yùn)動(dòng)特性的大氣薄層,11.5km), ?？寺╲.w.ekman）（瑞典）螺線(xiàn)的導(dǎo)出； 2. 1963年，美國(guó)氣象學(xué)家lorenz在研究熱對(duì)流的不穩(wěn)定問(wèn)題時(shí)，使用高截?cái)嗟淖V方法，由boussinesq流體的閉合方程組得到了一個(gè)完全確定的三階常微

5、分方程組，即著名的lorenz系統(tǒng)。6lorenz系統(tǒng)dx / dt = a (y - x) dy / dt = x (b - z) - y dz / dt = xy - c z 其中，a=10,(prandtl number); b=28(rayleigh number); c=8/3; (x,y,z)_0=(0.01;0.01;1e-10)705101520253035404550-30-20-10010203040508-20020-30-20-1001020300102030405091011franceshini 將navier-stokes方程截?cái)酁槲寰S的截譜模型如下：112345

6、221 33312441655142449357re53xxx xx xxxx xxxx xxxx xxxx x = -+= -+ = -+= -= - 12歐拉法折線(xiàn)法 常微分方程能直接進(jìn)行積分的是少數(shù)，而多數(shù)是借助于計(jì)算機(jī)來(lái)求常微分方程的近似解； 有限差分法是常微分方程中數(shù)值解法中通 常有效的方法； 建立差分算法的兩個(gè)基本的步驟： 1. 建立差分格式，包括：a. 對(duì)解的存在域剖分；b. 采用不同的算法可得到不同的逼近誤差截?cái)嗾`差（相容性）；c.數(shù)值解對(duì)真解的精度整體截?cái)嗾`差（收斂性）；d.數(shù)值解收斂于真解的速度；e. 差分算法舍人誤差（穩(wěn)定性）.132.差分格式求解將積分方程通過(guò)差分方程轉(zhuǎn)

7、化為代數(shù)方程求解，一般常用遞推算法。 在常微分方程差分法中最簡(jiǎn)單的方法是euler方法，盡管在計(jì)算中不會(huì)使用，但從中可領(lǐng)悟到建立差分格式的技術(shù)路線(xiàn)，下面將對(duì)其作詳細(xì)介紹:14差分方法的基本思想“就是以差商代替微商”23( )1111()( )( )( )( )( )()2!3!nnnu thu tu t hu t hut hut ho hn23( )1111()( )( )( )( )( )()2!3!nnnu thu tu t hu t hut hut ho hn考慮如下兩個(gè)taylor公式：（1）（2）從（1）得到：1()( )( )( )iiiu tu tu to hh15從（2）得到：

8、1()( )( )( )iiiu tu tu to hh211()()( )()2iiiu tu tu to hh從（1）-（2）得到：從（1）+（2）得到：2112()2 ( )()( )()iiiiu tu tu tu to hh16對(duì)經(jīng)典的初值問(wèn)題0( , )(0)duf t udtuu(0, )tt1212( ,)( ,)f t uf t ul uu滿(mǎn)足lipschitz條件保證了方程組的初值問(wèn)題有唯一解。17一、算法構(gòu)造：0tutnt0tit1iitth1. 在求解域上等距離分割：2. 在 有：1 ,iit t1()( )( , )( )iiu tu tf t uo hh1( ,)i

9、iiiuuf t uh1( ,)iiiiuuhf t u微分方程的精確解差分方程的精確解183. 應(yīng)用時(shí)采用如下遞推方式計(jì)算：0u1u2u3u0t1t2t4. 例題對(duì)初值問(wèn)題2(0)1yxyy (0,1)x5,n 0.2h 用euler法求解，用即，001111223334445(,)1.000,1.200;( ,)1.600,1.520;(,)2.320,1.984;(,)3.184,2.621;(,)4.221,3.465f xyyf x yyf xyyf xyyf xyy195. euler法的幾何意義0t0tit1iitth在遞推的每一步，設(shè)定( )iiu tu1()iu t1iu2(

10、)o h( )iu t( ),iiu tu過(guò)點(diǎn)( ,)iit u作 的切線(xiàn)，該切線(xiàn)的( )u t方程為：11( )()iiiiiuuu ttt即：1( ,)iiiiuuhf t u20二、誤差分析構(gòu)造算法后，這一算法在實(shí)際中是否可行呢？也就是說(shuō)是否使計(jì)算機(jī)仿真而不失真，這還需要進(jìn)一步分析。1. 局部截?cái)嗾`差-相容性為了分析分析數(shù)值方法的精確度，常常在( )iiuu t成立的假定下，估計(jì)誤差111()iiieu tu這種誤差稱(chēng)為“局部截?cái)嗾`差”，如圖。21()( )( )( )2iiihu tu thu tu2( )( , ( )( )2iiihu thf t u tu22111()( )()2

11、iiiheu tuuo h局部截?cái)嗾`差是以點(diǎn)it的精確解( )iu t為出發(fā)值，用數(shù)值方法推進(jìn)到下一個(gè)點(diǎn)1it而產(chǎn)生的誤差。21整體截?cái)嗾`差是以點(diǎn)0t的初始值0u為出發(fā)值，用數(shù)值方法推進(jìn)i+1步到點(diǎn)1it，所得的近似值 與精確值 的偏差：2.整體截?cái)嗾`差收斂性1iu1()iu t111()iiiu tu稱(chēng)為整體截?cái)嗾`差。1( ,)iiiiuuhf t u11()( )( , ( )iiiiiu tu thf t u te21 ( , ( )( ,)iiiiiih f t u tf t udh21iiihldh12210(1)1 (1)(1)(1) iiihldhhlhlhl2110(1)(1)

12、1iidhhlhlhl220(1)tltldheel特例，00若不計(jì)初始誤差，即則1(1)tlidhel1( )io h即3.舍入誤差穩(wěn)定性假設(shè)一個(gè)計(jì)算機(jī)僅表示4個(gè)數(shù)字（小數(shù)點(diǎn)后面）， 那么10.3333,310.11119計(jì)算20.373410.1112011119110.33340.33333xxx20.33341190.6667133xxxx23我們的要求是：最初產(chǎn)生的小誤差在以后的計(jì)算中雖然會(huì)傳遞下去，但不會(huì)無(wú)限制的擴(kuò)大，這就是穩(wěn)定性所描述的問(wèn)題。下面引進(jìn)穩(wěn)定性的概念：設(shè)由初值0u0v得到精確解 ，nu由初值得到精確解 ，nv若存在常數(shù)c和充分小的步長(zhǎng)0h使得00nnuvc uv則稱(chēng)數(shù)

13、值方法是穩(wěn)定的。tu00u0vntnunv242(0)1dyxydxyy計(jì)算例題(0,1)x其解析解為：12yxx = 0 0.2000 0.4000 0.6000 0.8000 1.0000y = 1.0000 1.2000 1.3733 1.5315 1.6811 1.8269250xy26三、改進(jìn)的euler法將微分方程( )(, ( )u tf t ut在區(qū)間1 ,iit t上積分，得到11()( )( , ( )iitiitu tu tf t u t dt用梯形法計(jì)算積分的

14、近似值，有1111( , ( ) ( , ( )(, ()2iitiiiitf t u t dtf t u tf tu t于是1111()( ) ( , ( )(, ()2iiiiiiu tu th f t u tf tu t1111 ( ,)(,)2iiiiiiuuh f t uf tu這是一個(gè)隱式格式，一般需要用迭代法來(lái)求 ，而用顯式的euler法提供初值。1iu2701( ,)iiiiuuhf t u1 111/2 ( ,)(,)kkiiiiiiuuhf t uf tu為了簡(jiǎn)化計(jì)算的過(guò)程，在此基礎(chǔ)上進(jìn)一步變?yōu)槿缦滤惴ǎ?11/2 ( ,)(,)iiiiiiuuhf t uf tu1( ,

15、)iiiiuuhf t u此式稱(chēng)為“改進(jìn)的euler法。接下來(lái)討論其幾何意義預(yù)估校正其局部截?cái)嗾`差為3()o h這個(gè)問(wèn)題將在下節(jié)討論。28tu0it1it0( ,)iimf t u( )u tiu1iu111(,)iimf tu122averagemmm1iu1()iu t290xyeuler法、改進(jìn)的euler法和解析解的比較2(0)1dyxydxyy(0,1)x12yx30四、（龍格-庫(kù)塔）runge-kutta方法簡(jiǎn)單的euler法是建立在taylor級(jí)數(shù)的一項(xiàng)展開(kāi)；改進(jìn)的eu

16、ler法是以?xún)身?xiàng)taylor級(jí)數(shù)為基礎(chǔ)建立的，如：23()( )( )( )()2iiiihu thu thu tu to h23()( )( )( ()( )/ /2()iiiiiu thu thu thu thu tho h()( )( )2iiiu thu tu th11(,)( ,)2iiiiif tuf t uuh如果我們截取taylor級(jí)數(shù)的更多項(xiàng)會(huì)得到什么樣的求解方法呢？?jī)蓚€(gè)德國(guó)數(shù)學(xué)家（c.runge & m.kutta）以這種思想為基礎(chǔ)建立了求解微分方程的龍格-庫(kù)塔方法。它是常微分方程數(shù)值解法中使用最為廣泛的方法之一。31一般地，一個(gè)k階的runge-kutta方法可用

17、下面的公式表示：11111( ,)(,)kiijjjiijjijijllluuw kkhf t ukhf tc h ua k2,3,jk其中，jw是待定的加權(quán)系數(shù)，2321,1,kk kc cc aa是待定的系數(shù)。euler法就是11,1kw的r-k法。其系數(shù)的確定如下：將1iu展開(kāi)成h的冪級(jí)數(shù)，并與微分方程的精確解1()iu t在點(diǎn)it的taylor展開(kāi)式相比較，使兩者的前1p項(xiàng)相同，這樣確定的r-k法，其局部截?cái)嗾`差為1()po h，根據(jù)所得關(guān)于待定系數(shù)的方程組，求出它們的值后代入公式，就成為一個(gè)p階r-k方法。32例題以二階r-k法為例說(shuō)明上述過(guò)程11122112211( ,)(,)ii

18、iiiiuuw kw kkhf t ukhf tc h ua k21( )( )2iiiihuuhu tu t2( ,)( ,)2iiiiihuhf t uf t u211( ,)( ,)( ,) ( ,)22iiitiiuiiiiuhf t uhf t uf t uf t u把12,k k代入1iu中，有112221( ,),( ,)iiiiiiiiuuwhf t uw hf tc h ua hf t u12221( ,) ( ,)( ,)( ,)iiiiitiiuiiuwhf t uw h f t uc hf t ua hff t u21222221() ( ,)( ,)( ,)iiiti

19、iuiiuh wwf t uh w c f t uw a ff t u33經(jīng)比較得到122222111212www cw a1222212112www cac取 為自由參數(shù)：2c21 2,12 3c 121 31 1(,)(0,1),( , ),( , )4 42 2w w從而得到不同的但都是二階的r-k方法，對(duì)應(yīng)的有中點(diǎn)法、heun（亨）法以及改進(jìn)的euler法。34基于相同的過(guò)程，通過(guò)比較五次taylor多項(xiàng)式，得到更加復(fù)雜的結(jié)果，給出了包含13個(gè)未知數(shù)的11個(gè)方程。得到多組系數(shù)，其中常用的是以下四階r-k法：1123412132431()6( ,)11(,)2211(,)22(,)iii

20、iiiiiiiuukkkkkhf t ukhf th ukkhf th ukkhf th uk改進(jìn)的euler法、r-k法以及解析解的比較：3501.71.8xy2(0)1dyxydxyy(0,1)x12yx36五、線(xiàn)性多步（linear multistep method）法1. 預(yù)備知識(shí)：插值多項(xiàng)式插值是離散函數(shù)逼近的重要方法，利用它可通過(guò)函數(shù)在有限個(gè)點(diǎn)處的取值狀況，估算出函數(shù)在其他點(diǎn)處的近似值。從幾何上理解：對(duì)一維而言，已知平面上n1個(gè)不同點(diǎn)，要尋找一條n次多項(xiàng)式曲線(xiàn)通過(guò)這些點(diǎn)。插值多項(xiàng)式一般常見(jiàn)

21、的是拉格朗日插值多項(xiàng)式。2. 氣象應(yīng)用不均勻站點(diǎn)上的氣象要素?cái)?shù)據(jù)均勻網(wǎng)格點(diǎn)上的數(shù)據(jù)插值 3. 拉格朗日插值多項(xiàng)式拉格朗日插值多項(xiàng)式逼近可能是求插值節(jié)點(diǎn)不均勻的插值多項(xiàng)式的最簡(jiǎn)單的方法。實(shí)驗(yàn)觀察結(jié)果或原始測(cè)量數(shù)據(jù)的分布通常是非均勻的。例如，四個(gè)點(diǎn)可以確定一個(gè)三次多項(xiàng)式，其拉格朗日形式為：23413431212131421232412312434313234414243()()()()()()( )()()()()()()()()()()()()()()()()()()xxxxxxxxxxxxp xffxxxxxxxxxxxxxxxxxxxxxxxxffxxxxxxxxxxxx374. adams

22、-bashforth（阿達(dá)姆斯貝雪福斯）公式首先，用以下四個(gè)點(diǎn)對(duì)( , ( )f t u t進(jìn)行三次langrage插值：111222333( ,( , ( ),(,(, (),(,(, (),(,(, ()nnnnnnnnnnnntf t u ttf tu ttf tu ttf tu t則1233123123333132()()()( , )( )( , ( )()()()()()()(, ()()()()nnnnnnnnnnnnnnnnnnnnnnttttttf x ttf t u tttttttttttttf tu ttttttt于是，有1113()( )( , ( )( )( )nnn

23、ntnnttntu tu tf t u t dtu tt dt容易算出，13123( )55 ()59 ()37 ()9 ()24nnxnnnnxhx dxf xf xf xf x例如，我們可以算得11233123()()()1()(2 )(3 )()()()6nnnnxxhnnnnnnxxnnnnnnxxxxxxdxxxh xxh xxh dxxxxxxxh*2*138433011 55()(2 )(3 )5566424hhth th th dthhh將（*2）代入（*1）得到adams-bashforth公式：111223355 (,)59 (,)37 (,)9 (,)24nnnnnnnn

24、nnhuuf x uf xuf xuf xu基于同樣的計(jì)算過(guò)程，可以得到另外一個(gè)計(jì)算公式：11111229 (,) 19 (,)5 (,)(,)24nnnnnnnnnnhuuf xuf x uf xuf xu這稱(chēng)為adams-moulto公式。111223355 (,)59 (,)37 (,)9 (,)24nnnnnnnnnnhuuf x uf xuf xuf xu11111229 (,) 19 (,)5 (,)(,)24nnnnnnnnnnhuuf xuf x uf xuf xu預(yù)估校正39偏微分方程數(shù)值解主講：王曰朋40一、區(qū)域的離散1則函數(shù)可表示為：( , )( ,)u x

25、 yu ih jk1,2,3i 1,2,3j 二、1.（一維）一、二階導(dǎo)數(shù)的有限差分近似表達(dá)式422.（二維）一、二階偏導(dǎo)數(shù)的有限差分近似(,)(,)( ,)( )ijijijx yu xh yu x yuo hxh1(,)ijjjiix yuuuxh2(,)(,)(,)()2ijijijx yu xh yu xh yuo hxh11(,)2ijjjiix yuuuxh211111111(,)1()22ijjjjjiiiix yuuuuux yhk 22()()o ho k433. 拋物型方程拋物型方程初條：精確解為（以熱傳導(dǎo)或磁擴(kuò)散方程為例）（以熱傳導(dǎo)或磁擴(kuò)散方程為例）22xutu)()0,

26、(xfxu初值問(wèn)題)0,(txd)(4)(exp41),(2ftxttxu不論初始分布如何集中，它總在瞬間影響于無(wú)窮遠(yuǎn)，雖該影響隨距離按指數(shù)衰減，然而它是以無(wú)限速度傳播。此乃拋物型方程解的特征。44三、熱傳導(dǎo)方程（拋物方程）1. 熱傳導(dǎo)方程的介紹2.222(0, )( , )0( ,0)( )uuatxutu l tu xf x離散化121122jjjjjiiiiiuuuuuakh0(0,)0juujk( , )0jnuu l t0( ,0)( )iiuu ihf ihf（1）向前差分格式：45計(jì)算：111(1 2 )jjjjiiiiusus usu這是一個(gè)顯式格式（四點(diǎn)格式）j1j i1i1

27、i每一層各個(gè)節(jié)點(diǎn)上的值是通過(guò)一個(gè)方程組求解得到的。這可以從下面的計(jì)算過(guò)程看出來(lái)。22kash461000121010002321100034321000112(1 2 )(1 2 )(1 2 )(1 2 )nnnnusus usuusus usuusus usuusus usu1 21 21 21 2ssssssss系數(shù)矩陣為4722( ,0)sin()(0, )(1, )01;00.1uutxu xxututxt 計(jì)算實(shí)例：4800.020.040.060.080.100.51-3-2-10123txu492. 向后差分格式111121122jjjjjiiiiiuuuuuakh11111(1

28、2 )jjjjiiiisus usuu當(dāng)知道第n層上的 時(shí)，要確定第n+1層上各點(diǎn)值 必須通過(guò)求解一個(gè)線(xiàn)性代數(shù)方程組。jiu1jiu11121011113212111432311111(12 )(12 )(12 )(12 )jjjjjjjjjjjjjjjjnnnnsus usuusus usuusus usuusus usuu50111122111112121212jjjjjjnnjjnnuussuusssssuussuu其矩陣表達(dá)式如下：j1j 1ii1i51這是一個(gè)古典四點(diǎn)向后差分格式。計(jì)算實(shí)例22( ,0)sin()(0, )(1, )01;00.1uutxu xxututxt 5200

29、.020.040.060.080.100.50.81txu533. crank-nicolson格式，亦稱(chēng)六點(diǎn)對(duì)稱(chēng)格式1211111112222()2jjjjjjjjiiiiiiiiuuauuuuuukhh1111111(12 )(1 2 )jjjjjjiiiiiisus ususus usu11121111/2/21/21/212jjjnjnussusssssussu54j1j 1i1i1211/2/21/21/2/21 2jjjnjnussusssssussu5500.020.040.060.080.100.50.81txu564. richar

30、dson格式11211222jjjjjiiiiiuuuuuakh11112 (2)jjjjjiiiiius uuuu1j 1i1ij1j i這是一個(gè)五點(diǎn)三層差分顯式格式57討論：假若由于 的作用，導(dǎo)致差分方程的近似解設(shè)為：jiu 于是，我們可得到差分格式的誤差方程如下：11112 (2)jjjjjiiiiisxtrichardson格式是不穩(wěn)定的。585. 穩(wěn)定性判別 von-neumann 穩(wěn)定性在判斷有限差分近似的穩(wěn)定性方法中，以von-neumann方法使用較為廣泛，它僅適用于線(xiàn)性常系數(shù)的有限差分近似。其過(guò)程如下：首先，要研究的差分方程可寫(xiě)為：101nnmj mmj mm nm na u

31、b u如，其次，對(duì)進(jìn)行變量分離：jiu2expnnjpjppuvixl111(1 2 )nnnnjjjjusus usu59expnnjjuvix最后將代入所考察的有限差分方程。111expexp(1 2 )expexpnnjjnnjjvi xsvi xs vi xsvi x1exp(1 2 )expnnnnvsvi hs vsvi h1(cossin)(1 2 )(cossin)1 2 (1 cos)nnvshihsvshihsh 定義為放大系數(shù)6021 4 sin2hs 下面說(shuō)明，在什么條件下能使11nnvv對(duì)所有的成立。21 1 4 sin12hs 從上式，我們看出，21 1 4 sin

32、2hs 212 sin2sh6112s . 交替顯隱式格式（）顯式預(yù)測(cè)隱式校正格式在n+1/2層上用古典顯式格式計(jì)算出過(guò)度值，再在n+1層上用12nju古典隱格式校正預(yù)測(cè)值，即：1/ 211211/ 21111122222nnnnnjjjjjnnnnnjjjjjuuuuuhuuuuuh62(2). 跳點(diǎn)格式首先將網(wǎng)格點(diǎn)（j, n）按j+n等于偶數(shù)或奇數(shù)分成兩組，分別稱(chēng)為偶數(shù)網(wǎng)點(diǎn)和奇數(shù)網(wǎng)點(diǎn)。從 到的計(jì)算過(guò)程中，先在偶數(shù)網(wǎng)格點(diǎn)上用古典顯式格式計(jì)算，再在奇數(shù)網(wǎng)點(diǎn)上用古典隱格式計(jì)算，即：nt1nt111211111122121nnnnnjjjjjnnnnnjjjjjuuuuuhnjuuuuuhnj偶

33、數(shù)奇 數(shù)63222220uuatx0uuatx（a）一階常系數(shù)線(xiàn)性雙曲型方程（b）二階常系數(shù)線(xiàn)性雙曲型方程（波動(dòng)方程） 其中a為常數(shù)為為 這些方程的定解條件，可僅有初始條件，也可以有初始條件和邊界條件。其中a為常數(shù) 同橢圓型方程與拋物型方程相比，雙曲型方程差分格式的性質(zhì)與定解問(wèn)題解析解的性質(zhì)有更密切的關(guān)系。64考慮0dtdxatudtdxxutudtduxxuatu0 由于u(x,t)沿x-t平面上方向?yàn)閐x/dt=a的直線(xiàn) xat=c（c為常數(shù)）的變化率為0，即故沿x-t平面上任一條斜率為1/a的直線(xiàn)xat=c，u(x,t)為常數(shù)。平行直線(xiàn)族xat=c就是方程(3.1)的特征線(xiàn)。(3.2)(

34、3.1)xxxu)()0 ,(65 利用特征線(xiàn)，可以求出初值問(wèn)題(3.1)、(3.2)的解：)(),(atxtxu),(00tx 由于u(x,t)在點(diǎn) 處的值依賴(lài)與(x)在點(diǎn) 的值，故初始線(xiàn)t=0上的點(diǎn) 稱(chēng)為解u(x,t)在點(diǎn) 的依賴(lài)區(qū)域。00atx 00atx ),(00tx 與拋物型方程求解類(lèi)似，對(duì)x-t平面進(jìn)行矩形網(wǎng)格剖分，x方向的步長(zhǎng)為h，t方向的步長(zhǎng)為，網(wǎng)點(diǎn) 簡(jiǎn)記為(j,k)。),(kjtx66對(duì)方程(3.1)，利用差商代替導(dǎo)數(shù)的方法，可得011huuauukjkjkjkj011huuauukjkjkjkj02111huuauukjkjkjkj, 2, 1, 0, 2, 1, 0kj

35、 前兩個(gè)格式的局部截?cái)嗾`差階為 ，分別稱(chēng)為左、右偏心格式。)(ho)(2ho 第三個(gè)格式的截?cái)嗾`差階為 ，它稱(chēng)為中心差分格式。kjkjkjurarauu)1 (11kjkjkjrauurau11)1 ()(2111kjkjkjkjuuaruu其中hr/即67 從差分格式依賴(lài)區(qū)域和微分方程依賴(lài)區(qū)域的關(guān)系，可以得到差分格式收斂的必要條件： 對(duì)于左偏心格式， cfl條件為：a0，且 。1ar1ar對(duì)于右偏心格式， cfl條件為：a0（或a0， 上網(wǎng)點(diǎn)a(j1,k)，b(j,k)，c(j+1,k)處的解值已經(jīng)算出，從點(diǎn)p(j,k+1)作特征線(xiàn)，它與線(xiàn)段ab交于點(diǎn) d。ktt p1kttktt cbad

36、由u(p)=u(d)，有)(2)(2)(auhahcuhahdu, 2, 1, 0, 2, 1, 0kj)1 ()1(21111kjkjkjuaruaru這樣，得到lax格式：當(dāng) ，lax格式穩(wěn)定，截?cái)嗾`差階為 。1ar)(22hho69 對(duì)方程(3.1)，利用特征線(xiàn)作二次插值，即可得到laxwendroff格式： 當(dāng) ， laxwendroff格式穩(wěn)定，它的截?cái)嗾`差階為 。1ar)(22hokjxkjkjkjkjurauuaruu2221112)(270應(yīng)該注意：邊值條件的給法與其它兩類(lèi)方程不同。 如果 a0，方程特征線(xiàn)向右傾，只能在 x變化區(qū)域的左邊界上給出邊界條件：)(), 0(1ttu

37、 如果 a0，方程特征線(xiàn)向左傾，只能在 x變化區(qū)域的右邊界上給出邊界條件，即使 x 的變化區(qū)域?yàn)? x d ，也只能給出邊界條件：)(),(2ttdua0xoy0 x0 ，考慮下面模型問(wèn)題：前面建立的幾個(gè)顯格式，都適用于這個(gè)問(wèn)題。ttttuxxxuxttxuatu0)(), 0(0)()0 ,(000下面建立隱格式。72連同初始條件與邊界條件：, 2, 1, 2, 1, 001111jkhuuauukjkjkjkjkjkjkjuaruararu111111該格式的局部截?cái)嗾`差階為 。)(ho令 ，格式可改寫(xiě)為hr/該格式可在0 x, t 內(nèi)所有網(wǎng)點(diǎn)上顯示地計(jì)算解之近似值。, 2, 1, 0,

38、2, 1,00kjuukkjj),(kj) 1,(kj) 1, 1(kj73 然后用中心差商逼近這些導(dǎo)數(shù)值，則可得到wendroff格式：在點(diǎn) 處，用)21,21(kjpp1kttktt chadgfbejhx hjx) 1( geptututu21fhpxuxuxu21, 2, 1, 2, 1, 0022111111111jkhuuhuuauuuukjkjkjkjkjkjkjkj）（）（74連同初始條件與邊界條件：11111()1kkkkjjjjaruuuuar 該格式的局部截?cái)嗾`差階為 ，且無(wú)條件穩(wěn)定。22()oh令 ，格式可改寫(xiě)為/rh該格式可在0 x, t 內(nèi)所有網(wǎng)點(diǎn)上顯示地計(jì)算解之近

39、似值。00,1, 2,0, 1, 2,kkjjuujk75 222220uuaxatx 為正常數(shù)考察 對(duì)x-t平面進(jìn)行矩形網(wǎng)格剖分，x方向的步長(zhǎng)為h，t方向的步長(zhǎng)為，網(wǎng)點(diǎn) 簡(jiǎn)記為(j,k)。),(kjtx 用二階中心差商代替(3.3)中的二階導(dǎo)數(shù)，則得到網(wǎng)點(diǎn)(j,k)處的差分方程： (3.3)(3.4)( ,0)( )( ,0)( )uu xxxtxt 222221kktjxjauuh12222111()2(1)kkkkkjjjjjua ruua ruu, 2, 1, 0, 2, 1, 0;/kjhr其中 ?；?1kk1k1jj1j76該格式穩(wěn)定的充分條件為 。1ar初始條件)()0 ,(xx

40、u離散：)()0 ,(txtu初始條件 離散：由0222022111)(21jxjtjjjuhauuu消去 ，得1jujjjjjararu)1 (22211221)(22ho 上述差分格式與初始條件的截?cái)嗾`差階均為 。jju0取為77ttxxuatu0, 10022222上述方法也可用于求解初邊值問(wèn)題： 10)()0 ,()()0 ,(xtxtuxxuttttuttu0)(), 1 ()(), 0(78 798081在條件下2122ssldsdsdydsdx求使得泛函達(dá)到最大的函數(shù) 。)(),(sysx2121),(ssdsdsdxydsdyxyxs 在長(zhǎng)度一定的所有平面封閉曲線(xiàn)中，求所圍面積

41、為最大的曲線(xiàn)。82當(dāng)求泛函在一個(gè)函數(shù)集合k中的極?。ɑ驑O大）問(wèn)題，則該問(wèn)題稱(chēng)為變分問(wèn)題?？疾靔(x)的極值情況。為實(shí)常數(shù)）fllfxlxxj, 0(2)(2設(shè)求 ，使rx 0)(min)(0 xjxjrx與求解方程 lx = f 等價(jià)。83對(duì)稱(chēng)正定nnnnnnaaaaaaaaaa212222111211求j(x)取極小值的駐點(diǎn)， 其中 ninjniiijiijnxbxxaxxxjxj1112121),()(設(shè)設(shè) 則j(x)可表示為：tnxxxx),(21tnbbbb),(21),(),(21)(xbxxaxj84設(shè)矩陣a對(duì)稱(chēng)正定，則下列兩個(gè)命題等價(jià)：求 ，使nrx 0)(min)(0 xjxj

42、nrx(a)(b) 是方程 的解0 xbxa：),(),(21)(xbxxaxj其中 。85 考察一根長(zhǎng)為 l 的弦，兩端固定在點(diǎn) a(0,0)和b(l,0)。當(dāng)沒(méi)有外力作用時(shí)，它的位置沿水平方向與x軸重合。設(shè)有強(qiáng)度為f(x)的外荷載垂直向下作用在弦上，于是弦發(fā)生形變。假定荷載很小，因而發(fā)生的形變也很小。用 u(x) 表示在荷載f(x)的作用下弦的平衡位置。u)0 , 0(a)0 ,(lbx86求弦的平衡位置歸結(jié)為求解兩點(diǎn)邊值問(wèn)題：設(shè)弦處于某一位置u=u(x)，可得到其總位能為 ldxufutuj02)2(21)(： 0)(0)0(0)()(luulxxfxut其中t是弦的張力。)(xuu 弦

43、的平衡位置 (記為 )將在滿(mǎn)足邊值條件 u(0)=0，u(l)=0 的一切可能位置中，使位能取極小值。弦的平衡位置 是下列變分問(wèn)題的解)(xuu)(min)(20ujujcu87 在數(shù)學(xué)上，要將某個(gè)微分方程的定解問(wèn)題轉(zhuǎn)化為一個(gè)變分問(wèn)題求解，必須針對(duì)已給的定解問(wèn)題構(gòu)造一個(gè)相應(yīng)的泛函，并證明定解問(wèn)題的解與泛函極值問(wèn)題的解等價(jià)。88考察二階常微分方程邊值問(wèn)題：0)(0)(),(buaubaxfqudxdupdxdlu),(),(21)(ufuuuj),(),(21)(ufuluujbababafudxdxquudxdxdupdxd221badxfuquup)2(2122dxquvdxdvdxdupv

44、uba),(引入泛函算子則89 與前述二階常微分方程邊值問(wèn)題相應(yīng)的變分問(wèn)題是),(),(21)(ufuuuj其中0)(,)()(,11aubahxuxubahe求 ，使1ehu )(min)(1ujujehubadxfuquup)2(2122900)(, 0)(),(buaubaxfqudxdupdxdlu設(shè) ， 是邊值問(wèn)題)(0icf 2*cu 的解，則 使 j(u) 達(dá)到極小值；*u12ehcu 反之，若 使 j(u) 達(dá)到極小值，則 是邊值問(wèn)題的解。*u),(),(21)(ufuuuj其中 是強(qiáng)制邊界條件， 是自然邊界條件，區(qū)別這兩類(lèi)邊界條件在用有限元方法求解邊值問(wèn)題時(shí)很重要。0)( b

45、u0)(au91對(duì)兩點(diǎn)邊值問(wèn)題：0)(, 0)(),(buaubaxfqudxdupdxdlu其中1*ehu ，且滿(mǎn)足變分方程： 設(shè) ，以v乘方程兩端，沿a,b積分，并利用 ，得變分方程1ehv0)(, 0)(bvavdxquvdxdvdxdupvuba),(0),(),(vfvu對(duì)任意1ehv0),(),(*vfvu在力學(xué)里， 表示虛功),(),(vfvu1ehv2*cu 設(shè) ，則 是邊值問(wèn)題解的充要條件是：*u92 。若將微分方程化為相應(yīng)的變分問(wèn)題或變分方程，則只需處理強(qiáng)加邊界條件，無(wú)需處理自然邊界條件（自然邊界條件已包含于變分問(wèn)題中泛函的構(gòu)造或已包含于給出的變分方程之中）。這一特點(diǎn)對(duì)研究

46、微分方程離散化方法及其數(shù)值解帶來(lái)了極大的方便。93模型方程其中g(shù)是平面有界區(qū)域。0),(),(|2222ugyxyxfyuxuu),(),(21)(ufuuuj),(),(21)(ufuuujgdxdyyvyuxvxuvu)(),(引入泛函算子則94與前述二階橢圓邊值問(wèn)題相應(yīng)的變分問(wèn)題是求 ，使)(10ghu )(min)()(10ujujghu其中0),(),(),(| ),()(110yxughyxuyxugh),(),(21)(ufuuuj),(),(21ufuu95 對(duì)，泛函是一樣的，只是邊界條件要作為強(qiáng)加邊值條件加在所取的函數(shù)類(lèi)上。 設(shè) ， 是二階橢圓邊值問(wèn)題的解，則 使 j(u)

47、達(dá)到極小值；)(0icf )(2*gcu *u)()(102ghgcu 反之，若 使 j(u) 達(dá)到極小值，則 是二階橢圓邊值問(wèn)題的解。*u),(),(21)(ufuuuj其中 對(duì)，二次泛函形式相對(duì)于第一邊值問(wèn)題有所改變，但函數(shù)類(lèi)的選取與邊界條件無(wú)關(guān)。96問(wèn)題000),(),(21|aaunuugyxyxfu其中 設(shè) ，以v乘方程兩端后在g上積分，并利用green公式，得變分方程)(1ghve2),(auvdsdxdyyvyuxvxuvug0),(),(vfvu0),(),(),(| ),()(111yxughyxuyxughe)(1ghve97在力學(xué)里， 表示虛功),(),(vfvu)(2g

48、cu 設(shè) 是邊值問(wèn)題的解，則對(duì)任意 ， 滿(mǎn)足變分方程。)(1ghveu 反之，若 ，且對(duì)任意 滿(mǎn)足變分方程，則 為邊值問(wèn)題的解。)()(12ghgcue)(1ghveu 與極小位能原理類(lèi)似，第一類(lèi)邊界條件為強(qiáng)加邊界條件，第二、三類(lèi)邊界條件為自然邊界條件。 虛功原理比極小位能原理應(yīng)用更廣。98：求解相應(yīng)的變分問(wèn)題或相應(yīng)的變分方程。 這兩種算法統(tǒng)稱(chēng)為ritz-galerkin方法。 以下用v表示 等sobolev空間，l表示微分算子，(u,v)為由l及邊值條件決定的雙線(xiàn)性泛函。110,hhhe 用有限維空間的函數(shù)代替變分問(wèn)題(或變分方程)中無(wú)限維空間的函數(shù)，從而在有限維函數(shù)空間中求變分問(wèn)題(或變分

49、方程)的近似解，并要求當(dāng)有限維空間的維數(shù)不斷增加時(shí)，有限維近似解逼近原變分問(wèn)題(或變分方程)的解。99由極小位能原理得出的變分問(wèn)題為：niiincu1求 ，使vu )(min)(ujujvu其中 ，),(),(21)(ufuuuj 設(shè) 是v 的n維子空間， 是 的一組基底(稱(chēng)為基函數(shù)) 。 中任一元素 可表示為nvn,21nvnvnu即選擇適當(dāng)?shù)?，使 取極小值。nccc,21)(nuj求 ，使nnvu )(min)(vjujnvvn：100展開(kāi))(nuj令nccc,21則 滿(mǎn)足解出 代入 ，則得), 1(nicinuniiincu1njfcjniiji, 2, 1),(),(1ninjnjj

50、jjijicfcc111),(),(21njcujjn, 2, 10)(),(),(21)(nnnnufuuuj101 根據(jù)最小位能原理構(gòu)造相應(yīng)于微分方程或物 理問(wèn)題的變分問(wèn)題； 取 作為 的一組基底，即用 近似代替無(wú)窮維空間v；), 1(nii,21nnspanvnv 根據(jù)二次函數(shù)取極值的必要條件，得到 nniiinvcu1ic中 所滿(mǎn)足的方程組： 求解關(guān)于 的線(xiàn)性代數(shù)方程組。icnjfcjniiji, 2, 1),(),(1102由虛功原理得出的變分方程為：niiincu1 設(shè) 是v 的n維子空間， 是 的一組基底(稱(chēng)為基函數(shù)) 。 中任一元素 可表示為nvn,21nvnvnu：0),()

51、,(vfvu0),(),(vfvun求 ，使對(duì) ， 滿(mǎn)足nnvu nvvnu103nvv由 的任意性，取 作為v ，則得n,21將 代入變分方程，則nniiinvcu1),(),(1vfvcniiinvv解出 代入 ，則得), 1(nicinuniiincu1njfcjniiji, 2, 1),(),(1104 根據(jù)虛功原理構(gòu)造相應(yīng)于微分方程或物理問(wèn) 題的變分方程； 取 作為 的一組基底，即用 近似代替無(wú)窮維空間v；), 1(nii,21nnspanvnv 求解關(guān)于 的線(xiàn)性代數(shù)方程組。icnjfcjniiji, 2, 1),(),(1nniiinvcu1icn,21 取 作為v ，將 代 入變

52、分方程，得到 滿(mǎn)足的方程組：105 基函數(shù)選取必須滿(mǎn)足強(qiáng)加邊界條件，因此 選取困難； 計(jì)算量、存儲(chǔ)量巨大； 方程組求解病態(tài)嚴(yán)重。 充分發(fā)揮了變分形式和ritz-galerkin方法的 優(yōu)點(diǎn)； 擺脫了傳統(tǒng)的基函數(shù)取法； 各種問(wèn)題的結(jié)構(gòu)程序格式統(tǒng)一。106 有限元方法基于變分原理， 又具有差分方法的一些特點(diǎn)，并且適于較復(fù)雜的區(qū)域和不同粗細(xì)的網(wǎng)格。 差分法解偏微分方程，解得的結(jié)果就是準(zhǔn)確解u在節(jié)點(diǎn)上的近似值； ritz-galerkin方法得到近似的解析解，但對(duì)一般區(qū)域，卻往往難以實(shí)現(xiàn)。 有限元方法與傳統(tǒng)ritz-galerkin方法的差別在于有限維函數(shù)空間的構(gòu)造方法。107考慮兩點(diǎn)邊值問(wèn)題：0)(

53、, 0)()(buaubxafqudxdupdxdlu 將區(qū)間a,b分割為n個(gè)子區(qū)間 。),.2 , 1(,1nixxii第i個(gè)單元記為 ，其長(zhǎng)度 。,1iiixxi1iiixxh)(),()()(1的多項(xiàng)式上為次數(shù)不超過(guò)在且mixuihxuxuvihehhh設(shè)則 稱(chēng)為試探函數(shù)空間， 稱(chēng)為試探函數(shù)。hvhhvu 108設(shè)在節(jié)點(diǎn)上試探函數(shù) 在節(jié)點(diǎn)上的一組值為hu最簡(jiǎn)單的試探函數(shù)空間 由分段線(xiàn)性函數(shù)組成。hvnuuuu,., 0210在第i個(gè)單元 上的線(xiàn)性插值函數(shù)為,1iiixxiiiiiiiihixuhxxuhxxxu11)(iiiiiiiiihixuxxxxuxxxxxu1111)(即 當(dāng) 時(shí)

54、， 的（線(xiàn)性）插值公式稱(chēng)為（線(xiàn)性）單元形狀函數(shù)。iix)(xuh109 把每個(gè)單元形狀函數(shù)合并起來(lái)，就得到整個(gè)區(qū)間a,b 上都有定義的函數(shù) ：)(xuhnnnnnnniiiiiiiiiiiiiihixuhxxuhxxixuhxxuhxxixuhxxuhxxixhxxuhxxxu111111111110011)(nx1ixix1ix0 x1x1uiu1iu1iunu110為使分段插值標(biāo)準(zhǔn)化，通常用仿射變換顯然iihxx 1 , 0把 變到 ，令,1iixxx1) 1 (0)0(0) 1 (1)0(1100nnnniiihixununxu)()()(110)(1)(10nn則iiiiiiihixu

55、hxxuhxxxu11)(變?yōu)榛?1 , 0)()()(110iihununu111定義基函數(shù)系)(xi,01, 2, 1,)(1111111iiiiiiiiiiixxxnixxxhxxxxxhxxx,0,)(1010110 xxxxxxhxxx,0,)(111nnnnnnnxxxxxxhxxx112)(),.,(),(10 xxxn 線(xiàn)性無(wú)關(guān)，它們可組成試探函數(shù)空間的基，常稱(chēng)為節(jié)點(diǎn)基函數(shù)。幾何形狀如圖)(0 x)(xia)(xn0 x1, 1nib1x1ixix1ix1nxnx任一試探函數(shù) 可表示為hhvu niiihbaxxuxu0,)()( 用這類(lèi)插值型基函數(shù)，可以構(gòu)造出適合各種邊界條件

56、的試探函數(shù)。113若借助前述放射變換iihxx1節(jié)點(diǎn)基函數(shù)可用變量表示為1.,210,)(,)()(110111nihxxxxxnhxxxxxnxiiiiiiiii，其它其它0,)()(101000hxxxxxnx別處0,)()(111nnnnnhxxxxxnx114(a) 把表達(dá)式 代入泛函；niiihxuu1)(b) 將泛函表達(dá)式中積分區(qū)間a,b變到0,1；(c) 由 達(dá)到極小值的條件)(hujnjuujjh,.,2 , 10)(得到含 的), 2, 1(,11njuuujjj), 2, 1(01, 11, 1njbuauauajjjjjjjjjj這兒00, 10 nnau),.,1(ni

57、ui(d) 解出有限元方程的數(shù)值解 ，就得到使二次泛函取極小的近似函數(shù)(有限元解)()(1xuxuiniih115有限元方程可用矩陣表示為buk其中稱(chēng)為總剛矩陣。nnnnaaaaaaak, 13222122111116 工程中形成有限元方程時(shí)，通常先在每個(gè)單元上形成單元矩陣（稱(chēng)為單元?jiǎng)偠染仃嚕?，然后由單元?jiǎng)偠染仃囆纬煽倓偠染仃嚕ǚQ(chēng)為總體合成）。(a) 把 按單元組織，則在第i個(gè)單元上，令)(huj)(,)(1,)(, 1)(1, 1)(iiiiiiiiiiiiiaaaak)()(1)(1)()(1,)(, 1iiiiiiiiiiiiiifffuuuaa其中 稱(chēng)為單元?jiǎng)偠染仃?。各元素可?jì)算得到。

58、)(ik117)(,)(1,)(, 1)(1, 1)(iiiiiiiiiiiiiaaaak 再把 擴(kuò)展成nn矩陣，使其第i1行、第i行和第i1列、第i列交叉位置的元素就是單元?jiǎng)偠染仃嚨乃膫€(gè)元素，其余全為零（ 只是第一行，第一列元素非零）。即)(ik)1(ktntnbbbbuuuu),.,(),.,(2121記則),(),(2121)(ubuukbuukuujttnniikk1)(其中 稱(chēng)為總剛矩陣。118(b) 由 達(dá)到極小值的條件)(hujnjuujjh,.,2 , 10)(),.,1(niui(c) 解出有限元方程的數(shù)值解 ，就得到使二次泛函取極小的近似函數(shù)(有限元解)()(1xuxuin

59、iih得到。buk119把表達(dá)式 代入變分方程niiihxuu1)(對(duì)前面的兩點(diǎn)邊值問(wèn)題，變分方程變?yōu)閚jdxfubajiniji,.,2 , 1),(1dxquvvupvuba)(),(其中該方程即為galerkin法形成的。120 在單元 上可構(gòu)造一、二、三及高次插值多項(xiàng)式，其：,1iiixxi ：在單元內(nèi)部增加一些插值節(jié)點(diǎn)。：在節(jié)點(diǎn)引進(jìn)一階、二階乃至更高階導(dǎo)數(shù)。 121：在每一個(gè)單元上是一次多項(xiàng)式，在單元節(jié)點(diǎn)處連續(xù)。：在單元的兩個(gè)端點(diǎn)取指定值。：在每一個(gè)單元上是二次多項(xiàng)式，在單元節(jié)點(diǎn)處連續(xù)。：在單元的兩個(gè)端點(diǎn)及單元中點(diǎn)取指定值。：在每一個(gè)單元上是三次多項(xiàng)式，在單元節(jié)點(diǎn)處連續(xù)。：在兩個(gè)端點(diǎn)

60、取指定的函數(shù)值和一階導(dǎo)數(shù)值。122 采用高次元，有限元方程形成的方法和線(xiàn)性元類(lèi)似，但工作量增加。一是計(jì)算積分的復(fù)雜性增加，二是矩陣的帶寬增加。 高次元的主要優(yōu)點(diǎn)是收斂階高，且提高了函數(shù)逼近的光滑性。123 假定區(qū)域g可以分割成有限個(gè)矩形的和，且每個(gè)小矩形（單元）的邊和坐標(biāo)軸平行。 1 , 0; 1 , 0iiyyyxxxii/ )(/ )(;yyyyxxxxrjjiiij通過(guò)仿射變換采用矩形剖分后，任一個(gè)矩形總可變成單位正方形 如果在 上造出單元形狀函數(shù)，就可得到試探函數(shù) 。而 上的形狀函數(shù)可通過(guò)先在 上造出形狀函數(shù)，再通過(guò)仿射變化而得到。iiijrijr)(xuh124ii 在上構(gòu)造形狀函數(shù)，也采用