常微分方程數(shù)值解法

常微分方程數(shù)值解法第一頁(yè)，共四十八頁(yè)，2022年，8月28日所謂數(shù)值解法,就是設(shè)法將常微分方程離散化,建立差分方程,給出解在一些離散點(diǎn)上的近似值.

a=x0<x1<x2<…<xn<…<xN=b其中剖分節(jié)點(diǎn)xn=a+nh,n=0,1,…,N,h稱為剖分步長(zhǎng).數(shù)值解法就是求精確解y(x)在剖分節(jié)點(diǎn)xn上的近似值yny(xn),n=1,2,…,N.假設(shè)初值問(wèn)題(8.1)的解y=y(x)唯一存在且足夠光滑.對(duì)求解區(qū)域[a,b]做剖分我們采用數(shù)值積分方法來(lái)建立差分公式.

§1.2構(gòu)造數(shù)值解法的基本思想在區(qū)間[xn,xn+1]上對(duì)方程(8.1)做積分,則有第二頁(yè)，共四十八頁(yè)，2022年，8月28日對(duì)右邊的積分應(yīng)用左矩形公式，則有第三頁(yè)，共四十八頁(yè)，2022年，8月28日梯形公式oxyab左矩形公式y(tǒng)=(x)右矩形公式中矩形公式第四頁(yè)，共四十八頁(yè)，2022年，8月28日對(duì)右邊的積分應(yīng)用左矩形公式，則有因此，建立節(jié)點(diǎn)處近似值yn滿足的差分公式稱之為Euler公式.稱為梯形公式.若對(duì)(8.2)式右邊的積分應(yīng)用梯形求積公式，則可導(dǎo)出差分公式第五頁(yè)，共四十八頁(yè)，2022年，8月28日利用Euler方法求初值問(wèn)題

解此時(shí)的Euler公式為稱為Euler中點(diǎn)公式或稱雙步Euler公式.若在區(qū)間[xn-1,xn+1]上對(duì)方程(8.1)做積分,則有對(duì)右邊的積分應(yīng)用中矩形求積公式，則得差分公式例1的數(shù)值解.此問(wèn)題的精確解是y(x)=x/(1+x2).第六頁(yè)，共四十八頁(yè)，2022年，8月28日分別取步長(zhǎng)h=0.2,0.1,0.05,計(jì)算結(jié)果如下第七頁(yè)，共四十八頁(yè)，2022年，8月28日hxnyny(xn)y(xn)-ynh=0.20.000.400.801.201.602.000.000000.376310.542280.527090.466320.406820.000000.344830.487800.491800.449440.400000.00000-0.03148-0.05448-0.03529-0.01689-0.00682h=0.10.000.400.801.201.602.000.000000.360850.513710.509610.458720.404190.000000.344830.487800.491800.449440.400000.00000-0.01603-0.02590-0.01781-0.00928-0.00419h=0.050.000.400.801.201.602.000.000000.352870.500490.500730.454250.402270.000000.344830.487800.491800.449440.400000.00000-0.00804-0.01268-0.00892-0.00481-0.00227第八頁(yè)，共四十八頁(yè)，2022年，8月28日Euler中點(diǎn)公式則不然,計(jì)算yn+1時(shí)需用到前兩步的值yn,yn-1,稱其為兩步方法,兩步以上的方法統(tǒng)稱為多步法.在Euler公式和梯形公式中,為求得yn+1,只需用到前一步的值yn,這種差分方法稱為單步法,這是一種自開(kāi)始方法.隱式公式中,每次計(jì)算yn+1都需解方程,要比顯式公式需要更多的計(jì)算量,但其計(jì)算穩(wěn)定性較好.在Euler公式和Euler中點(diǎn)公式中,需要計(jì)算的yn+1已被顯式表示出來(lái),稱這類差分公式為顯式公式,而梯形公式中,需要計(jì)算的yn+1隱含在等式兩側(cè),稱其為隱式公式.第九頁(yè)，共四十八頁(yè)，2022年，8月28日從數(shù)值積分的角度來(lái)看,梯形公式計(jì)算數(shù)值解的精度要比Euler公式好,但它屬于隱式公式,不便于計(jì)算.實(shí)際上,常將Euler公式與梯形公式結(jié)合使用:§2改進(jìn)的Euler方法和Taylor展開(kāi)方法

