偏微分方程數(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. 李榮華，馮國忱. 微分方程數(shù)值解. 北京：人民教育

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

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

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方程截斷為五維的截譜模型如下：112345

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

7、化為代數(shù)方程求解，一般常用遞推算法。 在常微分方程差分法中最簡單的方法是euler方法，盡管在計算中不會使用，但從中可領(lǐng)悟到建立差分格式的技術(shù)路線，下面將對其作詳細(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考慮如下兩個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對經(jīng)典的初值問題0( , )(0)duf t udtuu(0, )tt1212( ,)( ,)f t uf t ul uu滿足lipschitz條件保證了方程組的初值問題有唯一解。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)用時采用如下遞推方式計算：0u1u2u3u0t1t2t4. 例題對初值問題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過點(diǎn)( ,)iit u作 的切線，該切線的( )u t方程為：11( )()iiiiiuuu ttt即：1( ,)iiiiuuhf t u20二、誤差分析構(gòu)造算法后，這一算法在實(shí)際中是否可行呢？也就是說是否使計算機(jī)仿真而不失真，這還需要進(jìn)一步分析。1. 局部截斷誤差-相容性為了分析分析數(shù)值方法的精確度，常常在( )iiuu t成立的假定下，估計誤差111()iiieu tu這種誤差稱為“局部截斷誤差”，如圖。21()( )( )( )2iiihu tu thu tu2( )( , ( )( )2iiihu thf t u tu22111()( )()2

11、iiiheu tuuo h局部截斷誤差是以點(diǎn)it的精確解( )iu t為出發(fā)值，用數(shù)值方法推進(jìn)到下一個點(diǎn)1it而產(chǎn)生的誤差。21整體截斷誤差是以點(diǎn)0t的初始值0u為出發(fā)值，用數(shù)值方法推進(jìn)i+1步到點(diǎn)1it，所得的近似值 與精確值 的偏差：2.整體截斷誤差收斂性1iu1()iu t111()iiiu tu稱為整體截斷誤差。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若不計初始誤差，即則1(1)tlidhel1( )io h即3.舍入誤差穩(wěn)定性假設(shè)一個計算機(jī)僅表示4個數(shù)字（小數(shù)點(diǎn)后面）， 那么10.3333,310.11119計算20.373410.1112011119110.33340.33333xxx20.33341190.6667133xxxx23我們的要求是：最初產(chǎn)生的小誤差在以后的計算中雖然會傳遞下去，但不會無限制的擴(kuò)大，這就是穩(wěn)定性所描述的問題。下面引進(jìn)穩(wěn)定性的概念：設(shè)由初值0u0v得到精確解 ，nu由初值得到精確解 ，nv若存在常數(shù)c和充分小的步長0h使得00nnuvc uv則稱數(shù)

13、值方法是穩(wěn)定的。tu00u0vntnunv242(0)1dyxydxyy計算例題(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用梯形法計算積分的

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這是一個隱式格式，一般需要用迭代法來求 ，而用顯式的euler法提供初值。1iu2701( ,)iiiiuuhf t u1 111/2 ( ,)(,)kkiiiiiiuuhf t uf tu為了簡化計算的過程，在此基礎(chǔ)上進(jìn)一步變?yōu)槿缦滤惴ǎ?11/2 ( ,)(,)iiiiiiuuhf t uf tu1( ,

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

16、ler法是以兩項taylor級數(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級數(shù)的更多項會得到什么樣的求解方法呢？兩個德國數(shù)學(xué)家（c.runge & m.kutta）以這種思想為基礎(chǔ)建立了求解微分方程的龍格-庫塔方法。它是常微分方程數(shù)值解法中使用最為廣泛的方法之一。31一般地，一個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展開成h的冪級數(shù)，并與微分方程的精確解1()iu t在點(diǎn)it的taylor展開式相比較，使兩者的前1p項相同，這樣確定的r-k法，其局部截斷誤差為1()po h，根據(jù)所得關(guān)于待定系數(shù)的方程組，求出它們的值后代入公式，就成為一個p階r-k方法。32例題以二階r-k法為例說明上述過程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方法，對應(yīng)的有中點(diǎn)法、heun（亨）法以及改進(jìn)的euler法。34基于相同的過程，通過比較五次taylor多項式，得到更加復(fù)雜的結(jié)果，給出了包含13個未知數(shù)的11個方程。得到多組系數(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五、線性多步（linear multistep method）法1. 預(yù)備知識：插值多項式插值是離散函數(shù)逼近的重要方法，利用它可通過函數(shù)在有限個點(diǎn)處的取值狀況，估算出函數(shù)在其他點(diǎn)處的近似值。從幾何上理解：對一維而言，已知平面上n1個不同點(diǎn)，要尋找一條n次多項式曲線通過這些點(diǎn)。插值多項式一般常見

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

22、-bashforth（阿達(dá)姆斯貝雪福斯）公式首先，用以下四個點(diǎn)對( , ( )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基于同樣的計算過程，可以得到另外一個計算公式：11111229 (,) 19 (,)5 (,)(,)24nnnnnnnnnnhuuf xuf x uf xuf xu這稱為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初值問題)0,(txd)(4)(exp41),(2ftxttxu不論初始分布如何集中，它總在瞬間影響于無窮遠(yuǎn)，雖該影響隨距離按指數(shù)衰減，然而它是以無限速度傳播。此乃拋物型方程解的特征。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計算：111(1 2 )jjjjiiiiusus usu這是一個顯式格式（四點(diǎn)格式）j1j i1i1

27、i每一層各個節(jié)點(diǎn)上的值是通過一個方程組求解得到的。這可以從下面的計算過程看出來。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 計算實(shí)例：4800.020.040.060.080.100.51-3-2-10123txu492. 向后差分格式111121122jjjjjiiiiiuuuuuakh11111(1

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

29、.020.040.060.080.100.50.81txu533. crank-nicolson格式，亦稱六點(diǎn)對稱格式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這是一個五點(diǎn)三層差分顯式格式57討論：假若由于 的作用，導(dǎo)致差分方程的近似解設(shè)為：jiu 于是，我們可得到差分格式的誤差方程如下：11112 (2)jjjjjiiiiisxtrichardson格式是不穩(wěn)定的。585. 穩(wěn)定性判別 von-neumann 穩(wěn)定性在判斷有限差分近似的穩(wěn)定性方法中，以von-neumann方法使用較為廣泛，它僅適用于線性常系數(shù)的有限差分近似。其過程如下：首先，要研究的差分方程可寫為：101nnmj mmj mm nm na u

31、b u如，其次，對進(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 下面說明，在什么條件下能使11nnvv對所有的成立。21 1 4 sin12hs 從上式，我們看出，21 1 4 sin

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

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

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

35、 前兩個格式的局部截斷誤差階為 ，分別稱為左、右偏心格式。)(ho)(2ho 第三個格式的截斷誤差階為 ，它稱為中心差分格式。kjkjkjurarauu)1 (11kjkjkjrauurau11)1 ()(2111kjkjkjkjuuaruu其中hr/即67 從差分格式依賴區(qū)域和微分方程依賴區(qū)域的關(guān)系，可以得到差分格式收斂的必要條件： 對于左偏心格式， cfl條件為：a0，且 。1ar1ar對于右偏心格式， 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)作特征線，它與線段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)定，截斷誤差階為 。1ar)(22hho69 對方程(3.1)，利用特征線作二次插值，即可得到laxwendroff格式： 當(dāng) ， laxwendroff格式穩(wěn)定，它的截斷誤差階為 。1ar)(22hokjxkjkjkjkjurauuaruu2221112)(270應(yīng)該注意：邊值條件的給法與其它兩類方程不同。 如果 a0，方程特征線向右傾，只能在 x變化區(qū)域的左邊界上給出邊界條件：)(), 0(1ttu

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

39、似值。00,1, 2,0, 1, 2,kkjjuujk75 222220uuaxatx 為正常數(shù)考察 對x-t平面進(jìn)行矩形網(wǎng)格剖分，x方向的步長為h，t方向的步長為，網(wǎng)點(diǎ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 上述差分格式與初始條件的截斷誤差階均為 。jju0取為77ttxxuatu0, 10022222上述方法也可用于求解初邊值問題： 10)()0 ,()()0 ,(xtxtuxxuttttuttu0)(), 1 ()(), 0(78 798081在條件下2122ssldsdsdydsdx求使得泛函達(dá)到最大的函數(shù) 。)(),(sysx2121),(ssdsdsdxydsdyxyxs 在長度一定的所有平面封閉曲線中，求所圍面積

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

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

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

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

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

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

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

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

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

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

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

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

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

54、， 的（線性）插值公式稱為（線性）單元形狀函數(shù)。iix)(xuh109 把每個單元形狀函數(shù)合并起來，就得到整個區(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 線性無關(guān)，它們可組成試探函數(shù)空間的基，常稱為節(jié)點(diǎn)基函數(shù)。幾何形狀如圖)(0 x)(xia)(xn0 x1, 1nib1x1ixix1ix1nxnx任一試探函數(shù) 可表示為hhvu niiihbaxxuxu0,)()( 用這類插值型基函數(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其中稱為總剛矩陣。nnnnaaaaaaak, 13222122111116 工程中形成有限元方程時，通常先在每個單元上形成單元矩陣（稱為單元剛度矩陣），然后由單元剛度矩陣形成總剛度矩陣（稱為總體合成）。(a) 把 按單元組織，則在第i個單元上，令)(huj)(,)(1,)(, 1)(1, 1)(iiiiiiiiiiiiiaaaak)()(1)(1)()(1,)(, 1iiiiiiiiiiiiiifffuuuaa其中 稱為單元剛度矩陣。各元素可計算得到。

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

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

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