数值积分 论文.docx
《数值积分 论文.docx》由会员分享，可在线阅读，更多相关《数值积分 论文.docx（30页珍藏版）》请在冰豆网上搜索。
数值积分论文
第一章数值积分计算的重述1
1.1引言1
1.2问题重述2
第二章复化梯形公式2
2.1复化梯形公式的算法描述3
2.2复化梯形公式在C语言中的实现3
2.3测试结果4
第三章复化simpson公式5
3.1复化simpson公式的算法描述5
3.2复化simpson公式在C语言中的实现6
3.3测试结果7
第四章复化cotes公式8
4.1复化cotes公式的算法描述8
4.2复化cotes公式在C语言中的实现9
4.3测试结果10
第五章Romberg积分法11
5.1Romberg积分法的算法描述11
5.2Romberg积分法在C中的实现12
5.3测试结果13
第六章结果对比分析和体会14
参考文献16
附录17
数值积分
（一）
第一章数值积分计算的重述
1.1引言
数值积分是积分计算的重要方法，是数值逼近的重要内容，是函数插值的最直接应用，也是工程技术计算中常常遇到的一个问题。
在应用上，人们常要求算出具体数值，因此数值积分就成了数值分析的一个重要内容。
在更为复杂的计算问题中，数值积分也常常是一个基本组成部分。
在微积分理论中，我们知道了牛顿-莱布尼茨（Newton-Leibniz公式
其中
是被积函数
的某个原函数。
但是随着学习的深入，我们发现一个问题:
对很多实际问题，上述公式却无能为力。
这主要是因为:
它们或是被积函数没有解析形式的原函数,或是只知道被积函数在一些点上的值,而不知道函数的形式,对此,牛顿—莱布尼茨（Newton-Leibniz公式就无能为力了。
此外,即使被积函数存在原函数,但因找原函数很复杂,人们也不愿花费太多的时间在求原函数上,这些都促使人们寻找定积分近似计算方法的研究,特别是有了计算机后,人们希望这种定积分近似计算方法能在计算机上实现,并保证计算结果的精度,具有这种特性的定积分近似计算方法称为数值积分。
由定积分知识,定积分只与被积函数和积分区间有关,而在对被积函数做插值逼近时,多项式的次数越高,对被积函数的光滑程度要求也越高,且会出现Runge现象。
如
时,Newton-Cotes公式就是不稳定的。
因而,人们把目标转向积分区间,类似分段插值,把积分区间分割成若干小区间,在每个小区间上使用次数较低的Newton-Cotes公式,然后把每个小区间上的结果加起来作为函数在整个区间上积分的近似,这就是复化的基本思想。
本文主要研究的公式有:
复化梯形公式﹑复化Simpson公式﹑复化Cotes公式﹑Romberg积分法。
1.2问题重述
本文主要介绍微积分方程的复化解法。
通过运用复化梯形公式、复化Simpose公式、复化cotes公式和Romberg积分法这四种积分法方法，解出微分方程的近似解。
并进行误差分析和结果比较。
当积分区间[a,b]的长度较大,而节点个数n+1固定时，直接使用Newton-Cotes公式的余项将会较大，而如果增加节点个数,即n+1增加时，公式的舍入误差又很难得到控制，为提高公式的精度,又使算法简单易行,往往使用复化方法。
即将积分区间[a,b]分成若干个子区间，然后在每个小区间上使用低阶Newton-Cotes公式，最后将每个小区间上的积分的近似值相加。
将定积分的积分区间[ab]分割为n等份
各节点为
k=0,1,…n
在子区间
（k=0,1,1…..n-1上使用NewtonCotes公式将子区间分割为l等份,
节点为
记为
在子区间上作f（x的l阶Newton-Cotes求积公式
由积分的区间可加性,可得
复化求积公式
第二章复化梯形公式
2.1复化梯形公式的算法描述
复化求积公式
当L=1时可得复化梯形公式：
=
复化梯形公式=
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_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*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测试结果
用复化梯形有效数字四位求得的结果如下：
用复化梯形有效数字七位求得的结果如下：
由以上结果可以看出取两个不同的精度相对误差比较小，但计算次数大大的增加，复化梯形公式计算次数多。
第三章复化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（startpoint）、d_point
//longpowl（longdoublex,longdoubley
sum_num=pow（2,i-1;//thesameasT
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-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公式和积分过程都用计算机语言表示出来。
首先我们给出复化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*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;//thesameasT
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_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
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有效数字七位求得的结果如下：
由以上结果可以两次不同精度计算的结果相差相对前面的方法要大，计算次数增加了三次。
第五章Romberg积分法
5.1Romberg积分法的算法描述
Romberg方法也称为逐次分半加速法。
它是在梯形公式、辛卜生公式和柯特斯公式之间的关系的基础上，构造出一种加速计算积分的方法。
作为一种外推算法,它在不增加计算量的前提下提高了误差的精度。
在等距基点的情况下，用计算机计算积分值通常都采用把区间逐次分半的方法进行。
这样，前一次分割得到的函数值在分半以后仍可被利用,且易于编程。
对区间[a，b]，令h=b-a构造梯形值序列{T2K}。
T1=h[f（a+f（b]/2
把区间二等分，每个小区间长度为h/2=（b-a/2，于是
T2=T1/2+[h/22]f（a+h/2）
把区间四（22等分，每个小区间长度为h/22=（b-a）/4，于是
T4=T2/2+[h/2][f（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.2Romberg积分法在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;//系数
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++//thefirstbasenumber
{
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*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;
//checkusingthe"1"data
err_T=（R1-R0/255;
if（err_T<0
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，进行迭代，分别算出S1，C1，R1直到精确度达到设定要求，算出最后结果。
5.3测试结果
用romboerg有效数字四位求得的结果如下：
用romboerg有效数字七位求得的结果如下：
由以上结果可看出，用romboerg取不同的精度对T1，S1，C1，R1的结果影响大小不相同，T1影响最大，R1影响最小，迭代次数越多精度影响大小越小。
第六章结果对比分析和体会
通过对不同精度的测试发现复化梯形公式的计算量增加最快，而romberg达到一定的精度要求结果无法正常计算显示。
如下图所示当精度要求达到20时结果无法正常显示。
而其他可正常显示结果但是计算次数相对较大如复化梯形计算次数为三十三次，
由以上程序测试的数据结果的对比显示可知不同求积公式各有特点.梯形求积公式和Simpson求积公式虽然计算简单、使用方便,但是精度较差,但对于光滑性较差的被积函数有时比高精度方法更为有效。
尤其梯形公式对被积函数是周期函数的效果更为突出。
时,复化梯形公式和复化Simpson公式在保留了低阶公式的优点,又能获得较高的精度,因此在实际计算中应用的最为广泛。
利用二分技术得到的Romberg方法的算法简单,易于编程实现。
当节点加密提高积分近似程度时,前面计算的结果可以为后面所用,对减少计算量很有好处,并有比较简单的误差估计,能得到若干积分序列,如果在做收敛性控制时,同时检查各行、各列,对于不同性态的函数可以用其中最快的收敛序列来逼近积分。
参考文献
[1]李庆扬，王能超，易大义.数值分析[M].武汉.华中科技大学出版社，2006.7.
[2]清华大学、北京大学计算方法编写组.计算方法[M].北京.科学出版社，1980
[3]吕同斌《复化梯形公式及其应用》[期刊论文]《安徽水利水电职业技术学院学报》2002年4期
[4]溪梅成《数值分析方法》[M]合肥：
中国科学技术大学，2003
附录
1.复化梯形源程序
#include
#include
#definea5
#defineb1
#defineE0.00000005//即保留七位有效数字0.5*10^-7
//原函数
doublef_x（doublex
{
doubley;
y=exp（-（x*x;
return（y;
}
intpow（intx,inty
{
intz=1;
inti;
for（i=0;i
{
z=z*x;
}
returnz;
}
voidmain（
{
//计算T1,T2,T4,T8......
doubleT0=0,T1=0;
doubleerr_T=10;
inti=0,j=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（"i2^iT_2^i（T_2^i-T_2^（i-1/3\n";
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
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*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";
printf（"T1=%9.7f",T1;
printf（"\n";
printf（"\n";
}
2.复化simpose源代码
#include
#include
#definea5
#defineb1
#defineE0.00000005//即保留七位有效数字0.5*10^-7
//原函数
doublef_x（doublex
{
doubley;
y=exp（-（x*x;
return（y;
}
intpow（intx,inty
{
intz=1;
inti;
for（i=0;i
{
z=z*x;
}
returnz;
}
voidmain（
{
//计算T1,T2,T4,T8......
doubleT0=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（"i2^iT_2^i（T_2^i-T_2^（i-1/3\n";
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;//thesameasT
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-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;
}
}
//result
printf（"\n";
printf（"T1=%9.7f",T1;
printf（"\n";
printf（"\n";
}
3.复化cotes源代码
#include
#include
#definea5
#defineb1
#defineE0.00000005//即保留七位有效数字0.5*10^-7
//原函数
doublef_x（doublex
{
doubley;
y=exp（-（x*x;
return（y;
}
intpow（intx,inty
{
intz=1;
inti;
for（i=0;i
{
z=z*x;
}
returnz;
}
voidmain（
{
//计算T1,T2,T4,T8......
doubleT0=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*（7*f_x（1+32*f_x（2+12*f_x（3+32*f_x（4+7*f_
升级会员 