§2.1改進(jìn)的Euler方法第十頁(yè)，共四十八頁(yè)，2022年，8月28日

由迭代法收斂的角度看，當(dāng)

(是給定的精度要求)時(shí),取就可以保證迭代公式收斂,而當(dāng)h很小時(shí),收斂是很快的.而且,只要通常采用只迭代一次的算法:稱之為改進(jìn)的Euler方法.

這是一種單步顯式方法.改進(jìn)的Euler方法也可以寫(xiě)成第十一頁(yè)，共四十八頁(yè)，2022年，8月28日

y=y-2x/y,0x1的數(shù)值解,取步長(zhǎng)h=0.1.[精確解為y(x)=(1+2x)1/2.]例2

求初值問(wèn)題

y(0)=1

解

(1)利用Euler方法第十二頁(yè)，共四十八頁(yè)，2022年，8月28日計(jì)算結(jié)果如下:

(2)利用改進(jìn)Euler方法第十三頁(yè)，共四十八頁(yè)，2022年，8月28日nxnEuler方法yn改進(jìn)Euler法yn精確解y(xn)01234567891000.10.20.30.40.50.60.70.80.9111.11.1918181.2774381.3582131.4351331.5089661.5803381.6497831.7177791.78477011.0959091.1840961.2662011.3433601.4164021.4859561.5525151.6164761.6781681.73786911.0954451.1832161.2649911.3416411.4142141.4832401.5491931.6124521.6733201.732051第十四頁(yè)，共四十八頁(yè)，2022年，8月28日在節(jié)點(diǎn)xn+1的誤差y(xn+1)-yn+1,不僅與yn+1這一步計(jì)算有關(guān),而且與前n步計(jì)算值yn,yn-1,…,y1都有關(guān).為了簡(jiǎn)化誤差的分析,著重研究進(jìn)行一步計(jì)算時(shí)產(chǎn)生的誤差.即假設(shè)yn=y(xn),求誤差y(xn+1)-yn+1,這時(shí)的誤差稱為局部截?cái)嗾`差,它可以反映出差分公式的精度.§2.2差分公式的誤差分析如果單步差分公式的局部截?cái)嗾`差為O(hp+1),則稱該公式為p階方法.這里p為非負(fù)整數(shù).顯然,階數(shù)越高,方法的精度越高.研究差分公式階的重要手段是Taylor展開(kāi)式,一元函數(shù)和二元函數(shù)的Taylor展開(kāi)式為:第十五頁(yè)，共四十八頁(yè)，2022年，8月28日另外,在yn=y(xn)的條件下,考慮到y(tǒng)(x)=(x,y(x)),則有

y(xn)=(xn,y(xn))=(xn,yn)=n

y(xn)=第十六頁(yè)，共四十八頁(yè)，2022年，8月28日

yn+1=yn+h(xn,yn)對(duì)Euler方法，有

=yn+(xn,yn)h+O(h2)從而有:y(xn+1)-yn+1=O(h2)所以Euler方法是一階方法.再看改進(jìn)Euler方法,因?yàn)榭傻玫谑唔?yè)，共四十八頁(yè)，2022年，8月28日所以,改進(jìn)的Euler方法是二階方法.而從而有:y(xn+1)-yn+1=O(h3)§2.3Taylor展開(kāi)方法設(shè)y(x)是初值問(wèn)題(8.1)的精確解,利用Taylor展開(kāi)式可得第十八頁(yè)，共四十八頁(yè)，2022年，8月28日稱之為p階Taylor展開(kāi)方法.

