Ch數(shù)值計算方法之數(shù)值積分實用

會計學1Ch數(shù)值計算方法之數(shù)值積分實用1．術語和記號為了計算f(x)在區(qū)間[a,b]上的定積分近似值，我們通常的做法是，把積分區(qū)間[a,b]劃分為n等分，記h=(b-a)/n,x0=a,xk=a+kh,k=0,1,2,…,n,稱x0,x1,…,xn為[a,b]的一個等份分劃。假如x0,x1,…,xn為[a,b]的一個等份分劃那么求積公式(1)中的w0,w1,…,wn的選取僅僅只與n有關，從而可以簡化對求積公式的研究。結論，只要給出了一個如何確定（1）式中的諸w0,w1,…,wn的機制，我們就可以得到相應的對任何被積函數(shù)都有效的計算定積分方法。第1頁/共36頁2.求積公式的性質微積分學中我們曾研究過，定積分保持函數(shù)的線性關系不變，它的含義是，若f(x),g(x)都是[a,b]上的可積函數(shù)，則對任意實數(shù)u,v,我們有u·f(x)+v·g(x)也是[a,b]上的可積函數(shù)，而且不難驗證，求積公式也保持函數(shù)的線性關系不變，即第2頁/共36頁3.幾種常見的求積公式在后面的討論中，我們將經常用到下面一些非常簡單的求積公式，他們是中點公式、梯形公式和辛卜生公式。中點公式我是我們課程中強調的一個名詞，與求積公式（1）對比分析，可以認為它是這樣一種機制：把積分區(qū)間分為2等分，取w0=w2=0，w1=1所形成的求積公式。從幾何上看，它實際上是取區(qū)間中點的函數(shù)值與區(qū)間長的積作為定積分值，類似于用中位線乘以高來計算梯形的面積。第3頁/共36頁4．截斷誤差在求積公式中，我們使用的是近似等號，這是因為，對于一般的被積函數(shù)來說，利用這些公式計算所得的結果除了舍入誤差外，還有截斷誤差，因為定積分是用極限來定義的。有時為了進行誤差分析，我們可以把上面的（1）式寫成（1’）其中R[f]表示的就是截斷誤差?？疾烨懊娼o出的三個求積公式，如果被積函數(shù)是線性函數(shù)，那么利用中點公式或梯形公式所得到的結果就是準確值，否則一般不是。對于一般的非線性函數(shù)，感覺上辛卜生公式更好一些。為了刻劃求積公式對一般的被積函數(shù)的精確度，我們引進代數(shù)精度的概念。第4頁/共36頁5.代數(shù)精度的概念定義：一個求積公式如果對所有的次數(shù)不超過m的多項式嚴格相等，而對某些m+1次多項式不相等，則稱該公式具有代數(shù)精度m，或該公式的代數(shù)精度為m。利用求積公式的線性性，我們不難證明下面的結論。定理：如果求積公式對1，x,…,xm嚴格相等，而對xm+1不相等，則該公式的代數(shù)精度為m。作為課外練習，鼓勵大家給出完整證明。第5頁/共36頁6.基本結論我們可以利用上面的定理所給出的方法證明辛卜生公式的代數(shù)精度是3，而中點公式和梯形公式的代數(shù)精度是1?，F(xiàn)在我們可以對這三個公式作一個簡單的評價：中點公式和梯形公式的代數(shù)精度雖然都是1，但中點公式只計算一個點的函數(shù)值，而梯形公式卻要計算兩個點處的函數(shù)值，所以中點公式優(yōu)于梯形公式。與梯形公式相比，辛卜生公式只多計算一個點的函數(shù)值，但代數(shù)精度卻增加到3，顯然辛卜生公式更為優(yōu)越。第6頁/共36頁10．2牛頓-柯特斯求積公式牛頓-柯特斯求積公式就是利用Lagrange插值多項式導出的求積公式。把一般的函數(shù)的積分轉化為相應的插值多項式函數(shù)的積分也是我們學習插值法的基本目的之一。第7頁/共36頁1.利用插值多項式近似替代被積函數(shù)設f(x)為被積函數(shù)，[a,b]為積分區(qū)間，x0,x1,…,xn為[a,b]內的n+1個互異的點，記Ln(x)為相應的拉格朗日插值多項式，那么我們有f(x)=Ln(x)+Rn(x)第8頁/共36頁2.利用插值多項式導出求積公式提示：在上面給出的公式中，由于諸lk(x)都是多項式函數(shù)，所以諸wk都可以精確地計算出來。從而可以得到一般性的求積公式。第9頁/共36頁2.利用插值多項式導出求積公式（注釋）回顧上一章關于多項式插值的結論，由于任意次數(shù)不超過n的多項式與它的任意n+1個基點的插值多項式恒等，再由求積公式的代數(shù)精度的定義，我們立即得到：由n+1個基點的拉格朗日插值多項式所形成的求積公式的代數(shù)精度至少式n，為此，我們上面的wk改寫為w(n,k),k=0,1,…,n。第10頁/共36頁3.牛頓-柯特斯求積公式牛頓-柯特斯求積公式就是利用等距基點的拉格朗日插值多項式導出的求積公式。將積分區(qū)間[a,b]劃分為n等分，記h=(b-a)/n，取x0=a,xk=a+kh.k=0,1,…,n我們可以得到第11頁/共36頁3.牛頓-柯特斯求積公式（注釋）在牛頓-柯特斯公式中，我們稱c(n,k)為牛頓-柯特斯系數(shù)，一般可通過查表得到。崔國華教材p60列出了直到n=8的所有牛頓-柯特斯系數(shù)，應該說，實用意義不大。當n>8時牛頓-柯特斯公式并沒有實際意義。實際上，我們通常只用到n=1,2,4的情形，相應的公式分別稱為梯形公式，辛卜生公式，和柯特斯公式。對于現(xiàn)代的計算工具來說，有梯形公式和辛卜生公式也就夠用了。利用牛頓-柯特斯系數(shù)，我們可以方便地寫出牛頓-柯特斯求積公式：第12頁/共36頁4.梯形公式在牛頓-柯特斯求積公式中，如果我們取n=1,那么k可以取0和1。由此所形成的求積公式就是梯形公式。由于得到的結果是梯形面積公式，所以稱梯形公式。第13頁/共36頁5.辛卜生公式在牛頓-柯特斯求積公式中，如果我們取n=2,那么k可以取0，1，2，由此所形成的求積公式就是辛卜生公式。第14頁/共36頁6.柯特斯公式作為課外作業(yè)，大家可以取n=4,相應地k可以取0，1，2，3和4，仿照上面的方式，可以得到：從而可進一步寫出相應的求積公式，這就是柯特斯公式。在后面將要介紹的龍貝格求積算法中，我們將產生梯形序列，辛卜生序列，柯特斯序列和龍貝格序列，前三個序列都是基于牛頓-柯特斯公式產生的序列，而龍貝格序列則不是。第15頁/共36頁10.3變步長復化梯形公式法在上一節(jié)中我們介紹了牛頓-柯特斯公式以及它的特款，并且得到了牛頓-柯特斯公式的代數(shù)精度為n.，接下來就是如何利用這個公式。一般說來，在一個較大的積分區(qū)間上利用較高階的牛頓-柯特斯公式雖然可以得到較高的代數(shù)精度，但實際效果并不好，道理也不難理解。我們可以把積分區(qū)間劃分為若干等分，在每個子區(qū)間上利用較低階的求積公式，(當然元可以是牛頓-柯特斯公式，也可以不是，)并把這些積分值相加。按照這樣的思路所得到的求積公式統(tǒng)稱為復化型求積公式。第16頁/共36頁1.復化中點公式復化中點公式也許最不為人們所注意，以至在一般的教科書中還沒有這個名稱，我們在后面將會看到，對于求數(shù)值積分來說，它實際上是最有用的公式。把積分區(qū)間[a,b]劃分為n等分，記x0,x1,…,xn為等分點，記[xj-1,xj]為第j個子區(qū)間，zj為區(qū)間的中點，j=1,2,…,n，記h=(b-a)/n,記Mn為所有子區(qū)間上利用中點公式所求得的積分值的和，那么我們有（1）第17頁/共36頁1.復化中點公式(等價形式)現(xiàn)在我們把積分區(qū)間[a,b]劃分為2n等分,記y0,y1,…,y2n-1,y2n為等分點，此時我們可以把下標為奇數(shù)的點y2k-1,看成是區(qū)間[y2(k-1),y2k]的中點，k=1,2,…,n。(如圖)y0y1y2y3y4y2(n-1)y2n-1y2nx0x1x2xn-1xn這樣，我們也可以把復化中點公式寫成（2）提示：這兩種表示方法我們在后面都會用到，所以一定要記住。不管公式的形式如何，編寫求Mn的程序還是相同的。第18頁/共36頁1.復化中點公式(c語言代碼)在實際應用中，n通常取2的某個整數(shù)冪2k，我們可以把計算結果作為數(shù)組M的第k個元素，程序段命名為,C語言源程序文件可以按如下的方式編寫。#include<math.h>#defineF(X)4.0/(1.0+X*X)staticdoublea=0.0,b=1.0;doubleGetM(intk){doublex,y,step;intn,j;n=1;for(j=0;j<k;j++)n*=2;step=(b-a)/n;x=a+step/2;for(j=0;j<n;j++){y+=F(x);x+=step;}return(step*y);}第19頁/共36頁2.復化梯形公式把積分區(qū)間[a,b]劃分為n等分，記x0,x1,…,xn為等分點，記h=(b-a)/n,則有記Ij,j=1,2,…,n為f(x)在第j個子區(qū)間[xj-1,xj]應用梯形公式所求得的積分值，則有第20頁/共36頁2.復化梯形公式(續(xù))

