5數(shù)值積分與數(shù)值微分課件

1、PAGE PAGE 163第5章 數(shù)值積分與數(shù)值微分方法關于定積分計算，已經(jīng)有較多方法，如公式法、分步積分法等，但實際問題中，經(jīng)常出現(xiàn)不能用通常這些積分方法計算的定積分問題。怎樣把這些通常方法失效的定積分在一定精度下快速計算出來，特別是通過計算機編程計算出來就是本章研究的內(nèi)容。此外，怎樣根據(jù)函數(shù)在若干個點處的函數(shù)值去求該函數(shù)的導數(shù)近似值也是本章介紹的內(nèi)容。本章涉及的方法有Newton-Cotes求積公式、Gauss求積公式、復化求積公式、Romberg求積公式和數(shù)值微分。5.1 引 例人造地球衛(wèi)星軌道可視為平面上的橢圓。我國的第一顆人造地球衛(wèi)星近地點距離地球表面439km，遠地點距地球表面23

2、84km，地球半徑為6371km，求該衛(wèi)星的軌道長度。本問題可用橢圓參數(shù)方程來描述人造地球衛(wèi)星的軌道，式中a, b分別為橢圓的長短軸，該軌道的長度L就是如下參數(shù)方程弧長積分但這個積分是橢圓積分，不能用解析方法計算。5.2問題的描述與基本概念要想用計算機來計算，應對其做離散化處理。注意到定積分是如下和式的極限要離散化，做去掉極限號將取為具體的值為減少離散化帶來的誤差，將用待定系數(shù)代替于是就得到定義5.1 若存在實數(shù)且任取都有 （5.1）則稱式(5.1)為一個數(shù)值求積公式。式中稱為求積系數(shù)，稱為求積節(jié)點，而稱 （5.2）為求積余項或求積公式（5.1）的截斷誤差。從定義可以看到，數(shù)值求積公式依賴于求

3、積節(jié)點個數(shù)n、求積節(jié)點和求積系數(shù)，這三個量有一個發(fā)生變化，則產(chǎn)生不同的求積公式。定義5.2 若求積公式對所有不超過m次的多項式有求積余項，而對某一個m+1次多項式有，則稱該求積公式的代數(shù)精度為m。一般，一個求積公式的代數(shù)精度越大，則該求積公式越好。確定代數(shù)精度的方法依次取代入公式并驗證是否成立。若第一個使不成立的k值為m，則對應的代數(shù)精度為m-1。例 5.1確定求積公式的代數(shù)精度。解 取代入求積公式有易驗證，但，故本題求積公式代數(shù)精度為3。例 5.2確定下面求積公式的參數(shù)A，B，C，使它具有盡可能高的代數(shù)精度，并指出相應的代數(shù)精度。解 本題要先求出具體的求積公式，然后再判斷所求公式的代數(shù)精度。

4、公式有3個待定參數(shù)，h不是求積公式的參數(shù)，故利用3個條件得到的3個等式關系就可以解決求出具體求積公式的問題。依次取代入求積公式并取等號，有A+B+C=4hA-2C=0A+4C=163h解之得A=169h , B=43h , C=89 h故所求的求積公式為為確定其代數(shù)精度，再取代入求出的公式繼續(xù)計算，有，故所求的求積公式具有二階代數(shù)精度。5.3 插值型求積公式借助多項式插值函數(shù)來構造的求積公式稱為插值型求積公式。一般選用不同的插值公式就可以得到不同的插值型求積公式?；舅枷肜帽环e函數(shù) QUOTE fx 的插值函數(shù)x代替 QUOTE fx 做定積分的近似計算來構造求積公式。1.構造原理考慮 QU

5、OTE fx 在n個節(jié)點 QUOTE x1,x2,xn 上的n-1次Lagrange插值多項式 QUOTE Ln-1(x) 與 QUOTE fx 的余項，有這里 QUOTE lin-1x=k=1kinx-xkxi-xk, nx=k=1n(x-xk) 。兩邊取積分,有記 （5.3）則有 （5.4）若舍去，得求積公式求積系數(shù)的求積公式就是插值型求積公式。插值型求積公式的求積余項當為次數(shù)小于n次的多項式時，有 QUOTE fn(x)0 ，對應的。因此插值型求積公式的代數(shù)精度至少為n-1。若取，代入式(5.4)，可得插值型求積公式的求積系數(shù)之和為下面具體介紹常用的幾個插值型求積公式。2. Newton

