數(shù)值積分論文

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

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

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

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

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

6、化梯形公式在C語言中的實現(xiàn)復(fù)化梯形公式運用的程序如下： 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; 在這個函數(shù)中我們將復(fù)化梯形公式和積分過程都用計算機語言表示出來。首先我們給

8、出復(fù)化梯形公式，進行迭代，直到精確度達到設(shè)定要求，算出最后結(jié)果。2.3 測試結(jié)果用復(fù)化梯形有效數(shù)字四位求得的結(jié)果如下：用復(fù)化梯形有效數(shù)字七位求得的結(jié)果如下：由以上結(jié)果可以看出取兩個不同的精度相對誤差比較小，但計算次數(shù)大大的增加，復(fù)化梯形公式計算次數(shù)多。第三章 復(fù)化simpson公式3.1 復(fù)化simpson公式的算法描述復(fù)化求積公式當L=2時可得復(fù)化梯形公式： = 復(fù)化simpson公式= 3.2 復(fù)化simpson公式在C語言中的實現(xiàn)復(fù)化梯形公式運用的程序如下：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; 在這個函數(shù)中我們將復(fù)化simpose公式和積分過程都用計算機語言表示出來。首先我們給出復(fù)化simpose公式，進行迭代，直到精確度達到設(shè)定要求，算出最后結(jié)果。3.3 測試結(jié)果用復(fù)化simpose迭代取有效數(shù)字四位求得的結(jié)果如下：用復(fù)化simpose迭代取有效數(shù)字七位求得的結(jié)果如下：由以上結(jié)果可以看出兩次不同精度要求的計算可以看出不同精度計算計算次數(shù)相差較多。精度為四和七間計算次數(shù)相差了三次。第四章 復(fù)化cotes公式4.1 復(fù)化cotes公式的算法描述復(fù)化求積公式當L=4時可得復(fù)

12、化梯形公式： =復(fù)化cotes公式=4.2 復(fù)化cotes公式在C語言中的實現(xiàn)復(fù)化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; 在這個函數(shù)中我們將復(fù)化cotes公式和積分過程都用計

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

16、區(qū)間逐次分半的方法進行。這樣，前一次分割得到的函數(shù)值在分半以后仍可被利用,且易于編程 。對區(qū)間a， b，令h=b-a構(gòu)造梯形值序列T2K。                 T1=hf(a)+f(b)/2 把區(qū)間二等分，每個小區(qū)間長度為 h/2=(b-a)/2，于是  T2 =T1/2+h/22f(a+h/2）     把區(qū)間四(22)等分，每個小區(qū)間長度為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）每個小區(qū)間長度為（b-a）/ 2k ，由歸納法可得面所說的的第一個公式.（二）計算公式如下：整個程序就是循著這四個公式進行計算的。Sn,Cn, Rn 分別代表特例梯形積分,拋物線積分,龍貝格積分.當然,編程的時候統(tǒng)一處理即可。5.2 Romberg積分法在C中的實現(xiàn)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;/系數(shù)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; /完成計算后，準備下一次循環(huán) T0=T1; T1=0; S0=S1; S1=0; C0=C1; C1=0; R0=R1; R1=0; 在這個函數(shù)中我們將romboerg公式和的積分過程都用計算機語言表示出來。首先我們給出romboerg公式的T0，進行迭代，分別算出S1，C1，R1直到精確度達到設(shè)定要求，算出最后結(jié)果。5.3 測試結(jié)果用romboerg有效數(shù)字四位求得的結(jié)果如下：用romboerg有效數(shù)字七位求得的結(jié)果如下：由以上結(jié)果可看出，用romboerg取不同的精度對T1，S1

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

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

23、其應(yīng)用期刊論文 安徽水利水電職業(yè)技術(shù)學(xué)院學(xué)報 2002 年4 期4 溪梅成數(shù)值分析方法 M 合肥：中國科學(xué)技術(shù)大學(xué)，2003附錄1.復(fù)化梯形源程序#include <stdio.h>#include <math.h>#define a 5#define b 1#define E 0.00000005 /即保留七位有效數(shù)字0.5*10-7/原函數(shù)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("=數(shù)值積分_利用復(fù)化梯形公式=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.復(fù)化simpose源代碼#include <stdio.h>#include <math.h>#define a 5#define b 1#define E 0.00000005 /即保留七位有效數(shù)字0.5*10-7/原函數(shù)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("=數(shù)值積分_利用復(fù)化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.復(fù)化cotes源代碼#include <stdio.h>#include <math.h>#define a 5#define b 1#define E 0.00000005 /即保留七位有效數(shù)字0.5*10-7/原函數(shù)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("=數(shù)值積分_利用復(fù)化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 /即保留七位有效數(shù)字0.5*10-7/原函數(shù)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;/系數(shù)double s_point=0,d_point=0; double