数值积分(论文)

1、精品文档目录第一章数值积分计算的重述.1.1.1 引言1.1.2 问题重述2.第二章复化梯形公式2.2.1 复化梯形公式的算法描述32.2 复化梯形公式在C语言中的实现32.3 测试结果4.第三章复化simpson公式5.3.1 复化simpson公式的算法描述5.3.2 复化simpson公式在C语言中白实现63.3 测试结果6.第四章复化cotes公式8.4.1 复化cotes公式的算法描述8.4.2 复化cotes公式在C语言中的实现84.3 测试结果9.第五章Romber则分法.105.1 Romberg积分法的算法描述1.05.2 Romberg积分法在C中的实现1.15.3 测试结

2、果12第六章结果对比分析和体会134参考文献1.4附录1.4.精品文档数值积分je，2dx(一)第一章数值积分计算的重述1.1 引言数值积分是积分计算的重要方法，是数值逼近的重要内容，是函数插值的最直接应用，也是工程技术计算中常常遇到的一个问题。在应用上，人们常要求算出具体数值，因此数值积分就成了数值分析的一个重要内容。在更为复杂的计算问题中，数值积分也常常是一个基本组成部分。在微积分理论中，我们知道了牛顿-莱布尼茨(Newton-Leibniz)公式bf(x)dxF(F(aa其中F(x)是被积函数f(x)的某个原函数。但是随着学习的深入，我们发现一个问题：对很多实际问题，上述公式却无能为力。

3、这主要是因为：它们或是被积函数没有解析形式的原函数，或是只知道被积函数在一些点上的值，而不知道函数的形式，对此，牛顿莱布尼茨(Newton-Leibniz)公式就无能为力了。止匕外，即使被积函数存在原函数，但因找原函数很复杂，人们也不愿花费太多的时间在求原函数上，这些都促使人们寻找定积分近似计算方法的研究，特别是有了计算机后，人们希望这种定积分近似计算方法能在计算机上实现，并保证计算结果的精度，具有这种特性的定积分近似计算方法称为数值积分。由定积分知识，定积分只与被积函数和积分区间有关，而在对被积函数做插值逼近时，多项式的次数越高，对被积函数的光滑程度要求也越高，且会出现Runge现象。如n&

4、gt;7时,Newton-Cotes公式就是不稳定的。因而，人们把目标转向积分区间，类似分段插值，把积分区间分割成若干小区问，在每个小区间上使用次数较低的Newton-Cotes公式，然后把每个小区间上的结果加起来作为函数在整个区间上积分的近似，这就是复化的基本思想。本文主要研究的公式有：复化梯形公式、复化Simpson公式、复化Cotes公式、Romberg积分法。1.2 问题重述本文主要介绍微积分方程的复化解法。通过运用复化梯形公式、复化Simpose公式、复化cotes公式和Romberg积分法这四种积分法方法，解出微分方程的近似解。并进行误差分析和结果比较。当积分区间a,b的长度较大，

5、而节点个数n+1固定时，直接使用Newton-Cotes公式的余项将会较大，而如果增加节点个数，即n+1增加时，公式舍舍入误差又很难得到控制，为提高公式的精度，又使算法简单易行，往往使用复化方法。即将积分区间a,b分成若干个子区间，然后在每个小区间上使用低阶Newton-Cotes公式，最后将每个小区间上的积分的近似值相加。将定积分的积分区间ab分割为n等份b-a各下点为xk=a+kh,k=0,1,-nh=在子区间x.xk书(k=0,1,1；.n-1)上n使用NewtonCotes公式将子区间分割为l等份，节点为xk,xk*,xk2h,xk'=xk1记为xk,x1,x2,xl=xk1k

6、-k-k-lll在子区间上作f(x)的l阶Newton-Cotes求积公式af (x)d(x)：口)=(xk 1b广 G(l)f(xkD»尸f(xkPi -0i =0由积分的区间可加性，可得f(x)d(x)cf(x)dxkfxkn 4复化求积公式JI胃k=0n 4 1=h' '、Ci f 区.i) =I n k =0 i =0第二章复化梯形公式2.1 复化梯形公式的算法描述n1nJ1复化求积公式，I胃=卜;1二Cif(Xki)=InkOk=0i=0an11当L=1时可得复化梯形公式：f(x)d(x)%Tn=h££Cif(xkQbk=0i=0n1=h

7、-f(Xk)f(Xk1)k=02b-an1复化梯形公式=f(a)2vf(Xk)f(b)2n«2.2 复化梯形公式在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(startpoint)、d_pointsum_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

8、,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;/outputprintf("%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;elseT0=T1;T1=0;add_T=0;err_T=0;在这个函数中我们将复

9、化梯形公式和积分过程都用计算机语言表示出来。首先我们给出复化梯形公式，进行迭代，直到精确度达到设定要求，算出最后结果。2.3测试结果用复化梯形有效数字四位求得的结果如下：用复化梯形有效数字七位求得的结果如下由以上结果可以看出取两个不同的精度相对误差比较小，但计算次数大大的增加，复化梯形公式计算次数多。第三章复化simpson公式3.1 复化simpson公式的算法描述n1n11复化求积公式、11的大，二Cif(xk.i)=1nkfk=0i=0an41当L=2时可得复化梯形公式：Jf(x)d(x)%Sn=h££Ci(2)f(xk+)bk=0i=0n4 1咤”4%)f(xk 1

10、)n-1复化simpson公式二n-4f(xk1)2'f(xk)f(b)2k=03.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(startpoint)、d_point/longpowl(longdoublex,longdoubley)sum_num=pow(2,i-1);/thesameasTxishu=double(a-b)/sum_num/6;s

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

12、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;elseT0=T1;T1=0;add_T=0;err_T=0;在这个函数中我们将复化simpose公式和积分过程都用计算机语言表示出来。首先我们给出复化simpose公式，进行迭代，直到精确度达到设定要求，算出最后结果。3.3 测试结果用复化simpose迭代取有效数字四位求得的结果如下：ca"C:Document!andn

13、XPCLIENT*而'囊学软件专区DcbugK口口,MM I .1 zzzJ=数值积分/U用复化Simpson公式i2AiT_2Ai010.2455820512013901594380.13937155Tl-0.1393715Pressanykeytocontinue<T_2Ai-T_2A<i-l>>/30-0.00656349-0.000540920.00002371用复化simpose迭代取有效数字七位求得的结果如下:由以上结果可以看出两次不同精度要求的计算可以看出不同精度计算计算次数相差较多。精度为四和七间计算次数相差了三次。第四

14、章复化cotes公式4.1 复化cotes公式的算法描述n_1n_11复化求积公式1胃斗9'Cif(Xk.i)=Ink0k=0i卫an4当L=4时可得复化梯形公式：Jf(x)d(x)生Cn=h££C(4)f(xi)bk=0i=0k-n 二90 k£2复化cotesa bn 1公式二7f(a)90ny7f(xJ32f(xki)12f(xk2)32f(xkj)7f0)44432f(x1)12f(x2)32f(x3)7f(b)k4k4k44.2 复化cotes公式在C语言中的实现复化cotes公式运用的程序如下：T0=(a-b)*(7*f_x(1)+32*f_x

15、(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);/thesameasTxishu=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)-3

16、2*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(s_point+3*sd_point+(j-1)*d_point);add_T=add_T*xishu;T1=(T0-add_T)/2;err_T=(T1-T0)/63;/outputp

17、rintf("%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;elseT0=T1;T1=0;add_T=0;err_T=0;在这个函数中我们将复化cotes公式和积分过程都用计算机语言表示出来。首先我们给出复化cotes公式，进行迭代，直到精确度达到设定要求，算出最后结果。4.3测试结果用复化cotes有效数字四位求得的结果如下：用复化cotes有效数字七位求得的结果如下:-J.n| x|ca&#

18、39;L:DucumentsandSetting号JtPCLIENT桌面,麴学软件专区Debug、匚口pLexe"数值积分/用复化cute将公式i2rT_2Ai<T_2Ai-T_2A<i-l>>/3U10.140566270120.13847502-0.00003319240139402760.000000124160.13940279T1=0.1394028Pressanvkeytocontinue由以上结果可以两次不同精度计算的结果相差相对前面的方法要大，计算次数增加了三次。第五章Romberg积分法5.1 Romberg积分法的

19、算法描述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/2)把区间四(22)等分，每个小区间长度为h/22=(b-a)/4,于是T4=T2/2+h/2f

20、(a+h/4)+f(a+3h/4)把a,b2k等分，分点xi=a+(b-a)/2ki(i=0,1,22k)每个小区间长度为(b-a)/2k,由归纳法可得面所说的的第一个公式.(二)计算公式如下：整个程序就是循着这四个公式进行计算的。Sn,Cn,Rn分别代表特例梯形积分，抛物线积分，龙贝格积分.当然，编程的时候统一处理即可。5.2 Romberg积分法在C中的实现Romberg公式运用的程序如下：doubleT0=0,S0=0,C0=0,T1=0,S1=0,C1=0,R0=0,R1=0;doubleerr_T=10;inti=0,j=0;intsum_num=0;doublexishu=0;/系

21、数doubles_point=0,d_point=0;doubleadd_T=0;T0=(a-b)*(f_x(a)+f_x(b)/2;/n=1时的cotes公式即梯形公式。for(i=1;i<=100;i+)/thefirstbasenumbersum_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)*d_po

22、int);-add_T=add_T*xishu;T1=(T0+add_T)/2;add_T=0;/计itsiS1=4叮1/3-T0/3;if(i>=2)C1=16*S1/15-S0/15;if(i>=3)R1=64*C1/63-C0/63;/checkusingthe"1"dataerr_T=(R1-R0)/255;if(err_T<0)err_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

23、=0;在这个函数中我们将romboerg公式和的积分过程都用计算机语言表示出来。首先我们给出romboerg公式的T0,进行迭代，分别算出S1,C1,R1直到精确度达到设定要求，算出最后结果。5.3 测试结果用romboerg有效数字四位求得的结果如下：r"r:nnrun-ipntsanriSettinqC.MPClIENT'卓面'救学软件专区tDEbug'CpfJ1=数值积分加用ru皿加9公式($阶=i2Ai(T_2i-T_2*<i-l)>/3918.73575988BIl=B.14324283B1-8.1393715SCl=0.13939525

24、Bl=0,13946986err_Tl=0,0000B38flPressanykeytocontinue-IDI-Xl用romboerg有效数子七位求得的结果如下：c|HC;Doturnenf5andSctHng八XP(J.TFZ卓面,敏中软件专区，6hiiqXpp，数值积分环用nbcrg公式8介2i£T_2&1_2气IT)>/30.7357saes0=0,1403613181=0,13940080Cl=0,13940276Hl-0-13740287err_Tl=0,00000>B3Pressanykeytocontinue由以上结果可看出，用romboerg取不

25、同的精度对T1,S1,C1,R1的结果影响大小不相同，T1影响最大，R1影响最小，迭代次数越多精度影响大小越小。第六章结果对比分析和体会通过对不同精度的测试发现复化梯形公式的计算量增加最快，而romberg达到一定的精度要求结果无法正常计算显示。如下图所示当精度要求达到20时结果无法正常显示。而其他可正常显示结果但是计算次数相对较大如复化梯形计算次数为三十三次，由以上程序测试的数据结果的对比显示可知不同求积公式各有特点.梯形求积公式和Simpson求积公式虽然计算简单、使用方便，但是精度较差，但对于光滑性较差的被积函数有时比高精度方法更为有效。尤其梯形公式对被积函数是周期函数的效果更为突出。n

26、>7时,复化梯形公式和复化Simpson公式在保留了低阶公式的优点，又能获得较高的精度，因此在实际计算中应用的最为广泛。利用二分技术得到的Romberg方法的算法简单,易于编程实现。当节点加密提高积分近似程度时，前面计算的结果可以为后面所用，对减少计算量很有好处，并有比较简单的误差估计，能得到若干积分序列，如果在做收敛性控制时，同时检查各行、各列，对于不同性态的函数可以用其中最快的收敛序列来逼近积分。参考文献1李庆扬，王能超，易大义.数值分析M.武汉.华中科技大学出版社，2006.7.2清华大学、北京大学计算方法编写组.计算方法M.北京.科学出版社，19803吕同斌复化梯形公式及其应用期

27、刊论文安徽水利水电职业技术学院学报2002年4期4溪梅成数值分析方法M合肥：中国科学技术大学，2003附录1 .复化梯形源程序#include<stdio.h>#include<math.h>#definea5#defineb1#defineE0.00000005/即保留七位有效数字0.5*10A-7/原函数doublef_x(doublex).doubley;y=exp(-(x*x);return(y);)intpow(intx,inty)(intz=1;inti;for(i=0;i<y;i+)(z=z*x;)returnz;)voidmain()(/计算T1,

28、T2.T4.T8doubleT0=0,T1=0;doubleerr_T=10;inti=0j=0;intsum_num=0;doublexishu=0;doubles_point=0,d_point=0;doubleadd_T=0;T0=(a-b)*(f_x(a)+f_x(b)/2;/n=1时的cotes公式即梯形公式printf("n");printf("=数值积分一利用复化梯形公式=n)；printf("n");printf("i2AiT_2Ai(T_2Ai-T_2A(i-1)/3n");printf("n&q

29、uot;);printff'O%d%10.8f0",pow(2,0),T0);printf("n");for(i=1;i<=100;i+)计算sum_num>xishu、s_point(startpoint)、d_pointsum_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_poin

30、t+(j-1)*d_point);-add_T=add_T*xishu;T1=(T0+add_T)/2;err_T=(T1-T0)/3;/outputprintf("%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;elseT0=T1;T1=0;add_T=0;err_T=0;./resultprintf("n");printf("T1=%9.7f",T1);

31、printf("n");printf("n");2 .复化simpose源代码#include<stdio.h>#include<math.h>#definea5#defineb1#defineE0.00000005/即保留七位有效数字0.5*10A-7/原函数doublef_x(doublex)(.doubley;y=exp(-(x*x);return(y);)intpow(intx,inty)(intz=1;inti;for(i=0;i<y;i+)(z=z*x;)returnz;)voidmain()(/计算T1,T2,

32、T4,T8doubleT0=0,T1=0;doubleerr_T=10;inti=0,j=0;intsum_num=0;doublexishu=0;doubles_point=0,d_point=0,sd_point=0;doubleadd_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("i2AiT_2Ai(T

33、_2Ai-T_2A(i-1)/3n");printf("n");printf("0%d%10.8f0",pow(2,0),T0);printf("n");for(i=1;i<=100;i+)(计算sum_num>xishu、s_point(startpoint)、d_point/longpowl(longdoublex,longdoubley)sum_num=pow(2,i-1);/thesameasTxishu=double(a-b)/sum_num/6;s_point=double(b)+double(a-b

34、)/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-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;/outputprintf("%d%d%10.8

35、f%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;elseT0=T1;T1=0;add_T=0;err_T=0;./resultprintf("n");printf("T1=%9.7f",T1);printf("n");printf("n");3 .复化cotes源代码#include<stdio.h>#include<math.

36、h>#definea5#defineb1#defineE0.00000005/即保留七位有效数字0.5*10A-7/原函数doublef_x(doublex)(doubley;y=exp(-(x*x);return(y);)intpow(intx,inty)(intz=1;inti;for(i=0;i<y;i+)(z=z*x;)returnz;)voidmain()(/计算T1,T2,T4,T8doubleT0=0,T1=0;doubleerr_T=10;inti=0,j=0;intsum_num=0;doublexishu=0;doubles_point=0,d_point=0,

37、sd_point=0;doubleadd_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("=数值积分_禾1用复化cotes公式=n");printf("n");printf("i2AiT_2Ai(T_2Ai-T_2A(i-1)/3n");printf("n");printf("0%d%10.8f0",pow(2,0),T0);p

38、rintf("n");for(i=1;i<=100;i+)(计算sum_num>xishu、s_point(startpoint)、d_point/longpowl(longdoublex,longdoubley)sum_num=pow(2,i-1);/thesameasTxishu=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

39、_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(s_point+3*sd_point+(j-1)*d_point);一一一一add_

40、T=add_T*xishu;T1=(T0-add_T)/2;err_T=(T1-T0)/63;/outputprintf("%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;elseT0=T1;T1=0;add_T=0;err_T=0;./resultprintf("n");printf("T1=%9.7f",T1);printf("n");printf("n");1. romberg源代码#include<stdio.h>#include<math.h># definea5# defineb1# defineE0.00000005/即保留七位有效数字0.5*10A-7/原函数doublef_x(doublex)(.doubley;y=exp(-(x*x);return(y);)intpow(intx,inty)(intz=1;inti;for(i=0;i<y;i+)(z=z*x;)returnz;)vo