6、-Cotes求積公式1) n點的Newton-Cotes公式的構造將求積節(jié)點 QUOTE xi 取為a,b上的等距節(jié)點做積分變量變換: 則當 QUOTE xa,b 時,有 QUOTE t0,n-1 ,于是有插值型求積公式的求積系數(shù)為記 QUOTE ci(n-1)=1n-10n-1k=1k1nt-k+1i-kdt ，則有 QUOTE cin-1 常稱為Cotes系數(shù)，易驗證通常稱 （5.6）為n點的Newton-Cotes公式。由于求積節(jié)點 QUOTE xi 是等距的，因此也稱式（5.6）為等距節(jié)點求積公式。利用可以得出下面常用的Newton-Cotes公式A) 2 點的Newton-Cotes

7、公式 （5.7）這正是我們熟悉的梯形公式。B) 3點的Newton-Cotes公式為 （5.8）稱它為Simpson公式或拋物線公式。表5.1 部分Cotes系數(shù)n2345678912 1216 46 1618 38 38 18790 1645 215 1645 79019288 2596 25144 25144 2596 1928841840 935 9280 34105 9280 935 41840751172803577172801323172802989172802989172801323172803577172807511728098928350588828350-928283501

8、049628350-4540283501049628350-9282835058882835098928350例5.3 試分別用梯形公式和Simpson公式計算解 用梯形公式計算，有I1-02f0+f1=121+sin10.920 735 49用Simpson公式計算，有I1-06f0+4f12+f1=161+8sin12+sin10.946 145 882）n點Newton-Cotes公式的代數(shù)精度 定理5.1 當求積節(jié)點個數(shù)n為奇數(shù)時，對應的Newton-Cotes求積公式的代數(shù)精度至少為n。證明 由于是插值型求積公式，故有對 QUOTE xn 有記 QUOTE s=k=12l+1s+l-

9、k+1=m=-lls+m, m=l-k+1 ，易知，故s是奇函數(shù)，得 QUOTE Rxn=0 ，得證。3） 梯形公式與Simpson公式的余項引理 5.1 (積分中值定理)設，且在上不變號，則有梯形公式余項為在a,b不變號，有2階連續(xù)導數(shù)，由引理 5.1，有梯形公式余項 （5.9）拋物線公式的余項 （5.10）4）Newton-Cotes公式的數(shù)值穩(wěn)定性設計算函數(shù) QUOTE fxk 時產(chǎn)生舍入誤差為 QUOTE k 實際在計算機中參加計算的是 QUOTE fxk 的近似值故Newton-Cotes公式在計算機中產(chǎn)生的誤差為若記 QUOTE =max1kn|k| ，則有由Cotes表5.1,當

10、n8時, ， QUOTE ck(n-1)0 從而有說明此時計算舍入誤差可以控制，從而Newton-Cotes公式是數(shù)值穩(wěn)定的。但當n8時， QUOTE ckn-1 有正有負， QUOTE k=1n|ckn-1| 隨n增大而增大，從而導致舍入誤差增加。故n8時，Newton-Cotes公式是數(shù)值不穩(wěn)定的。因而一般不用n8的Newton-Cotes公式來做定積分計算。3. Gauss求積公式1）Gauss求積公式的構造與概念n點的插值型求積公式的代數(shù)精度至少是n-1 ,那么是否還能提高其代數(shù)精度呢？若能，其代數(shù)精度最大能是多少？為回答這個問題，觀察一下插值求積公式的構造方法，發(fā)現(xiàn)其至少具有n-1次