這就是我們的復化梯形公式。提示：雖然我們也可以直接編寫求Tn的計算機程序，但是沒有這個必要。第21頁/共36頁3.變步長復化梯形公式假設對某個n,我們利用復化梯形公式，也就是上面的（3）式，得到了Tn，如果它不滿足我們的精度要求，那么我們可以把每個子區(qū)間再對分一次，這相當于把積分區(qū)間劃分為2n等分。記y0,y1,y2,…,y2(n-1),y2n-1,y2n為等分點，記t=(b-a)/(2n),則有再利用（3）式即得第22頁/共36頁3.變步長復化梯形公式(續(xù))再利用上面的（2）式和（3）式，并注意到（3）式中的諸xk即為這里的y2k,k=1,2,…,n-1，所以我們有

（4）這就是變步長復化梯形公式。我們可以利用它形成一個自動調整精度的算法。第23頁/共36頁4.變步長復化梯形公式方法為了方便地利用（4）式形成一個算法，我們總是取n等于2的某個指數(shù)冪2k,并把Tn記為T[k],由n=2k可得2n=2k+1,從而有

T[k+1]=(T[K]+M[K])/2算法說明如下：T[0]=(b-a)*(f(a)+f(b)/2.0;For(K=0;k<20;k++){M[k]=GetM(k);T[k+1]=(M[k]+T[k])/2.0;if(fabs(M[k+1]-M[k])<EPS)break;}提示：如果取k=20，那么GetM(k)中計算函數(shù)值的次數(shù)超過1000000次。所以我們一般不讓k值超過20第24頁/共36頁12.4變步長復化辛卜生公式法變步長復化梯形公式法的核心思想是：假設對于n=2K我們計算出來了T[K],如果T[k]不滿足精度要求，則進一步把區(qū)間劃分為2K+1等份，也就是在原來的每個子區(qū)間補上中點的函數(shù)值，并利用T[k+1]=(T[k]+M[k])/2來計算T[k+1]現(xiàn)在我們也可以這樣考慮，利用原來每個子區(qū)間的端點函數(shù)值f(y2(k-1))和f(x2k)以及后來補算的中點的函數(shù)值f(y2k-1)采用辛卜生公式計算第k個子區(qū)間的積分值，再把他們相加結果計為S[K+1]。無論是理論分析還是直觀考慮，這樣做的效果應該更好一些,問題是能否也利用MK,TK來簡化我們的計算。第25頁/共36頁1.復化辛卜生公式假設x0,x1,x2,…,x2(k-1),x2k-1,x2k,…,x2(n-1),x2n-1,x2n把積分區(qū)間[a,b]劃分為2n等分，我們可以認為其中的x0,x2,…,x2k,…,x2n把區(qū)間劃[a,b]劃分為n等份，并且x2k-1就是第k個子區(qū)間[x2(k-1),x2k]的中點。記Ik,k=1,2,…,n為f(x)在第k個子區(qū)間[x2(k-1),x2k]應用辛卜生公式所求得的積分值，則有這就是我們的復化辛卜生公式,雖然我們也可以直接編寫求S2n的計算機程序，但是沒有這個必要。第26頁/共36頁2.關于復化辛卜生公式的結論假設x0,x1,x2,…,x2(k-1),x2k-1,x2k,…,x2(n-10,x2n-1,x2n把積分區(qū)間[a,b]劃分為2n等分，那么我們可以得到下面3個不同的積分值：利用x2k,k=0,1,…,n這n+1個點處的函數(shù)值和復化梯形公式計算出Tn;利用xj,j=0,1,…,2n這2n+1個點處的函數(shù)值和復化梯形公式計算出T2n;利用復化辛卜生公式計算出S2n;重要結論：對于上面的約定，我們有

S2n=(4T2n-Tn)/3第27頁/共36頁3．詳細推導

第28頁/共36頁4.變步長復化辛卜生公式加速算法我們仍然取n為2的整數(shù)冪的形式，即n=2k,并記Sn為S[k],那么我們有T[k+1]=(M[k]+T[k])/2S[k+1]=(4T[K+1]-T[K])/3.0我們可以在變步長復化梯形公式法基礎上稍加改進，即可得到一個相應的加速算法也就是簡單地利用T[K+1]和T[K]的線性組合直接得到S[k+1]，額外的代碼和計算量可忽略不計。第29頁/共36頁5.算法說明變步長復化辛卜生公式加速算法是一個實用的方法，（而不是像前面介紹的方法只是過渡性的方法）所以要求大家都會編程。第30頁/共36頁12.5龍貝格求積算法前面我們介紹過變步長復化梯形公式法，利用它可以形成一個自動調整精度的算法。接下來又介紹了如何利用復化辛卜生公式實現(xiàn)變步長復化梯形公式法的加速，基本思想是利用T[k+1]和T[k]的線性組合來產生S[k+1],我們稱為進行一次外推，顯然，它比起直接計算S[k+1]優(yōu)越的多。類似地