……

……

……因此，可建立節(jié)點(diǎn)處近似值yn滿足的差分公式其中第十九頁(yè)，共四十八頁(yè)，2022年，8月28日所以,此差分公式是p階方法.由于Taylor展開(kāi)方法涉及很多復(fù)合函數(shù)(x,y(x))的導(dǎo)數(shù)的計(jì)算,比較繁瑣,因而很少直接使用,經(jīng)常用它為多步方法提供初始值.然而,Taylor展開(kāi)方法給出了一種構(gòu)造單步顯式高階方法的途徑.

Euler方法可寫(xiě)為可見(jiàn),公式的局部截?cái)嗾`差為:y(xn+1)-yn+1=O(hp+1).§3Runge-Kutta方法

§3.1Runge-Kutta方法的構(gòu)造第二十頁(yè)，共四十八頁(yè)，2022年，8月28日構(gòu)造差分公式

……

……改進(jìn)的Euler方法可寫(xiě)為其中i,i,ij為待定參數(shù).若此公式的局部截?cái)嗾`差為第二十一頁(yè)，共四十八頁(yè)，2022年，8月28日由于

yn+1=yn+h1n+h2(n+hxn+hnyn)+O(h3)O(hp+1),稱公式為p階Runge-kutta方法,簡(jiǎn)稱p階R-K方法.對(duì)于p=2的情形,應(yīng)有

=yn+h(1+2)n+h22(xn+nyn)+O(h3)所以,只要令1+2=1,2=1/2,2=1/2(8.4)第二十二頁(yè)，共四十八頁(yè)，2022年，8月28日一般地,參數(shù)由(8.4)確定的一族差分公式(8.3)統(tǒng)稱為二階R-K方法.稱之為中點(diǎn)公式,或可寫(xiě)為若取=1,則得1=2=1/2,=1,此時(shí)公式(8.3)就是改進(jìn)的Euler公式;若取1=0,則得2=1,==1/2,公式(8.3)為高階R-K公式可類似推導(dǎo).下面列出常用的三階、四階R-K公式.第二十三頁(yè)，共四十八頁(yè)，2022年，8月28日四階標(biāo)準(zhǔn)R-K公式

三階R-K公式第二十四頁(yè)，共四十八頁(yè)，2022年，8月28日

解四階標(biāo)準(zhǔn)R-K公式為例3

用四階標(biāo)準(zhǔn)R-K方法求初值問(wèn)題

y=y-2x/y,0x1

y(0)=1的數(shù)值解,取步長(zhǎng)h=0.2.計(jì)算結(jié)果如下:第二十五頁(yè)，共四十八頁(yè)，2022年，8月28日nxnyny(xn)nxnyny(xn)0120.00.20.41.001.18321.34171.001.18321.34163450.60.81.01.48331.61251.73211.48321.61251.7321也可以構(gòu)造隱式R-K方法,其一般形式為稱之為p級(jí)隱式R-K方法,同顯式R-K方法一樣確定參數(shù).如第二十六頁(yè)，共四十八頁(yè)，2022年，8月28日是二級(jí)二階隱式R-K方法,也就是梯形公式.但是p級(jí)隱式R-K方法的階可以大于p,例如,一級(jí)隱式中點(diǎn)公式為或?qū)憺樗嵌A方法.§3.2變步長(zhǎng)Runge-Kutta方法以p階R-K方法為例討論.設(shè)從xn以步長(zhǎng)h計(jì)算y(xn+1)的近似值為

,局部截?cái)嗾`差為其中,C是與h無(wú)關(guān)的常數(shù).第二十七頁(yè)，共四十八頁(yè)，2022年，8月28日如果將步長(zhǎng)減半,取h/2為步長(zhǎng),從xn經(jīng)兩步計(jì)算得到y(tǒng)(xn+1)的近似值記為

,其局部截?cái)嗾`差為于是有從而,得到事后誤差估計(jì)可見(jiàn),當(dāng)成立時(shí),可取