11、代數(shù)精度的結論是在限定求積節(jié)點 QUOTE xi 的前提下得出的，若讓求積節(jié)點 QUOTE xi 也可以自由取值，則就給提高代數(shù)精度創(chuàng)造了條件。為使問題討論更具一般性，這里考慮帶權的定積分求積公式 （5.11）式中 QUOTE x 是已知的非負函數(shù)，為區(qū)間a,b上的權函數(shù)，a,b QUOTE a,b 可以取為 QUOTE 。顯然在式（5.11）中，若 QUOTE x1 ，就是前面討論的求積公式。定理5.2 求積公式(5.11)的代數(shù)精度最大為 QUOTE 2n-1 。 證明 設 QUOTE x1,x2,xn 是求積公式(5.11)的任意一組求積節(jié)點，用此節(jié)點構造插值型求積公式(5.11)，并令

12、取2n次多項式代入公式(5.11) QUOTE fx=wn2x , 考慮它的求積余項，有因為 QUOTE wn2x 是 QUOTE 2n 次多項式，由代數(shù)精度定義得的代數(shù)精度不大于 QUOTE 2n-1 。為證求積公式的代數(shù)精度可以為 QUOTE 2n-1 ，設 QUOTE f(x) 是任意一個 QUOTE 2n-1 次多項式，用 QUOTE wnx 去除 QUOTE f(x) ，由多項式除法有 QUOTE fx=qxwnx+rx,fxi=rxi (i=1,2,n) 式中 QUOTE qx，rx 都是次數(shù)不大于 QUOTE n-1 的多項式，于是有(5.12)因為上面的求積公式(5.11)是具

13、有n個節(jié)點的插值型求積公式，故其代數(shù)精度不小于 QUOTE n-1 ，從而有要讓式（5.11）成為等式，由式（5.12）應有 （5.13）式（5.13）要求對固定的n次多項式和任意次數(shù)不超過 QUOTE n-1 的多項式 QUOTE qx 都成立，這樣可以用 QUOTE qx 的這種任意性，選擇求積節(jié)點 QUOTE x1,x2,xn。 。由正交多項式理論可知：使式（5.13）成立的點 QUOTE x1,x2,xn 是存在唯一的，且都在 QUOTE a,b 內(nèi)，它就是在 QUOTE a,b 關于權 QUOTE x 的 QUOTE n 次正交多項式的零點。于是選取這些點作為求積公式的求積節(jié)點，并構

14、造對應的插值型求積公式，就得到具有 QUOTE 2n-1 次迭代精度的求積公式了。關于插值型求積公式的結論定理5.3 QUOTE n 點插值型求積公式的代數(shù)精度至少是 QUOTE n-1 ，至多為 QUOTE 2n-1 。定義5.3 若點的求積公式具有 QUOTE 2n-1 次代數(shù)精度，則稱該求積公式為Gauss型求積公式，對應的求積節(jié)點 QUOTE xi 和求積系數(shù)分別稱為Gauss點和Gauss系數(shù)。例5.4確定參數(shù) QUOTE x1,x2,A1,A2 ，使求積公式成為Gauss求積公式。解 注意到被積函數(shù)中有一因式與求積公式右端無明顯的關系，可將其視為權函數(shù)。為確定四個參數(shù)，依次取代入公

15、式并將近似號取為等號，得聯(lián)立方程組解出由于本題求積節(jié)點個數(shù)為 QUOTE n=2 ，其最高代數(shù)精度 QUOTE 2n-1=3 ，故所求 QUOTE x1,x2,A1,A2 的求積公式代數(shù)精度至少是3，故它是Gauss公式。2) Gauss型求積公式的求積余項設 QUOTE f(x)C2na,b ,取 QUOTE f(x) 在Gauss點 QUOTE x1,x2,xn 上的 QUOTE 2n-1 次Hermite插值多項式 QUOTE H2n-1(x) ,由Hermite插值余項公式，有兩邊積分，有由Gauss型求積公式的代數(shù)精度為 QUOTE 2n-1 及積分中值定理有，Gauss求積余項為3

