数值积分(论文)

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第六章结果对比分析和体会14参考

2、文献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 复化梯形公式在C

6、语言中的实现复化梯形公式运用的程序如下： 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(s_point+(j-1)*

7、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; 在这个函数中我们将复化梯形公式和积分过程都用计算机语言表示出来。首先我们给出复化梯形公式，进行迭代，直到精确度达到设定要求，算出最后结果。2.3 测试结果用复化梯形有效

8、数字四位求得的结果如下：用复化梯形有效数字七位求得的结果如下：由以上结果可以看出取两个不同的精度相对误差比较小，但计算次数大大的增加，复化梯形公式计算次数多。第三章 复化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(i=1;i=100;i+) /计算sum_num、xishu、s_point(start po

9、int）、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_point-

10、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_T=0; 在这个函数中我们将复化simpose公式和积分过程都用计算机

11、语言表示出来。首先我们给出复化simpose公式，进行迭代，直到精确度达到设定要求，算出最后结果。3.3 测试结果用复化simpose迭代取有效数字四位求得的结果如下：用复化simpose迭代取有效数字七位求得的结果如下：由以上结果可以看出两次不同精度要求的计算可以看出不同精度计算计算次数相差较多。精度为四和七间计算次数相差了三次。第四章 复化cotes公式4.1 复化cotes公式的算法描述复化求积公式当L=4时可得复化梯形公式： =复化cotes公式=4.2 复化cotes公式在C语言中的实现复化cotes公式运用的程序如下：T0=(a-b)*(7*f_x(1)+32*f_x(2)+12*

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_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

13、_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(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

14、 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公式和积分过程都用计算机语言表示出来。首先我们给出复化cotes公式，进行迭代，直到精确度达到设定要求，算出最后结果。4.3 测试结果用复化cotes有效数字四位求得的结果如下：用复化cotes有效数字七位求得的结果如下：由以上结果可以两次不同精度计算的结果相差相对前面的方法要

15、大，计算次数增加了三次。第五章 Romberg积分法5.1 Romberg积分法的算法描述Romberg方法也称为逐次分半加速法。它是在梯形公式、辛卜生公式和柯特斯公式之间的关系的基础上，构造出一种加速计算积分的方法。作为一种外推算法, 它在不增加计算量的前提下提高了误差的精度。在等距基点的情况下，用计算机计算积分值通常都采用把区间逐次分半的方法进行。这样，前一次分割得到的函数值在分半以后仍可被利用,且易于编程 。对区间a， b，令h=b-a构造梯形值序列T2K。 T1=hf(a)+f(b)/2 把区间二等分，每个小区间长度为 h/2=(b-a)/2，于是 T2 =T1/2+h/22f(a+h

16、/2）把区间四(22)等分，每个小区间长度为h/22 =（b-a）/4，于是 T4 =T2/2+h/2f（a+h/4）+f（a+3h/4).把a，b2k等分，分点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,C1=0,R0

17、=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_point=dou

18、ble(a-b)/pow(2,i-1); for(j=1;j=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_T0) err_T=(-1)*err_T; if(err_T=4) break; /完成计算后，准备下一次循环 T0=T1; T1=0; S0=S1; S1=0; C0=C1; C1=0; R0=R1; R1=0; 在这个函数中我们将romboerg公式和的积分过程都用计算机语言表示出来。首先我们给出romboerg公式的T0，进行迭代，

19、分别算出S1，C1，R1直到精确度达到设定要求，算出最后结果。5.3 测试结果用romboerg有效数字四位求得的结果如下：用romboerg有效数字七位求得的结果如下：由以上结果可看出，用romboerg取不同的精度对T1，S1，C1，R1的结果影响大小不相同，T1影响最大，R1影响最小，迭代次数越多精度影响大小越小。第六章 结果对比分析和体会通过对不同精度的测试发现复化梯形公式的计算量增加最快，而romberg达到一定的精度要求结果无法正常计算显示。如下图所示当精度要求达到20时结果无法正常显示。而其他可正常显示结果但是计算次数相对较大如复化梯形计算次数为三十三次，由以上程序测试的数据结果

20、的对比显示可知不同求积公式各有特点.梯形求积公式和Simpson求积公式虽然计算简单、使用方便, 但是精度较差, 但对于光滑性较差的被积函数有时比高精度方法更为有效。尤其梯形公式对被积函数是周期函数的效果更为突出。时,复化梯形公式和复化Simpson公式在保留了低阶公式的优点, 又能获得较高的精度, 因此在实际计算中应用的最为广泛。利用二分技术得到的Romberg方法的算法简单, 易于编程实现。当节点加密提高积分近似程度时, 前面计算的结果可以为后面所用, 对减少计算量很有好处, 并有比较简单的误差估计, 能得到若干积分序列, 如果在做收敛性控制时, 同时检查各行、各列, 对于不同性态的函数可

21、以用其中最快的收敛序列来逼近积分。参考文献1 李庆扬，王能超，易大义.数值分析M.武汉.华中科技大学出版社，2006.7.2 清华大学、北京大学计算方法编写组.计算方法M .北京.科学出版社，19803 吕同斌复化梯形公式及其应用期刊论文 安徽水利水电职业技术学院学报 2002 年4 期4 溪梅成数值分析方法 M 合肥：中国科学技术大学，2003附录1.复化梯形源程序#include #include #define a 5#define b 1#define E 0.00000005 /即保留七位有效数字0.5*10-7/原函数double f_x(double x) double y; y

22、=exp(-(x*x); return(y);int pow(int x,int y) int z=1; int i; for(i=0;iy;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; double add_T=0; T0=(a-b)*(f_x(a)+f_x(b)/2;/n=1时的cotes公式即梯形公式 printf(n); pr

23、intf(=数值积分_利用复化梯形公式=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_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)/

24、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; T1=0; add_T=0; err_T=0; /result printf(n); p

25、rintf(T1=%9.7f,T1); printf(n); printf(n);2.复化simpose源代码#include #include #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;iy;i+) z=z*x; return z;void main() /计算T1,T2,T4,T8. double T0=0

26、,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公式即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

27、%10.8f 0,pow(2,0),T0); printf(n); for(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

28、_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_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;

29、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);3.复化cotes源代码#include #include #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

30、i; for(i=0;iy;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）cotes公式 printf(n); printf(=数值积分_利用复

31、化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_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/90; s_point=do

32、uble(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_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_poin

33、t+(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 %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 #include #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);