北京科技大學(xué)計(jì)算方法課件10第十章高斯求積公式及常微分方程初值問題

1、北京科技大學(xué)數(shù)理學(xué)院北京科技大學(xué)數(shù)理學(xué)院 衛(wèi)宏儒衛(wèi)宏儒W科學(xué)與工程計(jì)算科學(xué)與工程計(jì)算高斯求積公式高斯求積公式 常見的四個(gè)正交多項(xiàng)式常見的四個(gè)正交多項(xiàng)式P213-217勒讓德（勒讓德（Legendre）多項(xiàng)式）多項(xiàng)式 切比雪夫（切比雪夫（Chebyshev）多項(xiàng)式）多項(xiàng)式 拉蓋爾（拉蓋爾（Laguerre）多項(xiàng)式）多項(xiàng)式 埃爾米特（埃爾米特（Hermite）多項(xiàng)式）多項(xiàng)式 常微分方程初值問題的數(shù)值解法常微分方程初值問題的數(shù)值解法 1 1、基本概念和定理：、基本概念和定理：一階常微分方程初值問題是：一階常微分方程初值問題是： y =f(x,y) (1.1) y(x0)=y0 (1.2) 其中其中f

2、是已知的是已知的xoy平面上某個(gè)區(qū)域平面上某個(gè)區(qū)域D上連續(xù)函數(shù)，上連續(xù)函數(shù)，式（式（1.1）是微分方程，有無窮多解，式（）是微分方程，有無窮多解，式（1.2）是）是確定解的初始條件。確定解的初始條件。 如一元函數(shù)如一元函數(shù)y(x)對一切對一切axb b 滿足滿足（1 1） （x,y(x)Dx,y(x)D（2 2） y(xy(x0 0)=y)=y0 0（3 3） y y 存在，且存在，且y y (x)=f(x,y(x)(x)=f(x,y(x)則稱則稱y(x)y(x)是初值問題（是初值問題（1.11.1）、（）、（1.21.2）在）在a,ba,b上上的解。的解。 關(guān)于初值問題解的存在、唯一及對初始

3、條關(guān)于初值問題解的存在、唯一及對初始條件的連續(xù)依賴性，有下列定理：件的連續(xù)依賴性，有下列定理： 定理定理1:1: 設(shè)設(shè)f(x,y)f(x,y)是在是在D=(x,y)|D=(x,y)| axb, b, cycy d 上的連續(xù)函數(shù)，其中上的連續(xù)函數(shù)，其中a a，b b為有限實(shí)為有限實(shí)數(shù)，而且數(shù)，而且f(x,y)f(x,y)滿足對滿足對y y的的lipschitzlipschitz條件，條件，則對則對(x(x0 0,y,y0 0) ) D,D,初值問題（初值問題（1.11.1），（），（1.21.2）在在a,ba,b的解存在且唯一。的解存在且唯一。 若若y(x)y(x)是式（是式（1.11.1），（

4、），（1.21.2）的解，從方程）的解，從方程（1.11.1） 兩邊積分，再利用式（兩邊積分，再利用式（1.21.2） 可得可得積分積分方程方程 反之，若反之，若y(x) y(x) 滿足積分方程滿足積分方程 （1.41.4），可驗(yàn)證），可驗(yàn)證它滿足（它滿足（1.11.1） 和（和（1.21.2），所以），所以 （1.41.4）式與）式與初值問題（初值問題（1.11.1），（），（1.21.2） 等價(jià)等價(jià), ,這說明可用積這說明可用積分方程構(gòu)造初值問題的數(shù)值解法。分方程構(gòu)造初值問題的數(shù)值解法。 ).()(,()()(4100 xxdttytfxyxy定義定義3：若一種數(shù)值方法的：若一種數(shù)值方法的

5、局部截?cái)嗾`局部截?cái)嗾`差差O(hp+1)，則稱相應(yīng)數(shù)值方法是，則稱相應(yīng)數(shù)值方法是p階方階方法，其中法，其中p為正整數(shù)。為正整數(shù)。定義定義4：設(shè)：設(shè)y(x)是初始問題（是初始問題（1.1）的精）的精確解，確解，yn表示用某種數(shù)值方法算出的數(shù)表示用某種數(shù)值方法算出的數(shù)值解，值解， en= y(xn) - yn稱為該方法在稱為該方法在xn的的整體截?cái)嗾`差整體截?cái)嗾`差。 為了研究數(shù)值方法的絕對穩(wěn)定為了研究數(shù)值方法的絕對穩(wěn)定性，下面給出常系數(shù)線性差分方性，下面給出常系數(shù)線性差分方程的有關(guān)概念。程的有關(guān)概念。定義定義5 5： 方程方程111100n kkn knnyaya ya y 1110()0kknkr

6、ara ra r 1 12 2njnjnjnjk kyc rc rc r 0,1,jk 例例討論線性多步法的絕對穩(wěn)定性條件討論線性多步法的絕對穩(wěn)定性條件 得到相應(yīng)的齊次線性差分方程：得到相應(yīng)的齊次線性差分方程：00,0,0(,),.kkjkjjnjjjnnnjjyhfff xy 其中為常數(shù)不全為零將將應(yīng)用于實(shí)驗(yàn)方程應(yīng)用于實(shí)驗(yàn)方程yy 0()0kjjnjjhy 其對應(yīng)的特征方程為：其對應(yīng)的特征方程為：0()0kjjjjhr 0()0kjjnjjh 1knni iib r 2、數(shù)值解法的構(gòu)造途徑、數(shù)值解法的構(gòu)造途徑 （1）差商代替導(dǎo)數(shù)）差商代替導(dǎo)數(shù) 設(shè)初值問題設(shè)初值問題(1.1(1.1）的準(zhǔn)確解）

7、的準(zhǔn)確解y(x)y(x)在節(jié)在節(jié)點(diǎn)點(diǎn)x xn n之值為之值為y(xy(xn n),),記記y(xy(xn n) )的近似值為的近似值為y yn n ，又記，又記f fn n=f=f（x xn n ， y yn n ）, ,則初值問則初值問題題(1.1(1.1）離散化為）離散化為: : hxyhxyhxyhxyxykkkkhk)()()()(lim)( 0()()()(, ()(,)()()(,)0,1,2,.nnnnnnnnnnny xhy xy xf xy xf xyhy xhy xhf xyn即有得計(jì)算公式1(,)nnnnyyhf xy它稱為（向前）它稱為（向前）歐拉歐拉(Euler)公式

8、。（公式。（類類似地可以用向后差商、中心差商代替導(dǎo)似地可以用向后差商、中心差商代替導(dǎo)數(shù)產(chǎn)生相應(yīng)的歐拉數(shù)產(chǎn)生相應(yīng)的歐拉(Euler)公式公式） （2）數(shù)值積分法）數(shù)值積分法把把 y =f(x,y) 在在xn,xn+1積分，得積分，得對右端的定積分用數(shù)值積分方法做離對右端的定積分用數(shù)值積分方法做離散化，可得計(jì)算公式，如用散化，可得計(jì)算公式，如用矩形公式矩形公式可得歐拉公式，若用可得歐拉公式，若用梯形公式梯形公式可得改可得改進(jìn)的歐拉公式，它也稱為梯形公式：進(jìn)的歐拉公式，它也稱為梯形公式：11()()( , ( )nnxnnxy xy xf x y x dx111( (,)(,)2nnnnnnhyf

9、xyyf xy1111()()( , ( )( (, ()(, ()2nnxnnxnnnny xy xf x y x dxhf xy xf xy x （3）Taylor展開法 設(shè)設(shè) f(x,y) f(x,y) 充分光滑充分光滑, ,將將y(xy(xn+1n+1) )在在x x n n點(diǎn)作點(diǎn)作TaylorTaylor展展開開: : y (x y (x n+1 n+1)=y(x)=y(xn n)+hy)+hy (x(xn n)+(h)+(h2 2/2!)y/2!)y”(x xn n)+O(h)+O(h3 3) )取其關(guān)于取其關(guān)于h h 的的線性部分線性部分，并用，并用y yn n 代替代替 y (

10、x y (xn n),),就得就得到到EulerEuler公式公式。易知易知EulerEuler公式的局部截?cái)嗾`差為公式的局部截?cái)嗾`差為 T T1 1= (h= (h2 2/2!)y/2!)y ”(x(xn n)+O(h)+O(h3 3) = O(h) = O(h2 2) )改進(jìn)歐拉法的預(yù)改進(jìn)歐拉法的預(yù)- -校公式校公式1111(,)(,)(,)2nnnnnnnnnnyyhf xyhyyf xyf xy11(,)(,)2nnnnnnnnyyhf xyf xyhf xyEuler公式的幾何意義公式的幾何意義abY=y(x)xy0例題：用例題：用EulerEuler公式和改進(jìn)的公式和改進(jìn)的Eule

11、rEuler公式分別公式分別求下列初值問題的數(shù)值解求下列初值問題的數(shù)值解( (取步長取步長h=0.1h=0.1計(jì)算計(jì)算到到y(tǒng) y3 3) )： y y =-2xy =-2xy2 2 y(0)=1y(0)=1解：由歐拉公式解：由歐拉公式 y y n+1 n+1=y=y n n+hf(x+hf(xn n,y,y n n)=y)=y n n - 2hx - 2hxn ny y n n2 2計(jì)算如下計(jì)算如下y y 1 1 = y = y 0 0 - 2hx - 2hx0 0y y 0 02 2 =1-20.101=1-20.1012 2 =1=1y y 2 2 = y = y 1 1 - 2hx -

12、2hx1 1y y 1 12 2 =1-20.10.11=1-20.10.112 2 =0.98=0.98y y 3 3= y= y2 2 - 2hx - 2hx2 2y y 2 22 2 =0.98-20.10.20.98=0.98-20.10.20.982 2 =0.9416=0.9416用改進(jìn)歐拉法的預(yù)用改進(jìn)歐拉法的預(yù)-校公式計(jì)算如下：校公式計(jì)算如下： 21221112(22)2nnnnnnnnnnyyhx yhyyx yxy計(jì)算如下計(jì)算如下y 1 = 0.99; y 2 = 0.9614; y 3= 0.9173精確解精確解y(0.1)=0.99,y(0.2)=0.9614;y(0.3

13、)=0.9173 可見改進(jìn)歐拉公式比歐拉公式精度高?？梢姼倪M(jìn)歐拉公式比歐拉公式精度高。3、Runge-Kutta 方法 Runge-Kutta 方法是一種高精度的單步方法是一種高精度的單步法，簡稱法，簡稱R-K法。得到高精度方法的一個(gè)法。得到高精度方法的一個(gè)直接想法是利用直接想法是利用Taylor展開。展開。 假設(shè)式假設(shè)式 y =f(x,y) (axb) 中的中的 f(x,y) 充分光滑，將充分光滑，將y(xn+1)在在x n點(diǎn)作點(diǎn)作Taylor展開展開: （1）基本思想）基本思想對照標(biāo)準(zhǔn)形式對照標(biāo)準(zhǔn)形式 y n+1=y n+h(xn,y n,h) 。若若取取(x,y,h)=y(x)+(h/2

14、!)y(x)+.+(hp-1/p!)y(p)(x)并以并以y n代替代替y(xn)，則得到一個(gè)，則得到一個(gè)p階近似公式階近似公式 y n+1=y n+h(xn ,y n ,h) (n=0,1,2,.) (*) R-KR-K方法不是直接使用方法不是直接使用TaylorTaylor級數(shù)級數(shù), ,而是利用它的思想而是利用它的思想, ,即計(jì)算即計(jì)算f(x,y)f(x,y)在不在不同結(jié)點(diǎn)的函數(shù)值同結(jié)點(diǎn)的函數(shù)值, ,然后作這些函數(shù)值的然后作這些函數(shù)值的線性組合線性組合, ,構(gòu)造近似公式構(gòu)造近似公式, ,式中有一些式中有一些可 供 選 擇 的 參 數(shù) 。 將 近 似 公 式 與可 供 選 擇 的 參 數(shù) 。

15、 將 近 似 公 式 與TaylorTaylor展開式相比較展開式相比較, ,使前面的若干項(xiàng)使前面的若干項(xiàng)密合密合, ,從而使近似公式達(dá)到一定的精度。從而使近似公式達(dá)到一定的精度。 下面以二級二階下面以二級二階R-KR-K方法為例說明方法為例說明這一方法的基本思想。這一方法的基本思想。 在在 xn , , xn+1 上上, ,取取f(x,y)f(x,y)在兩個(gè)點(diǎn)的在兩個(gè)點(diǎn)的函數(shù)值作線性組合函數(shù)值作線性組合, ,即得到二級即得到二級R-KR-K方法方法: : y y n+1=y=y n+h(c+h(c1 1K K1 1+c+c2 2K K2 2) ) K K1 1=f(=f(xn,y,y n)

16、) （* * *）K K2 2=f(=f(xn+a+a2 2h,yh,y n+b+b2121hKhK1 1) )其中其中c c1 1,c,c2 2,a,a2 2,b,b2121為待定參數(shù)。對照式為待定參數(shù)。對照式( (* *) )有：有： (x,y,h)=c1 1f(x,y)+c2 2f(x+a2 2h,y+b2121hf(x,y)（2）二級二階）二級二階R-K方法方法(xn,y(xn);h)=(c1 1+c2 2)y(xn)+c2 2(a2 2hfx+b2 21 1hfyf)+O(h2 2)因?yàn)橐驗(yàn)閥(xn+1)在在xn處的處的Taylor展開為展開為 y(xn+1)=y(xn)+hy(xn

17、)+(h2/2!)y(xn)+O(h3)由由顯式單步法在顯式單步法在xn+1的局部截?cái)嗾`差的局部截?cái)嗾`差定義有：定義有： Tn+1=y(y(xn+1 )-y( )-y(xn )-h )-h(xn ,y(xn ),h) =h(1-c1-c2)y(xn )+h2(1/2-a2c2)fx+(1/2-c2b21)fyf+O(h3)顯然，若要求顯然，若要求Tn+1=O(h3),則應(yīng)有則應(yīng)有 c1+c2=1 c2a2=1/2 c2b21=1/2 當(dāng)當(dāng) =1時(shí)時(shí),c1=0,c2=1，得，得 yn+1= yn+hK2 n=0,1,.N-1 K1=f(xn,yn) K2=f(xn+h/2,yn+hK1/2)這就

18、是變形的這就是變形的歐拉方法或中點(diǎn)方法歐拉方法或中點(diǎn)方法。 二級二級R-KR-K方法是顯式單步式，每前進(jìn)一步方法是顯式單步式，每前進(jìn)一步需要需要計(jì)算兩個(gè)函數(shù)值計(jì)算兩個(gè)函數(shù)值。由上面的討論可知，由上面的討論可知，適當(dāng)選擇四個(gè)參數(shù)適當(dāng)選擇四個(gè)參數(shù)c c1 1,c,c2 2,a,a2 2,b,b2121, ,可使每步計(jì)可使每步計(jì)算兩次函數(shù)值的二階算兩次函數(shù)值的二階R-KR-K方法達(dá)到二階精度。方法達(dá)到二階精度。能否在計(jì)算函數(shù)值次數(shù)不變的情況下能否在計(jì)算函數(shù)值次數(shù)不變的情況下, ,通過通過選擇四個(gè)參數(shù)選擇四個(gè)參數(shù), ,使得二階使得二階R-KR-K方法的精度再提方法的精度再提高呢高呢? ? 答案是否定的

19、答案是否定的。無論四個(gè)參數(shù)怎樣選擇，。無論四個(gè)參數(shù)怎樣選擇，都不能使公式都不能使公式( (* * *) )提高到三階。提高到三階。 這說明每一步計(jì)算兩個(gè)函數(shù)值的二階這說明每一步計(jì)算兩個(gè)函數(shù)值的二階R-KR-K方法最高階為二階。若要獲得更高階得數(shù)值方法最高階為二階。若要獲得更高階得數(shù)值方法，就必須增加計(jì)算函數(shù)值的次數(shù)。方法，就必須增加計(jì)算函數(shù)值的次數(shù)。 仿照二級仿照二級R-KR-K方法，在方法，在 xn , , xn+1 上上, ,取取f f在在m m個(gè)個(gè)點(diǎn)的函數(shù)值做線性組合，即得到點(diǎn)的函數(shù)值做線性組合，即得到m m級級R-KR-K方法：方法：（3） m級顯式級顯式Runge-Kutta 方法方

20、法11111,2,3,mnnrrrnnrrnrnrsssyyhc KKfxyKfxhayhb Krm112341234123226,11,2211,22,nnnnnnnnnnhyyKKKKKfxyKfxh yhKfxh yhKfxh yKKhK 前面已經(jīng)看到前面已經(jīng)看到, ,二級、四級二級、四級R-KR-K方法可分別達(dá)到方法可分別達(dá)到最高階數(shù)二階、四階，但是最高階數(shù)二階、四階，但是N N級級R-KR-K方法的最高階卻方法的最高階卻不一定是不一定是N N階階。N N表示表示R-KR-K方法的級數(shù)表示公式中計(jì)方法的級數(shù)表示公式中計(jì)算函數(shù)值算函數(shù)值f f的次數(shù)。的次數(shù)。ButcherButcher給

21、出了給出了R-KR-K方法計(jì)算函方法計(jì)算函數(shù)值數(shù)值f f的次數(shù)與階數(shù)之間的關(guān)系表，如下：的次數(shù)與階數(shù)之間的關(guān)系表，如下：計(jì)算計(jì)算f f的次數(shù)的次數(shù) 1 2 3 1 2 3 4 4 5 6 7 5 6 7方法的最高階數(shù)方法的最高階數(shù) 1 21 2 3 3 4 44 54 5 6 6由表可見，四級以下由表可見，四級以下R-KR-K的方法其最高階數(shù)與計(jì)算的方法其最高階數(shù)與計(jì)算f f的次數(shù)一致，對的次數(shù)一致，對m m階階R-KR-K公式，當(dāng)公式，當(dāng)m4m4，雖然計(jì)算，雖然計(jì)算f f的的次數(shù)增加，但是方法階數(shù)不一定增加。因此次數(shù)增加，但是方法階數(shù)不一定增加。因此四級四四級四階階R-KR-K公式是應(yīng)用最為廣泛的公式。公式是應(yīng)用最為廣泛的公式。4 4、絕對穩(wěn)定性問題、絕對穩(wěn)定性問題(1)Euler方法的絕對穩(wěn)定性 將將Euler方法應(yīng)用于實(shí)驗(yàn)方程得到方法應(yīng)用于實(shí)驗(yàn)方程得到：1(,)(1)nnnnnnnyyhf x yyhyh y 這是一個(gè)齊次線性差分方程，這是一個(gè)齊次線性差分方程，其對應(yīng)的特征方程為：其對應(yīng)的特征方程為： 10rh 由定理由定理2 2可知，當(dāng)特征根可知，當(dāng)特征根 1rh 滿足滿足 11rh 時(shí)時(shí), , Euler方法是絕對