16、) Gauss型求積公式的數(shù)值穩(wěn)定性對任意i, 取 QUOTE fx=lin-12x, lin-1x 是關于Gauss點 QUOTE x1,x2,xn 的 QUOTE n-1 次Lagrange插值函數(shù)，有 QUOTE lin-1xj=ij ,由Gauss公式的代數(shù)精度為 QUOTE 2n-1 ，而 QUOTE fx 是 QUOTE 2n-2 次多項式，有由i的任意性可知，Gauss求積系數(shù)。 QUOTE Ai0 若在Gauss公式中取 QUOTE fx1 ，可得類似Newton-Cotes穩(wěn)定性處理方法，有在計算機計算時的舍入誤差這說明Gauss型求積公式是穩(wěn)定的。在Gauss公式中，不同的

17、權函數(shù) QUOTE x 和不同積分區(qū)間，對應不同形式的Gauss公式，但最基本和常用的Gauss公式有4種。4） 常用的Gauss型求積公式Gauss-Legendre求積公式權函數(shù) QUOTE x1 ，積分區(qū)間為 QUOTE -1,1 , Gauss點為 QUOTE n 次Legendre正交多項式的零點，Gauss-Legendre求積公式為Gauss-Legendre求積余項為Gauss-Legendre求積公式的Gauss點 QUOTE xi 與系數(shù) QUOTE Ai 表 5.2nx kA k nx kA k 10.000 000 02.000 000 040.861 136 30.3

18、47 854 820.577 350 31.000 000 00.339 981 00.652 145 230.000 000 00.888 888 950.000 000 00.568 888 90.774 596 70.555 555 60.906 179 80.236 926 90.538 469 30.478 628 7Gauss-Legendre求積公式可以計算任何有限積分區(qū)間的定積分，計算之前先作變量代換將積分區(qū)間 QUOTE a,b 變到 QUOTE -1，1 ，有然后再對用Gauss-Legendre求積公式。Gauss-Chebyshev求積公式權函數(shù) QUOTE x=11-

19、x2 ,求積區(qū)間為 QUOTE -1，1 ，Gauss點為 QUOTE n 次Chebyshev正交多項式的零點，Gauss-Chebyshev求積公式為式中Gauss點 QUOTE xi 與系數(shù)Gauss-Chebyshev求積余項Gauss-Laguerre求積公式權函數(shù) QUOTE x ,積分區(qū)間為 QUOTE 0，+ ，Gauss點為 QUOTE n 次Laguerre正交多項式的零點，Gauss-Laguerre求積公式為Gauss-Laguerre求積余項為Gauss-Laguerre求積公式的Gauss點和系數(shù)也有事先計算好的表。（4）Gauss-Hermite求積公式權函數(shù) Q

20、UOTE x=e-x2 ,積分區(qū)間(-,+), Gauss點為n次Hermite正交多項式的零點，Gauss-Hermite求積公式為Gauss-Hermite求積余項為Gauss-Hermite求積公式的Gauss點及系數(shù) QUOTE Ak 有表。例5.5用兩點Gauss公式求定積分的近似值。解 本題為有限區(qū)間的定積分，可用兩點Gauss-Legendre 求積公式計算。做積分換元，將其化為-1,1上的定積分，即令有本題準確值為I=3.140 52208，可見精度很高。5.4復化求積公式Newton-Cotes公式在n8時數(shù)值不穩(wěn)定，因此不能用增加求積節(jié)點的方法來提高計算精度。實用中常用復化

21、求積公式來求積區(qū)間a,b上的定積分，以獲得滿足給定計算精度要的定積分值。常用的復化求積公式有復化梯形公式和復化Simpson公式。 基本思想將求積區(qū)間a,b分成若干個小區(qū)間，然后在每個小區(qū)間上采用數(shù)值穩(wěn)定的Newton-Cotes公式求小區(qū)間上的定積分，最后把所有小區(qū)間上的計算結果相加起來作為原定積分的近似值。復化梯形公式1）復化梯形公式的構造原理取等距節(jié)點 QUOTE xk=a+kh, h=b-an(k=0,1,n) 將積分區(qū)間a,b n等分,在每個小區(qū)間 QUOTE xk,xk+1(k=0,1,n-1) 上用梯形公式做近似計算，就有得求積公式復化梯形公式： (5.14)2）復化梯形公式的余

