第十一章常微分方程邊值問(wèn)題的數(shù)值解法

1、第十一章常微分方程邊值問(wèn)題的數(shù)值解法工程技術(shù)與科學(xué)實(shí)驗(yàn)屮提出的大量問(wèn)題是常微分方程邊值問(wèn)題.本章將研究常微分方程 邊值問(wèn)題的數(shù)值求解方法.主要介紹三種邊界條件下的定解問(wèn)題和兩大類求解邊值問(wèn)題的數(shù) 值方法，打靶法算法和有限差分方法.11.1引言在很多實(shí)際問(wèn)題屮都會(huì)遇到求解常微分方程邊值問(wèn)題.考慮如下形式的二階常微分方程, a < x < b，(11.1.1)在如下三種邊界條件下的定解問(wèn)題：第一種邊界條件：ya) = a, yb) = (3(11.1.2)第二種邊界條件:第三種邊界條件:/=0c，yb) = pya)-a)y(a) = a(11. 1.2)其中 6r() 2 0, z?

2、o 2 0，tz0 + /?()0.常微分方程邊值問(wèn)題有很多不同解法，本節(jié)僅介紹打靶方法和有限差分方法.11.2打靶法對(duì)于二階非線性邊值m題/ =a<x<b, y(a) = af y(b) = /?.(11.2. 1)打靶法近似于使用初值求解的情況.我們需要利用一個(gè)如下形式問(wèn)題初值解的序列:”二 /(x，w，vv,)，a < x <b, vv(6z) = a, wa) = v(11.2.2)引進(jìn)參數(shù)7以近似原邊界值問(wèn)題的解.選擇參數(shù) = ,以使：lim vv(6, vk) = y=(3，(11.2.3)其屮w（x，va.）定義為初值問(wèn)題（11.2.2）在v = va.時(shí)

3、的解，同時(shí)少（x）定義為邊值問(wèn)題 （11.2. 1）的解.首先定義參數(shù)沿著如下初值問(wèn)題解的曲線，可以求出點(diǎn)（人60對(duì)應(yīng)的初始正視圖 w* = /（x, w, a < x < b，vitz） = a, wa） = v（n 2 4）如果w（6,v（>）不嚴(yán)格收斂于/?，那么我們選擇v,等值以修正近似值，直到w<z?，v（）嚴(yán)格逼近/?為了取得合適的參數(shù)va.,現(xiàn)在假定邊值問(wèn)題（11.2.1）有唯一解，如果w（x，v）定義為初始問(wèn)題（11.2.2）的解，那么v可由下式確定：(11.2.5)vv（/?, v） - /? = 0由于這是一個(gè)非線性方程，我們可以利用newton法求

4、解.首先選擇初始值,然后由下式生成序列(11.2.6)(字)udvv, = v,.,-(，-)_)，此處今)(,) =宇dw.dvdv同吋要求求得（）（6, vg），因?yàn)関v（6，v）的表達(dá)式未知，所以求解這個(gè)有一點(diǎn)難度dv我們只能得到這么一系列的值 ma vq），w（/?, v,）, w（z?, v2）,w（z?, va._, ）假如我們?nèi)缦赂膶?xiě)初值問(wèn)題（11.2.2）,使其強(qiáng)調(diào)解對(duì)x和v的依賴性= /（又，w（jv，v），vv（x，v），a < x <b，w（“，v） = 66 w（“，v） = v，（11.2.7）保留初始記號(hào)以顯式與*的微分相關(guān).既然要求當(dāng)v = 時(shí)（+）（

5、6，v）的值，那么我們dv耑要求出表達(dá)式（11.2.7）關(guān)于v的偏導(dǎo)數(shù).過(guò)程如下：（x，v）=拿（x，w（x，v），w （x, v） dvdv= -（x，w（x，v），wz（x, v）字 + 孚-（x，w（x，v），v）（x, v）oxdv owdv+（%，w（x，v），w（x，v）（x，v）owdv又因?yàn)?跟v相互獨(dú)立，所以當(dāng)上式如下；dw z df .,z 、3wdfv) = 7-(x, w(>，v)，w(x，v) 了 (x，v) + 40, w(x，v)，wx, v)(x, v)dvawovawav(11.2.8)初始條件為.3w-(6z,v)=0,-(v)= l avdv如果簡(jiǎn)單

