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

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

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

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

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

6、化梯形公式在C語(yǔ)言中的實(shí)現(xiàn)復(fù)化梯形公式運(yùn)用的程序如下： T0=(a-b)*(f_x(a)+f_x(b)/2;/n=1時(shí)的cotes公式即梯形公式for(i=1;i<=100;i+) /計(jì)算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; 在這個(gè)函數(shù)中我們將復(fù)化梯形公式和積分過程都用計(jì)算機(jī)語(yǔ)言表示出來。首先我們給

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

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

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

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

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

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

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

23、其應(yīng)用期刊論文 安徽水利水電職業(yè)技術(shù)學(xué)院學(xué)報(bào) 2002 年4 期4 溪梅成數(shù)值分析方法 M 合肥：中國(guó)科學(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() /計(jì)算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時(shí)的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+) /計(jì)算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() /計(jì)算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+) /計(jì)算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() /計(jì)算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+) /計(jì)算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() /計(jì)算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