22、項記(5.15)故復化梯形公式的求積余項 由此可知,復化梯形公式的代數(shù)精度是1。若 QUOTE |fx|M2 ，對給定計算精度 QUOTE ，令得出說明利用復化求積公式能得到計算誤差小于的定積分近似值。2. 復化Simpson公式復化Simpson公式的構造原理取等距節(jié)點 QUOTE xk=a+kh, h=b-an(k=0,1,n) 將積分區(qū)間a,b n等分,在每個小區(qū)間 QUOTE xk,xk+1(k=0,1,n-1) 上用Simpson公式做近似計算，再累加起來就有式中 QUOTE xk+12=xk+h2 ，得復化Simpson公式 (5.16)2）復化Simpson公式的余項記有復化Si

23、mpson公式的求積余項從復化simpson公式的求積余項可以看出復化simpson公式的代數(shù)精度是3，它在代數(shù)精度和計算精度上都比復化梯形公式好。復化Simpson公式也稱為復化拋物線公式。例5.6分別用復化梯形公式和復化Simpson公式計算，要求誤差不超過 QUOTE 0.510-3 。解 為較快得結果，積分區(qū)間分割數(shù)按 QUOTE 2n(n=0,1,) 進行。數(shù)值計算結果列表，其中 QUOTE Rn 代表求積余項。N復化梯形公式復化Simpson公式TnRnSnRn2-17.389 2595.32-11.592 840-0.47822-13.336 0231.27-11.984 944

24、-0.85410-123-12.382 1620.312-12.064 209-0.61410-224-12.148 0040.77710-1-12.069 951-0.39510-325-12.089 7420.19410-126-12.075 1940.48510-227-12.071 5580.12110-228-12.070 6490.30310-3本題積分的準確值為I=-12.070 346 316 4，可見復化梯形公式和復化Simpson公式能求出精度較高的解。例5.7考慮用復化Simpson公式計算要使誤差小于0.510-6，那么求積區(qū)間0,1應分成多少個子區(qū)間。以此計算積分近似

25、值。解 復化Simpson公式的求積余項為 QUOTE Rf, Sn =-1-02590h4f4() QUOTE (0,1) 式中 QUOTE h=1-0n=1n, fx=sinxx 。為估計誤差，要計算 QUOTE f4(x) 。注意到，故由此得從而有解出 QUOTE n10472 ，故要求出滿足計算精度要求的定積分值，只要將0,1分成4個子區(qū)間即可，此時可算出5.5 Romberg 求積方法Romberg 求積方法是對復化梯形公式用加速技術得到的一種求積方法，它也稱為逐次分半加速收斂法。基本思想將Richardson 外推算法應用于復化梯形公式中，用產(chǎn)生的加速數(shù)列來求定積分值。1.構造原理

26、定理 5.4 設函數(shù) QUOTE F1(h) 逼近數(shù) QUOTE F* 的余項為 QUOTE F*-F1h=1hP1+2hP2+khPk+ (0P1P2) 式中 QUOTE Pk,k 都是與h無關的常數(shù)，且 QUOTE k1 時，,則由 ，q為任意常數(shù)定義的函數(shù) QUOTE F2h 也逼近 QUOTE F* ，且有式中 QUOTE k1(k2) 都是與h無關的非零常數(shù)。證明 用qh 代替余項式（5.18）的變元h，有(5.19)用 QUOTE -qP1 乘式（5.18）并與式(5.19)相加，整理后，有(5.20)因為0q 1，故 QUOTE 1-qP10 ，用 QUOTE 1-qP1 同除式

27、(5.20)，有令 QUOTE F2h=F1qh-qP1F1h1-qP1 QUOTE k1=kqPk-qP11-qP1 即得所證。由微積分收斂階概念可知 QUOTE F2h 的收斂階 QUOTE P2 比的收斂階 QUOTE P1 大，故 QUOTE F2h 比 QUOTE F1h 逼近的速度快。通常稱如上的方法為Richardson 外推法。顯然這種外推可以不斷做下去以獲得逼近更快的函數(shù)，一般有 QUOTE Fm+1h=Fmqh-qPmFmh1-qPm (m=1,2,3) QUOTE F*-Fm+1h=m+1(m)hPm+1+m+1(m)hPm+2+ 這樣，每用一次Richardson外推，