6、地用z(x，v)定義(¥)(x，v)，并且假定x和v的微分階翻轉(zhuǎn)，(11.2.8)轉(zhuǎn)化 dv為初值fuj題:z = (x，vv，vv)z + 7；(x，w，vv)z，a < x <b，z(a) = 0，z(“)= 1dvvdw，(11.2.9)因此，對(duì)于(11. 2. 2)和(11. 2. 9)式的每次迭代需要求解兩個(gè)初伉問(wèn)題.那么從(11. 2. 6)式可得:(11.2. 10)事實(shí)上，這些初值問(wèn)題很難精確求解，而將這些解近似為一個(gè)初值問(wèn)題的解.同樣，我們可以按以上步驟考慮對(duì)于三階非線性邊值問(wèn)題的打靶法算法.對(duì)于三階非線性邊值w題= f(x,y9yy a <x &

7、lt;b, yci) = a, ya) = p,yb) = /. (11.2.11)轉(zhuǎn)變形式：= ,(x，w，w，v/)，a < x < b，= a, wa) = /?, vz(7) = v, (11. 2. 12) 選擇參數(shù)v = va.,以使：lim vk) = yb) = (5，(11.2. 13)炎一>oo其中定義為初值問(wèn)題(11.2. 12)在v = va.時(shí)的解，同時(shí)yoo定義為邊值問(wèn)題 (11.2. 11)的解.定義參數(shù)v(),沿著如下初值問(wèn)題解的曲線，可以求出點(diǎn)(人60對(duì)應(yīng)的初始正視圖 v/* =/(x，vv，w，vz)，a < x < b，w(a

8、) = a, wa) =w*(z) = v . (11.2.14)如果w'(z7,&)不嚴(yán)格收斂于/?,那么我們選擇&等值以修正近似值，直到>/(/),v(>)嚴(yán)為了取得合適的參數(shù)現(xiàn)在假定邊值問(wèn)題（11.2.11）有唯一解，如果vv（x，v）定義為初始問(wèn)題（11.2.12）的解，那么可由下式確定：w(6, v)-夕=0 (1l2. 15)利用newton法求解.首先選擇初始值,然后由下式生成序列= -1 -:廣-1 卜灼，此處($)(/?人-1) = $ (么 w-1),(1l2. 16)dv同時(shí)要求求得因?yàn)?/（6,v）的表達(dá)式未知，所以只能得到這么一系列

9、dv的值 w（z?，v0），wb, v,）, wb, v2）,wb, vk_x ）o改寫(xiě)初值問(wèn)題（11.2.12）,使其強(qiáng)調(diào)解對(duì)x和v的依賴性,(x，w(x，v)，wx, v),v),a < x < b，w(6r，v) = 66 wa,v) = (5, w"(o，v) = vdv(11.2. 17)要求當(dāng)v =的值，那么我們耑要求出表達(dá)式（11.2.17）關(guān)于v的偏導(dǎo)dv數(shù).所以當(dāng)a < x s z?上成如下;dv3vv dv dw dv dw dv(11.2. 18)初始條件為:3vvdwo,-(t7,v) = 0,(tz,v) = ldvdvdv對(duì)于第三種邊界條

10、件，注意到y(tǒng)（a）=，取s得到y(tǒng) = .f(xfyy)，/(11.2. 19)沒(méi)它的解為奴x，4,代入第二個(gè)邊界條件得到y(tǒng)b) - fioy(b) = y(m-= f二凡當(dāng)汊=從卩、即力所求的解.這樣，同第一種邊界條件一樣，可以利用打靶方法求解.例11.1利用打靶法求解兩點(diǎn)邊值問(wèn)題/-(/)2 + 又=。v(0) = 0，x0.4) = 1.8解：a)為應(yīng)用打靶方法，需要假定初值/(o)來(lái)先求解初值問(wèn)題，取初始參數(shù)0.4-0j；(0) = 0,/(0) = 2令少=%，/ = w2，則上述初值問(wèn)題變?yōu)?，由?biāo)準(zhǔn)的r-k方法，取a = 0.2,可得沐=2.540取 a = 0.2ux =w2，w1

11、(0) = l u2 = xt> u2 (0) = 2(2)再令z, =1，解ux =w”w1(0) = l = u x, u2 (0) = 1仍由標(biāo)準(zhǔn)的k-k方法，可得我=1.497<1.8(3)令z2 = 1 +1-2再解1.497-2.540 w/ = u2,u (0) = 1u2 -ux,u2(0) = 1 >29(1.8-1.497) = 1.29,取 a = 0.11,得為=1.7094(4)令z3 = 1.29 +(1.29-1)(1.8-1.7094)1.7094-1.4971.414,得a =1.81978(5)類似有，z4= 1.38,a = 1.7887

12、9, z5 = 1.3960, = 1.79995因?yàn)閨#5 -1.8| < 0.00005 ,可以認(rèn)為z = 1.3960即為所求。11.3有限差分方法有限差分方法常用來(lái)求解常微分方程邊伉問(wèn)題，下而就線性微分方程為例進(jìn)行討論.設(shè) 微分方程具有形式:y(a) - a, y(b) = (i.(11. 3. 2)不失一般性，將求解區(qū)間a?等份，分點(diǎn)為x0=6f，xj =a + h , , xz = tz + ih h ci(/ = 0, 1，n),步長(zhǎng)/7 = _,這里xz稱為結(jié)點(diǎn). n在內(nèi)部結(jié)點(diǎn) x,上，記 yxi) = yi p(xi) = pi, q(xi) = qi, f (x,.)

13、 = /,并分別用一階和二階中心差分來(lái)代替一階導(dǎo)數(shù)和二階導(dǎo)數(shù)，八+ 2)，2/7，('.) =”公 + +o(h2).將上面各式代入方程(11. 3.1),略去0(/72)項(xiàng)，得到關(guān)系式：0；+i 2乃 +z-i)+ a(z+i -z-i)+ /?2 =/?2z-將上式改寫(xiě)成19辦'z-i + c'z+i =一5,7 乂，(, =1，”-1)，(11. 3. 3)在上式中/z2l h1 hai =1-了厶,=2(2幾)，c' =2(1+2凡) (z=1 ，w_l)，(1l3.4) 利用邊界條件y(a) = a, y(b) = f可以得到如下形式的方程組：ay =

14、 f其中(11.3.5)a2c2a3(11.3.6)an-2f'b'a/20_fn- _+0_cn-p_(11.3.7)系數(shù)矩陣(11. 3.6)為n-1階三對(duì)角陣，當(dāng)系數(shù)矩陣非奇異時(shí)，方程組(11. 3. 5)有唯一解, 可以利用追趕法求解.對(duì)于第二種邊界條件：(11.3.8)/() = 0c, yb) = p和笫三種邊界條件：(11. 3. 9)/-aqy(a) = a, yb)-/30y(b) = b j其中a。20, bq >0,a0 + bq > 0.由于邊界條件中包含了導(dǎo)數(shù)，故邊界條件也必須用差商近似來(lái)表示.因?yàn)椴荒芾脜^(qū)間k,川以外的點(diǎn)，所以在引進(jìn)(即

15、/或凡(即/(/?)的近似式時(shí)，就不能利用屮心差商.如果只要求誤差是o(a)，則可以用最簡(jiǎn)單的近似公式vo乃1h(11.3. 10)要使誤差達(dá)到o(a2),可以利用newton等距插值公式,得到如下近似公式1h-全心。一滅2 + 4; 3y02h廠4凡-i +凡-22h(11.3. 11)(11. 3. 12)將這些近似式代入邊值條件中，再和方程(11. 3. 4)聯(lián)立，就可以得到對(duì)應(yīng)的差分方程組. 例11.2利用差分法求解兩點(diǎn)邊值問(wèn)題p"-(l + x2)產(chǎn) 1 y(0) = l，xl) = 3解：將方程離散后寫(xiě)成矩陣形式'-2.0101 1 '火2"-0

16、.99'0.011一 2.010411 -2.016410.011 -2.0181_入-2.99®a = 0.05,可以得到，凡三對(duì)角方程，下表11-1給fli了計(jì)算結(jié)果。表11h収a = 0.05吋計(jì)算結(jié)果xa = 0.1,近似值兄，a = 0.05近似值凡少=精確值0. 11. 07461. 07441. 07430.21. 17011. 16961. 16950.31. 28771. 28711. 28690.41. 42941. 42861. 42840. 51. 59761. 59681. 59650.61. 79581. 79501. 79470.72. 0285

17、2. 02762. 02740.82. 30132. 30062. 30040.92. 62192. 62152. 621411.4應(yīng)用實(shí)例與mat lab11.4.1 matlab關(guān)于常微分方程邊值問(wèn)題數(shù)值解法的命令 常微分方程初、邊值問(wèn)題的符號(hào)解法命令形式 1: dsolveg eqution,）功能：求常微分方程eqution的解。命令形式 2: dsolve（, eqution, condl, cond2，var，）功能：求常微分方程eqution的滿足初始條件的特解。其屮eqution足求解的微分方程或微分方程組，condl, cond2是初始條件或邊值，var是變量.例11. 3求

18、解兩點(diǎn)邊值|'uj題：ixi) = o y(5) = 0解mat lab命令為y=dsolve(4x*d2y-3*dy=xa2，y(l)=0, y(5)=0，x)執(zhí)行結(jié)果：y = 1 / 3 * xa3+125 / 468+31 / 468 * xa411.4.2應(yīng)用實(shí)例當(dāng)流體以均勻的流速縱掠一平壁時(shí)，巾于流體粘性的影響，在靠近殖而鄰近會(huì)形成很大 的速度梯度，這就是速度邊界層.用來(lái)描述速度邊界層的動(dòng)量方程經(jīng)過(guò)相似變換付以化為如 下常微分方程兩點(diǎn)邊界值m題(著名的布拉修斯邊界層方程)(11. 3. 13)邊界條件為(11.3. 14)/(0) = ±7v，/(0) = 0，/(

19、+oo) = 0這里/(0)表示壁面噴注(或吸入).根據(jù)給定的噴注速度值，/(0)二+au常數(shù))表示吸入，/(0) = -n (常數(shù))表示噴注.為了利用打靶法求解該問(wèn)題的解，假設(shè)(這里0是一個(gè)待定常數(shù)，需要在 計(jì)算屮估計(jì)).表11-1給出了數(shù)值計(jì)©方法。當(dāng)7 = 0吋取=、/(0>0和= 這里夕是一個(gè)常數(shù)(估計(jì)值)，這樣就能根據(jù)式(11.3. 13) - (11.3.14),求得接著，令7具有某一增量，分別汁算/、/、/"和/",繼續(xù)進(jìn)行下去一直到邊 界外邊緣(在平壁的情況下它相當(dāng)于77 = 5.0).此時(shí)，應(yīng)該得到/卜>1.如果這個(gè)要求 不能滿足，就

20、要修改原來(lái)的估計(jì)值，直至使/= l這個(gè)邊界條件得到滿足為止.，表u-2給山了數(shù)值計(jì)算方法表11-2布拉修斯a程式的數(shù)值計(jì)算a法7fffr0n0式(1.31)0. 1n+0. 1*00+0. 1氺 p夕十 0.1*/"(0)式(1.31)0.2n+0. 1*0+0. 1*(0+0. 1*#)0+0. 1*/?+0. 1*/?+() 1*/"(0)邊1*廣(0)+0. 1*.廣(0)式(1.31)攀礬拳攀拳拳拳拳攀50.99參番番番番oo1.0?在程序中為了能夠方便的變換#"系數(shù)，我們把方程中/廣系數(shù)設(shè)為epsilo, 從而把方程改寫(xiě)為：廣 + epsilo*/

21、9;*/'/r=0(11.3.13)mat lab 命令:?itblasius (0, 0, 0, 0. 5, 0. 00001, 0. 1, 5, 0. 001) beta=0. 435860得出每一步的迭代結(jié)果：matlab 命令:?blasius(0, 0, 0. 435860, 0. 1,5, 0. 001, 1, 1, 1, 1, 1)beta=0. 435860表11-2給出了當(dāng);v = 0情況在區(qū)間0, 5的部分計(jì)算結(jié)果.mat lab程序function itaconv=blasius (n, m，beta, dita, max it a, tol，options)fo

22、rmat long;%(f"+epsilo*f*f=0)%書(shū)上是epsilo=0. 5%epsilo=l:epsilo=0. 5;%if options(5)=1fprintf ( nbeta=%fnttitattfttf ttf ttf ， nn，beta); endfoundo;xo=n;yo=m;zo=bcta;wo=-xo*zo*epsilo;for ita=0 : dita : maxita switch options(1) case 1multiplicr=o. 1;% 推薦 otherwisefprintf (n出錯(cuò)了！ nn); breakendswitch opt

23、ions(2) case 1x=xo+mu ltipli er*y 0; %推薦 otherwisefprintf (n出錯(cuò)了！ nn ); breakendswitch options(3) case 1y=yo*mu ltipli er*zo; % 推薦 otherwisefprintf(11出錯(cuò)了! nn)； breakendswitch options(4) case 1z=zo+mu ltipli er*wo; %推薦 otherwisefprintf ( n出錯(cuò)了! nn); breakendw=-x*z*epsilo;%if options(5)=1fprintf ( t%2.

24、5ft%2. 5ft%2. 5ft%2. 5ft%2. 5fn，ita, xo, yo, zo, wo); end%if abs(l-y)tol found=l; itaconv=ita; breakendx0=x;zo-z;w0=w;endaif"founditaconv=-l;endfunction itblasius (n, in, minbeta, maxbeta, dbeta, dita, maxita, tol)found二0;fprintf ( n白拉修斯方程數(shù)值計(jì)算迭代程序nn正在計(jì)算.nn)；for beta=minbeta:dbeta:maxbetaita=bla

25、sius(n, m, beta，dita, maxita, tol, 1，1, 1, 1, 0); if ita>-lfprintf ( 11滿足白拉修斯收斂條件的beta已經(jīng)找到，beta %fnn , beta);found=found+l;break%若需找到所有beta,請(qǐng)注釋掉上一行的break endendif found=0fprintf (* n在給定的范圍(%f. .%f)內(nèi)無(wú)法找到滿白拉修斯收斂條件的 beta. nn，minbeta, maxbeta): end下面表11-3給出了當(dāng);v = 0情況在區(qū)間0, 5的部分計(jì)算結(jié)果 表11-3:當(dāng);v = 0部分計(jì)算結(jié)果

26、.0000.332060.20. 006640. 066410. 331990.40. 026560. 132770. 331470.60. 059740. 198940. 330080.80. 106110. 264710. 3273910. 165570. 329790. 323011.20. 237950. 395780.316591.40. 322980. 456270. 307871.60. 420320.516760. 296671.80. 529520. 574770. 2829320. 650030. 629770. 266752.20. 78120. 681320. 2483

27、52.40. 92230. 728990. 228092.61.072520. 772460. 206462.81.230990.811520. 1840131.396820. 846050.161363.21.569110. 876090. 139133.41.746960. 901770. 117883.61.929540. 923330. 098093.82. 116050. 941120. 0801342. 305760. 955520. 064244.22.498060. 966960. 050524.42. 692380. 975870. 038974.62.888260. 982690. 029484.83. 085340. 987790. 0218753. 283290. 991550.01591小結(jié)本章研宂了常微分方程邊值問(wèn)題的數(shù)值求解方法.主要介紹了三種邊界條件下的定解 w題，通常主要有兩大類求解a法，即:打靶法算法和有限差分方法.打靶算法對(duì)于線性或非線性邊界值問(wèn)題均適應(yīng).邊值問(wèn)題的解可以視為一個(gè)“彈道曲 線”，打靶算法是通過(guò)引進(jìn)初值參數(shù)以近似原邊界值問(wèn)題，即將原微