1、目錄第一章 數值積分計算的重述11.1引言11.2問題重述2第二章 復化梯形公式22.1 復化梯形公式的算法描述32.2 復化梯形公式在C語言中的實現32.3 測試結果4第三章 復化simpson公式53.1 復化simpson公式的算法描述53.2 復化simpson公式在C語言中的實現63.3 測試結果7第四章 復化cotes公式84.1 復化cotes公式的算法描述84.2 復化cotes公式在C語言中的實現94.3 測試結果10第五章 Romberg積分法115.1 Romberg積分法的算法描述115.2 Romberg積分法在C中的實現125.3 測試結果13第六章 結果對比分析和

2、體會14參考文獻16附錄17數值積分（一）第一章 數值積分計算的重述1.1引言數值積分是積分計算的重要方法，是數值逼近的重要內容，是函數插值的最直接應用，也是工程技術計算中常常遇到的一個問題。在應用上，人們常要求算出具體數值，因此數值積分就成了數值分析的一個重要內容。在更為復雜的計算問題中，數值積分也常常是一個基本組成部分。在微積分理論中，我們知道了牛頓-萊布尼茨(Newton-Leibniz)公式其中是被積函數的某個原函數。但是隨著學習的深入，我們發現一個問題: 對很多實際問題，上述公式卻無能為力。這主要是因為:它們或是被積函數沒有解析形式的原函數,或是只知道被積函數在一些點上的值,而不知道

3、函數的形式,對此,牛頓萊布尼茨(Newton-Leibniz)公式就無能為力了。此外,即使被積函數存在原函數,但因找原函數很復雜,人們也不愿花費太多的時間在求原函數上,這些都促使人們尋找定積分近似計算方法的研究,特別是有了計算機后,人們希望這種定積分近似計算方法能在計算機上實現,并保證計算結果的精度,具有這種特性的定積分近似計算方法稱為數值積分。由定積分知識,定積分只與被積函數和積分區間有關,而在對被積函數做插值逼近時,多項式的次數越高,對被積函數的光滑程度要求也越高,且會出現Runge現象。如時,Newton-Cotes公式就是不穩定的。因而,人們把目標轉向積分區間,類似分段插值,把積分區間

4、分割成若干小區間,在每個小區間上使用次數較低的Newton-Cotes公式,然后把每個小區間上的結果加起來作為函數在整個區間上積分的近似,這就是復化的基本思想。本文主要研究的公式有: 復化梯形公式復化Simpson公式復化Cotes公式Romberg積分法。1.2 問題重述本文主要介紹微積分方程的復化解法。通過運用復化梯形公式、復化Simpose公式、復化cotes公式和Romberg積分法這四種積分法方法，解出微分方程的近似解。并進行誤差分析和結果比較。當積分區間a,b的長度較大,而節點個數n + 1固定時，直接使用Newton-Cotes公式的余項將會較大，而如果增加節點個數,即n + 1

5、增加時，公式的舍入誤差又很難得到控制，為提高公式的精度,又使算法簡單易行,往往使用復化方法。即將積分區間a,b分成若干個子區間，然后在每個小區間上使用低階Newton-Cotes公式，最后將每個小區間上的積分的近似值相加。將定積分的積分區間a b分割為n等份各節點為 ,k=0,1,n 在子區間(k=0,1,1.n-1)上使用Newton Cotes公式將子區間分割為l等份,節點為記為在子區間上作f(x)的l階Newton-Cotes求積公式由積分的區間可加性,可得復化求積公式第二章 復化梯形公式2.1 復化梯形公式的算法描述復化求積公式當L=1時可得復化梯形公式： = 復化梯形公式=2.2 復