28、就使逼近階提高一個等級，從而達到加速的目的。分析記 稱 QUOTE I=abf(x)dx 為復化梯形公式的T形值，可以證明：此式說明 QUOTE T2h 逼近 QUOTE I 的階為 QUOTE O(h4) 。利用Richardson 外推法對做加速。因為,有 QUOTE T2h=T1qh-q2T1h1-q2 QUOTE T2h 逼近的階變?yōu)?QUOTE O(h4) 。若取 QUOTE q=12 ，則有同理對 QUOTE T2h 再做一次Richardson 加速，有逼近I的階為。一般，有經(jīng)Richardson加速的求定積分的序列為 (5.21) QUOTE Tmh=4m-1Tm-1h2-Tm

29、-1(h)4m-1-1 QUOTE T2h 逼近 QUOTE I I的階變?yōu)?QUOTE O(h4) 。 注意到 QUOTE T1h2=T2n ，直接計算可知 QUOTE T2h 就是復化Simpson公式的 QUOTE Sn ，即有 QUOTE Sn=4T2n-Tn4-1 (5.22)類似地可得 QUOTE Cn=42S2n-Sn42-1 (5.23) QUOTE Rn=43C2n-Cn43-1 (5.24)這樣可以用復化梯形公式計算出序列T2k。再逐次用公式（5.22），（5.23），（5.24）對其進行加工得到收斂更快的序列S2k，C2k，R2k等等。通常將序列T2k，S2k，C2k和R

30、2k依次稱為梯形序列、Simpson序列、Cotes序列和Romberg序列。若繼續(xù)外推，舍入誤差增大將使新序列體現(xiàn)不出更明顯的加速效果，因此一般只外推到Romberg序列。用Romberg序列求定積分的算法稱為Romberg求積方法。3Romberg求積方法的計算過程kT2kS2kC2kR2k0 eq oac(,1)T1 eq oac(,3)S1 eq oac(,6)C1 eq oac(,10)R11 eq oac(,2)T2 eq oac(,5)S2 eq oac(,9)C2 eq oac(,14)R22 eq oac(,4)T4 eq oac(,8)S4 eq oac(,13)C43 e

31、q oac(,7)T8 eq oac(,12)S84 eq oac(,11)T16因為 QUOTE Tn 與 QUOTE T2n 有相重的點，計算時還可以利用下面關系式來減小計算量。在給定計算誤差界 QUOTE 后，可用R2k+1-R2k來作為終止Romberg計算的條件。例5.8 用Romberg算法計算計算結果精確到0.510-6。解 由Romberg求積方法，得計算結果如下表kT2kS2kC2kR2k00.920 73550.946 14590.946 08300.946 083110.939 79330.946 08690.946 08310.946 083020.944 51350.

32、946 08340.946 083030.945 69090.946 083040.9845 9850由 ，故有。5.6 數(shù)值微分根據(jù)函數(shù)在若干個點處的函數(shù)值去求該函數(shù)的導數(shù)近似值稱為數(shù)值微分，所求導數(shù)的近似值常稱為數(shù)值導數(shù)。設是已知的個點處的函數(shù)值，下面分別討論兩種求數(shù)值微分的方法。1.利用次多項式插值函數(shù)求數(shù)值導數(shù)用插值函數(shù)代替被插值函數(shù)可以用來計算被插函數(shù)的近似值；用插值函數(shù)代替定積分的被積函數(shù)可以用來計算定積分的近似值；用插值函數(shù)的導數(shù)能否用來求被插函數(shù)導數(shù)值？定義5.4 設是的次插值多項式，稱用插值函數(shù)的導數(shù)來計算被插函數(shù)導數(shù)的公式 為的插值型求導公式。由插值理論，有特別有一階數(shù)值導數(shù)的余項關系在非節(jié)點處的余項由于涉及對未知中值函數(shù)的求導，使余項不易描述和控制，但在節(jié)點處有此余項