.否則,應(yīng)將步長(zhǎng)再次減半進(jìn)行計(jì)算.第二十八頁(yè)，共四十八頁(yè)，2022年，8月28日求解初值問(wèn)題的單步顯式方法可一統(tǒng)一寫(xiě)為如下形式

yn+1=yn+h(xn,yn,h)(8.5)對(duì)于Euler方法,有§4單步方法的收斂性和穩(wěn)定性

§4.1單步方法的收斂性

y=(x,y),axb

y(a)=其中(x,y,h)稱為增量函數(shù).(x,y,h)=(x,y)對(duì)于改進(jìn)的Euler方法,有(x,y,h)=1/2[(x,y)+(x+h,y+h(x,y))]第二十九頁(yè)，共四十八頁(yè)，2022年，8月28日設(shè)y(x)是初值問(wèn)題(8.1)的解,yn是單步法(8.5)產(chǎn)生的近似解.如果對(duì)任意固定的點(diǎn)xn,均有y(xn),則稱單步法(8.5)是收斂的.可見(jiàn),若方法(8.5)是收斂的,則當(dāng)h0時(shí),整體截?cái)嗾`差en=y(xn)-yn將趨于零.

定理8.1

設(shè)單步方法(8.5)是p1階方法,增量函數(shù)(x,y,h)在區(qū)域axb,-<y<+,0hh0上連續(xù),且關(guān)于y滿足Lipschitz條件,初始近似y0=y(a)=,則方法(8.5)是收斂的,且存在與h無(wú)關(guān)的常數(shù)C,使

|y(xn)-yn|Chp

證明因?yàn)閱尾椒椒?8.5)是p階方法,則y(x)滿足定義8.1第三十頁(yè)，共四十八頁(yè)，2022年，8月28日其中,局部截?cái)嗾`差|Rn(h)|C1hp+1,記en=y(xn)-yn,則有利用Lipschitz條件得

y(xn+1)=y(xn)+h(xn,y(xn),h)+Rn(h)遞推得到注意到

en+1=en+h[(xn,y(xn),h)-(xn,yn,h)]+Rn(h)

|en+1|(1+hL)|en|+C1hP+1

1+hLehL,(1+hL)nenhLeL(b-a)第三十一頁(yè)，共四十八頁(yè)，2022年，8月28日由于e0=y(a)-y0=0,所以有則有設(shè)(x,y)連續(xù)且關(guān)于y滿足Lipschitz條件,對(duì)于Euler方法,由于(x,y,h)=(x,y),故Euler方法是收斂的.對(duì)于改進(jìn)的Euler方法,由(x,y)的Lipschitz條件有

|en||e0|eL(b-a)+C1hp/L(eL(b-a)-1)

|en|C1hp/L(eL(b-a)-1)=Chp

|(x,y,h)-(x,y*,h)|1/2|(x,y)-(x,y*)|+1/2|(x+h,y+h(x,y))-(x+h,y*+h(x,y*))|1/2L(2+hL)|y-y*|則當(dāng)hh0時(shí),關(guān)于y滿足常數(shù)為1/2L(1+h0L)的Lipschitz第三十二頁(yè)，共四十八頁(yè)，2022年，8月28日條件,因此改進(jìn)的Euler方法是收斂的.可類似驗(yàn)證各階R-K方法是收斂的.§4.2單步方法的穩(wěn)定性

定義8.2

對(duì)于初值問(wèn)題(8.1),取定步長(zhǎng)h,用某個(gè)差分方法進(jìn)行計(jì)算時(shí),假設(shè)只在一個(gè)節(jié)點(diǎn)值yn上產(chǎn)生計(jì)算誤差,即計(jì)算值yn=yn+,如果這個(gè)誤差引起以后各節(jié)點(diǎn)值ym(m>n)的變化均不超過(guò),則稱此差分方法是絕對(duì)穩(wěn)定的.討論數(shù)值方法的穩(wěn)定性,通常僅限于典型的試驗(yàn)方程