6、化梯形公式在C語言中的實現復化梯形公式運用的程序如下： T0=(a-b)*(f_x(a)+f_x(b)/2;/n=1時的cotes公式即梯形公式for(i=1;i<=100;i+) /計算sum_num、xishu、s_point(start point）、d_point sum_num=pow(2,i-1); xishu=double(a-b)/sum_num; s_point=double(b)+double(a-b)/pow(2,i); d_point=double(a-b)/pow(2,i-1); for(j=1;j<=sum_num;j+) add_T=add_T+f_x

7、(s_point+(j-1)*d_point); add_T=add_T*xishu; T1=(T0+add_T)/2; err_T=(T1-T0)/3; /output printf("%d %d %10.8f %10.8f",i,pow(2,i),T1,err_T); printf("n"); if(err_T<=0) err_T=(-1)*err_T; if(err_T<=E) break; else T0=T1; T1=0; add_T=0; err_T=0; 在這個函數中我們將復化梯形公式和積分過程都用計算機語言表示出來。首先我們給

8、出復化梯形公式，進行迭代，直到精確度達到設定要求，算出最后結果。2.3 測試結果用復化梯形有效數字四位求得的結果如下：用復化梯形有效數字七位求得的結果如下：由以上結果可以看出取兩個不同的精度相對誤差比較小，但計算次數大大的增加，復化梯形公式計算次數多。第三章 復化simpson公式3.1 復化simpson公式的算法描述復化求積公式當L=2時可得復化梯形公式： = 復化simpson公式= 3.2 復化simpson公式在C語言中的實現復化梯形公式運用的程序如下：T0=(a-b)*(f_x(a)+4*f_x(a+b)/2)+f_x(b)/6;/n=2的cotes公式即simpson公式 for

9、(i=1;i<=100;i+) /計算sum_num、xishu、s_point(start point）、d_point /long powl (long double x, long double y) sum_num=pow(2,i-1); /the same as T xishu=double(a-b)/sum_num/6; s_point=double(b)+double(a-b)/pow(2,i); d_point=double(a-b)/pow(2,i-1); sd_point=double(a-b)/sum_num/4; for(j=1;j<=sum_num;j+)

10、 add_T=add_T+2*f_x(s_point+(j-1)*d_point)-4*f_x(s_point-sd_point+(j-1)*d_point)-4*f_x(s_point+sd_point+(j-1)*d_point); add_T=add_T*xishu; T1=(T0-add_T)/2; err_T=(T1-T0)/15; /output printf("%d %d %10.8f %10.8f",i,pow(2,i),T1,err_T); printf("n"); if(err_T<=0) err_T=(-1)*err_T; i

11、f(err_T<=E) break; else T0=T1; T1=0; add_T=0; err_T=0; 在這個函數中我們將復化simpose公式和積分過程都用計算機語言表示出來。首先我們給出復化simpose公式，進行迭代，直到精確度達到設定要求，算出最后結果。3.3 測試結果用復化simpose迭代取有效數字四位求得的結果如下：用復化simpose迭代取有效數字七位求得的結果如下：由以上結果可以看出兩次不同精度要求的計算可以看出不同精度計算計算次數相差較多。精度為四和七間計算次數相差了三次。第四章 復化cotes公式4.1 復化cotes公式的算法描述復化求積公式當L=4時可得復

12、化梯形公式： =復化cotes公式=4.2 復化cotes公式在C語言中的實現復化cotes公式運用的程序如下： T0=(a-b)*(7*f_x(1)+32*f_x(2)+12*f_x(3)+32*f_x(4)+7*f_x(5)/90;/四階（n=4）cotes公式。 for(i=1;i<=100;i+) sum_num=pow(2,i-1); /the same as T xishu=double(a-b)/sum_num/90; s_point=double(b)+double(a-b)/pow(2,i); d_point=double(a-b)/pow(2,i-1); sd_poi

13、nt=double(a-b)/sum_num/8; for(j=1;j<=sum_num;j+) add_T=add_T-2*f_x(s_point+(j-1)*d_point)-32*f_x(s_point-sd_point+(j-1)*d_point)+20*f_x(s_point-2*sd_point+(j-1)*d_point)-32*f_x(s_point-3*sd_point+(j-1)*d_point)-32*f_x(s_point+sd_point+(j-1)*d_point)+20*f_x(s_point+2*sd_point+(j-1)*d_point)-32*f_x

14、(s_point+3*sd_point+(j-1)*d_point); add_T=add_T*xishu; T1=(T0-add_T)/2; err_T=(T1-T0)/63; /output printf("%d %d %10.8f %10.8f",i,pow(2,i),T1,err_T); printf("n"); if(err_T<=0) err_T=(-1)*err_T; if(err_T<=E) break;else T0=T1; T1=0; add_T=0; err_T=0; 在這個函數中我們將復化cotes公式和積分過程都用計

15、算機語言表示出來。首先我們給出復化cotes公式，進行迭代，直到精確度達到設定要求，算出最后結果。4.3 測試結果用復化cotes有效數字四位求得的結果如下：用復化cotes有效數字七位求得的結果如下：由以上結果可以兩次不同精度計算的結果相差相對前面的方法要大，計算次數增加了三次。第五章 Romberg積分法5.1 Romberg積分法的算法描述Romberg方法也稱為逐次分半加速法。它是在梯形公式、辛卜生公式和柯特斯公式之間的關系的基礎上，構造出一種加速計算積分的方法。 作為一種外推算法, 它在不增加計算量的前提下提高了誤差的精度。在等距基點的情況下，用計算機計算積分值通常都采用把

16、區間逐次分半的方法進行。這樣，前一次分割得到的函數值在分半以后仍可被利用,且易于編程 。對區間a， b，令h=b-a構造梯形值序列T2K。                 T1=hf(a)+f(b)/2 把區間二等分，每個小區間長度為 h/2=(b-a)/2，于是  T2 =T1/2+h/22f(a+h/2）     把區間四(22)等分，每個小區間長度為h/22 =（b-a）

17、/4，于是  T4 =T2/2+h/2f（a+h/4）+f（a+3h/4).把a，b 2k等分，分點xi=a+（b-a）/ 2k ·i （i =0，1，2 · · · 2k）每個小區間長度為（b-a）/ 2k ，由歸納法可得面所說的的第一個公式.（二）計算公式如下：整個程序就是循著這四個公式進行計算的。Sn,Cn, Rn 分別代表特例梯形積分,拋物線積分,龍貝格積分.當然,編程的時候統一處理即可。5.2 Romberg積分法在C中的實現Romberg公式運用的程序如下：double T0=0,S0=0,C0=0,T1=0,S1=0,C

18、1=0,R0=0,R1=0; double err_T=10; int i=0,j=0; int sum_num=0; double xishu=0;/系數double s_point=0,d_point=0; double add_T=0; T0=(a-b)*(f_x(a)+f_x(b)/2;/n=1時的cotes公式即梯形公式。 for(i=1;i<=100;i+)/the first base number sum_num=pow(2,i-1); xishu=double(a-b)/sum_num; s_point=double(b)+double(a-b)/pow(2,i); d

19、_point=double(a-b)/pow(2,i-1); for(j=1;j<=sum_num;j+) add_T=add_T+f_x(s_point+(j-1)*d_point); add_T=add_T*xishu; T1=(T0+add_T)/2; add_T=0; /計算 S1 S1=4*T1/3-T0/3; if(i>=2) C1=16*S1/15-S0/15; if(i>=3) R1=64*C1/63-C0/63; /check using the "1" data err_T=(R1-R0)/255; if(err_T<0) err

20、_T=(-1)*err_T; if(err_T<E)&&(i>=4) break; /完成計算后，準備下一次循環 T0=T1; T1=0; S0=S1; S1=0; C0=C1; C1=0; R0=R1; R1=0; 在這個函數中我們將romboerg公式和的積分過程都用計算機語言表示出來。首先我們給出romboerg公式的T0，進行迭代，分別算出S1，C1，R1直到精確度達到設定要求，算出最后結果。5.3 測試結果用romboerg有效數字四位求得的結果如下：用romboerg有效數字七位求得的結果如下：由以上結果可看出，用romboerg取不同的精度對T1，S1

21、，C1，R1的結果影響大小不相同，T1影響最大，R1影響最小，迭代次數越多精度影響大小越小。第六章 結果對比分析和體會 通過對不同精度的測試發現復化梯形公式的計算量增加最快，而romberg達到一定的精度要求結果無法正常計算顯示。如下圖所示當精度要求達到20時結果無法正常顯示。而其他可正常顯示結果但是計算次數相對較大如復化梯形計算次數為三十三次，由以上程序測試的數據結果的對比顯示可知不同求積公式各有特點.梯形求積公式和Simpson求積公式雖然計算簡單、使用方便, 但是精度較差, 但對于光滑性較差的被積函數有時比高精度方法更為有效。尤其梯形公式對被積函數是周期函數的效果更為突出。 時,復化梯形

22、公式和復化Simpson公式在保留了低階公式的優點, 又能獲得較高的精度, 因此在實際計算中應用的最為廣泛。利用二分技術得到的Romberg方法的算法簡單, 易于編程實現。當節點加密提高積分近似程度時, 前面計算的結果可以為后面所用, 對減少計算量很有好處, 并有比較簡單的誤差估計, 能得到若干積分序列, 如果在做收斂性控制時, 同時檢查各行、各列, 對于不同性態的函數可以用其中最快的收斂序列來逼近積分。參考文獻1 李慶揚，王能超，易大義.數值分析M.武漢.華中科技大學出版社，2006.7.2 清華大學、北京大學計算方法編寫組.計算方法M .北京.科學出版社，1980 3 呂同斌復化梯形公式及

23、其應用期刊論文 安徽水利水電職業技術學院學報 2002 年4 期4 溪梅成數值分析方法 M 合肥：中國科學技術大學，2003附錄1.復化梯形源程序#include <stdio.h>#include <math.h>#define a 5#define b 1#define E 0.00000005 /即保留七位有效數字0.5*10-7/原函數double f_x(double x) double y; y=exp(-(x*x); return(y);int pow(int x,int y) int z=1; int i; for(i=0;i<y;i+) z=z*

24、x; return z;void main() /計算T1,T2,T4,T8. double T0=0,T1=0; double err_T=10; int i=0,j=0; int sum_num=0; double xishu=0; double s_point=0,d_point=0; double add_T=0; T0=(a-b)*(f_x(a)+f_x(b)/2;/n=1時的cotes公式即梯形公式 printf("n"); printf("=數值積分_利用復化梯形公式=n"); printf("n");printf(&q

25、uot;i 2i T_2i (T_2i-T_2(i-1)/3 n"); printf("n"); printf("0 %d %10.8f 0",pow(2,0),T0); printf("n"); for(i=1;i<=100;i+) /計算sum_num、xishu、s_point(start point）、d_point sum_num=pow(2,i-1); xishu=double(a-b)/sum_num; s_point=double(b)+double(a-b)/pow(2,i); d_point=dou

26、ble(a-b)/pow(2,i-1); for(j=1;j<=sum_num;j+) add_T=add_T+f_x(s_point+(j-1)*d_point); add_T=add_T*xishu; T1=(T0+add_T)/2; err_T=(T1-T0)/3; /output printf("%d %d %10.8f %10.8f",i,pow(2,i),T1,err_T); printf("n"); if(err_T<=0) err_T=(-1)*err_T; if(err_T<=E) break; else T0=T1;

27、 T1=0; add_T=0; err_T=0; /result printf("n"); printf("T1=%9.7f",T1); printf("n"); printf("n");2.復化simpose源代碼#include <stdio.h>#include <math.h>#define a 5#define b 1#define E 0.00000005 /即保留七位有效數字0.5*10-7/原函數double f_x(double x) double y; y=exp(-(x

28、*x); return(y);int pow(int x,int y) int z=1; int i; for(i=0;i<y;i+) z=z*x; return z;void main() /計算T1,T2,T4,T8. double T0=0,T1=0; double err_T=10; int i=0,j=0; int sum_num=0; double xishu=0; double s_point=0,d_point=0,sd_point=0; double add_T=0; T0=(a-b)*(f_x(a)+4*f_x(a+b)/2)+f_x(b)/6;/n=2的cotes公

29、式即simpson公式 printf("n"); printf("=數值積分_利用復化simpson公式=n"); printf("n"); printf("i 2i T_2i (T_2i-T_2(i-1)/3 n"); printf("n"); printf("0 %d %10.8f 0",pow(2,0),T0); printf("n"); for(i=1;i<=100;i+) /計算sum_num、xishu、s_point(start poi

30、nt）、d_point /long powl (long double x, long double y) sum_num=pow(2,i-1); /the same as T xishu=double(a-b)/sum_num/6; s_point=double(b)+double(a-b)/pow(2,i); d_point=double(a-b)/pow(2,i-1); sd_point=double(a-b)/sum_num/4; for(j=1;j<=sum_num;j+) add_T=add_T+2*f_x(s_point+(j-1)*d_point)-4*f_x(s_poi

31、nt-sd_point+(j-1)*d_point)-4*f_x(s_point+sd_point+(j-1)*d_point); add_T=add_T*xishu; T1=(T0-add_T)/2; err_T=(T1-T0)/15; /output printf("%d %d %10.8f %10.8f",i,pow(2,i),T1,err_T); printf("n"); if(err_T<=0) err_T=(-1)*err_T; if(err_T<=E) break; else T0=T1; T1=0; add_T=0; err_

32、T=0; /result printf("n"); printf("T1=%9.7f",T1); printf("n"); printf("n");3.復化cotes源代碼#include <stdio.h>#include <math.h>#define a 5#define b 1#define E 0.00000005 /即保留七位有效數字0.5*10-7/原函數double f_x(double x) double y; y=exp(-(x*x); return(y);int pow

33、(int x,int y) int z=1; int i; for(i=0;i<y;i+) z=z*x; return z;void main() /計算T1,T2,T4,T8. double T0=0,T1=0; double err_T=10; int i=0,j=0; int sum_num=0; double xishu=0; double s_point=0,d_point=0,sd_point=0; double add_T=0; T0=(a-b)*(7*f_x(1)+32*f_x(2)+12*f_x(3)+32*f_x(4)+7*f_x(5)/90;/四階（n=4）cote

34、s公式 printf("n"); printf("=數值積分_利用復化cotes公式=n"); printf("n"); printf("i 2i T_2i (T_2i-T_2(i-1)/3 n"); printf("n"); printf("0 %d %10.8f 0",pow(2,0),T0); printf("n"); for(i=1;i<=100;i+) /計算sum_num、xishu、s_point(start point）、d_poin

35、t /long powl (long double x, long double y) sum_num=pow(2,i-1); /the same as T xishu=double(a-b)/sum_num/90; s_point=double(b)+double(a-b)/pow(2,i); d_point=double(a-b)/pow(2,i-1); sd_point=double(a-b)/sum_num/8; for(j=1;j<=sum_num;j+) add_T=add_T-2*f_x(s_point+(j-1)*d_point)-32*f_x(s_point-sd_po

36、int+(j-1)*d_point)+20*f_x(s_point-2*sd_point+(j-1)*d_point)-32*f_x(s_point-3*sd_point+(j-1)*d_point)-32*f_x(s_point+sd_point+(j-1)*d_point)+20*f_x(s_point+2*sd_point+(j-1)*d_point)-32*f_x(s_point+3*sd_point+(j-1)*d_point); add_T=add_T*xishu; T1=(T0-add_T)/2; err_T=(T1-T0)/63; /output printf("%d

37、 %d %10.8f %10.8f",i,pow(2,i),T1,err_T); printf("n"); if(err_T<=0) err_T=(-1)*err_T; if(err_T<=E) break;else T0=T1; T1=0; add_T=0; err_T=0; /result printf("n"); printf("T1=%9.7f",T1); printf("n"); printf("n");1. romberg源代碼#include <stdi

38、o.h>#include <math.h>#define a 5#define b 1#define E 0.00000005 /即保留七位有效數字0.5*10-7/原函數double f_x(double x) double y; y=exp(-(x*x); return(y);int pow(int x,int y) int z=1; int i; for(i=0;i<y;i+) z=z*x; return z;void main() /計算T1,T2,T4,T8. double T0=0,S0=0,C0=0,T1=0,S1=0,C1=0,R0=0,R1=0; double err_T=10; int i=0,j=0; int sum_num=0; double xishu=0;/系數double s_point=0,d_point=0; double