y=y其中是復(fù)數(shù)且Re()<0.在復(fù)平面上,當(dāng)方法穩(wěn)定時(shí)要求變量h的取值范圍稱為方法的絕對(duì)穩(wěn)定域,它與實(shí)軸的交集稱為絕對(duì)穩(wěn)定區(qū)間.第三十三頁(yè)，共四十八頁(yè)，2022年，8月28日將Euler方法應(yīng)用于方程y=y,得到設(shè)在計(jì)算yn時(shí)產(chǎn)生誤差n,計(jì)算值yn=yn+n,則n將對(duì)以后各節(jié)點(diǎn)值計(jì)算產(chǎn)生影響.記ym=ym+m,mn,由上式可知誤差m滿足方程m=(1+h)m-1=…=(1+h)m-nn,mn對(duì)隱式單步方法也可類似討論.如將梯形公式用于方程y=y,則有

yn+1=yn+h/2(yn+yn+1)

yn+1=(1+h)yn

可見(jiàn),若要|m|<|n|,必須且只須|1+h|<1,因此Euler法的絕對(duì)穩(wěn)定域?yàn)閨1+h|<1,絕對(duì)穩(wěn)定區(qū)間是-2<Re()h<0.解出yn+1得

第三十四頁(yè)，共四十八頁(yè)，2022年，8月28日類似前面分析,可知絕對(duì)穩(wěn)定區(qū)域?yàn)橛捎赗e()<0,所以此不等式對(duì)任意步長(zhǎng)h恒成立,這是隱式公式的優(yōu)點(diǎn).一些常用方法的絕對(duì)穩(wěn)定區(qū)間為方法方法的階數(shù)穩(wěn)定區(qū)間Euler方法梯形方法改進(jìn)Euler方法二階R-K方法三階R-K方法四階R-K方法122234(-2,0)(-,0)(-2,0)(-2,0)(-2.51,0)(-2.78,0)第三十五頁(yè)，共四十八頁(yè)，2022年，8月28日

解因y0=1,計(jì)算得y10=1024,而y(1)=9.35762310-14.例4

考慮初值問(wèn)題

y=-30y,0x1

y(0)=1取步長(zhǎng)h=0.1,利用Euler方法計(jì)算y10y(1).[y(x)=e-30x]這是因?yàn)閔=-3不屬于Euler方法的絕對(duì)穩(wěn)定區(qū)間.若取h=0.01,計(jì)算得y100=3.23447710-16.若取h=0.001,計(jì)算得y1000=5.91199810-14.若取h=0.0001,計(jì)算得y10000=8.94505710-14.若取h=0.00001,計(jì)算得y100000=9.315610-14.第三十六頁(yè)，共四十八頁(yè)，2022年，8月28日單步顯式方法的穩(wěn)定性與步長(zhǎng)密切相關(guān),在一種步長(zhǎng)下是穩(wěn)定的差分公式,取大一點(diǎn)步長(zhǎng)就可能是不穩(wěn)定的.收斂性是反映差分公式本身的截?cái)嗾`差對(duì)數(shù)值解的影響;穩(wěn)定性是反映計(jì)算過(guò)程中舍入誤差對(duì)數(shù)值解的影響.只有即收斂又穩(wěn)定的差分公式才有實(shí)用價(jià)值.§5線性多步方法由于在計(jì)算yn+1時(shí),已經(jīng)知道yn,yn-1,…,及(xn,yn),(xn-1,yn-1),…,利用這些值構(gòu)造出精度高、計(jì)算量小的差分公式就是線性多步法.§5.1利用待定參數(shù)法構(gòu)造線性多步方法

r+1步線性多步方法的一般形式為第三十七頁(yè)，共四十八頁(yè)，2022年，8月28日當(dāng)-10時(shí),公式為隱式公式,反之為顯式公式.參數(shù)i,i的選擇原則是使方法的局部截?cái)嗾`差為

y(xn+1)-yn+1=O(h)r+2

選取參數(shù),0,1,2,使三步方法

yn+1=yn+h(0n+1n-1+2n-2)

這里,局部截?cái)嗾`差是指,在yn-i=y(xn-i),i=0,1,…,r的前提下,誤差y(xn+1)-yn+1.為三階方法.

例5

解設(shè)yn=y(xn),yn-1=y(xn-1),yn-2=y(xn-2),則有

第三十八頁(yè)，共四十八頁(yè)，2022年，8月28日n=(xn,y(xn))=y(xn)

y(xn+1)=y(xn)+hy(xn)+1/2h2y(xn)+1/6h3y(xn)于是有若使:y(xn+1)-yn+1=O(h4),只要,0,1,2滿足:n-1=(xn-1,y(xn-1))=y(xn-1)=y(xn-h)

=y(xn)-hy(xn)+1/2h2y(xn)-1/6h3y(4)(xn)+O(h4)n-2=y(xn)-2hy(xn)+2h2y(xn)-4/3h3y(4)(xn)+O(h4)

yn+1=y(xn)+h(0+1+2)y(xn)-h2(1+22)y(xn)

+h3(1/21+22)y(xn)-h4/6(1+82)y(4)(xn)+O(h5)

+1/24h4y(4)(xn)+O(h5)第三十九頁(yè)，共四十八頁(yè)，2022年，8月28日=1,0+1+2=1,1+22=-1/2,1+42=1/3于是有三步三階顯式差分公式設(shè)pr(x)是函數(shù)(x,y(x))的某個(gè)r次插值多項(xiàng)式,則有解之得:

yn+1=yn+h/12(23n-16n-1+5n-2)因?yàn)椤?.2利用數(shù)值積分構(gòu)造線性多步方法其中第四十頁(yè)，共四十八頁(yè)，2022年，8月28日選取不同的插值多項(xiàng)式pr(x),就可導(dǎo)出不同的差分公式.下面介紹常用的Adams公式.設(shè)已求得精確解y(x)在步長(zhǎng)為h的等距節(jié)點(diǎn)xn-r,…,xn上的近似值yn-r,…,yn,記k=(xk,yk),利用r+1個(gè)數(shù)據(jù)(xn-r,n-r),…,(xn,n)構(gòu)造r次Lagrange插值多項(xiàng)式由此,可建立差分公式

1.Adams顯式公式其中第四十一頁(yè)，共四十八頁(yè)，2022年，8月28日由此,可建立差分公式由于

hrj

則有稱之為r+1步Adams顯式公式.

第四十二頁(yè)，共四十八頁(yè)，2022年，8月28日下面列出幾個(gè)帶有局部截?cái)嗾`差主項(xiàng)的Adams顯式公式

r=0yn+1=yn+hn+(1/2)h2y(xn)

2.Adams隱式公式

r=1yn+1=yn+(h/2)(3n-n-1)+(5/12)h3y(xn)

r=2yn+1=yn+(h/12)(23n-16n-1+5n-2)+(3/8)h4y(4)(xn)

r=3yn+1=yn+(h/24)(55n-59n-1+37n-2-9n-3)+(251/720)h5y(5)(xn)如果利用r+1個(gè)數(shù)據(jù)(xn-r+1,n-r+1),…,(xn+1,n+1)構(gòu)造r次Lagrange插值多項(xiàng)式pr(x),則可導(dǎo)出數(shù)值穩(wěn)定性好的隱式公式,稱為Adams隱式公式,其一般形式為第四十三頁(yè)，共四十八頁(yè)，2022年，8月28日其中系數(shù)為下面列出幾個(gè)帶有局部截?cái)嗾`差主項(xiàng)的Adams隱式公式

r=0yn+1=yn+h