龙贝格求积分.docx
《龙贝格求积分.docx》由会员分享，可在线阅读，更多相关《龙贝格求积分.docx（14页珍藏版）》请在冰豆网上搜索。
龙贝格求积分
数值计算方法大作业
注:
由朱福利邹素云共同完成
龙贝格求积公式
【简介】
龙贝格求积公式也称为逐次分半加速法。
它是在梯形公式、辛卜生公式和柯特斯公式之间的关系的基础上，构造出一种加速计算积分的方法。
作为一种外推算法,它在不增加计算量的前提下提高了误差的精度.
在等距基点的情况下，用计算机计算积分值通常都采用把区间逐次分半的方法进行。
这样，前一次分割得到的函数值在分半以后仍可被利用,且易于编程。
【算法】
对区间[a，b]，令h=b-a构造梯形值序列{T2K}。
T1=h[f（a）+f（b）]/2
把区间二等分，每个小区间长度为h/2=（b-a）/2，于是
T2=T1/2+[h/2]f（a+h/2）
把区间四
（2）等分，每个小区间长度为h/2=（b-a）/4，于是
T4=T2/2+[h/2][f（a+h/4）+f（a+3h/4）.....................
把[a，b]2等分，分点xi=a+（b-a）/2·i（i=0，1，2···2k）每个小区间长度为（b-a）/2.
例：
I=∫0（4/1+X）dx
解按上述五步计算，此处f（x）=4/（1+x）a=0b=1f（0）=4f
（1）=2
由梯形公式得
T1=1/2[f（0）+f
（1）]=3
计算f（1/2）=16/5用变步长梯形公式得
T2=1/2[T1+f（1/2）]=3.1
由加速公式得
S1=1/3（4T2-T1）=3.133333333
求出f（1/4）f（3/4）进而求得
T4=1/2{T2+1/2[f（1/4）+f（3/4）]}
=3.131176471
S2=1/3（4T4-T2）=3.141568628
C1=1/15（16S2-S1）=3.142117648
计算f（1/8）f（3/8）f（5/8）f（7/8）进而求得
T8=1/2{T4+1/4[f（1/8）+f（3/8）+f（5/8）+f（7/8）]}
=3.138988495
S4=1/3（4T3-T4）=3.141592503
C2=1/15（16S4-S2）=3.141594095
R1=1/63（64C2-C1）=3.141585784
把区间再二分，重复上述步骤算得
T16=3.140941613S8=3.141592652
C4=3.141592662R2=3.141592640
程序步骤
1.输入积分上线
2.输入积分下线
3.输入区间等分数
4.输入要求的函数
5.计算出所求函数的积分，分别是
●复化梯形求积结果
●辛普森求积结果
●科特斯求积结果
●龙贝格求积结果
流程图
源程序
//MyTask.cpp:
Definestheentrypointfortheconsoleapplication.
//
#include"stdafx.h"
#include"math.h"
#include
#defineIS_INT（x）（（x）==int（（x）*10/10））
voidfuHuaTiXing（intk,inta,intb）;//复化梯形求积里面T1跟T2
doublef（doublex）;//原函数sinx/x
doublefuHuaTiXing（intk,doubley,inta,intb）;//复化梯形求积里面T1跟T2
doubleleiJia（inta,intb,inti）;//复化梯形求积里面的累加
voidsimpson（intk）;//k代表要计算的simpson的个数
voidcotes（intk）;//k>=3
voidromberg（intk）;
doubleTTT[512]={0.0};
doubleTT[128]={0.0};
doubleS[64]={0.0};
doubleC[32]={0.0};
doubleR[16]={0.0};
intflag;
intmain（intargc,char*argv[]）
{
inta,b,k;//a积分下线b积分上限k等分数
//k最大支持256等分
int_k;
doubleresult;
intloop_k;
intg_quit=0;
intg_quit1=0;
//欢迎界面
printf（"\t\t****************************************************\n"）;
printf（"\t\t****************************************************\n"）;
printf（"欢迎使用RomBerg方法计算积分\n"）;
printf（"本程序支持最小4等分最大256等分\n"）;
printf（"\t\t****************************************************\n"）;
printf（"\t\t****************************************************\n"）;
printf（"想继续运行吗，继续请按1，退出请按0\n"）;
scanf（"%d",&g_quit）;
if（1==g_quit）
{
while
（1）
{
printf（"pleaseinput3numberswhitchmeansthreeconstent.\n"）;
printf（"pleaseinputthefirstnumberwhitchmeans积分下限\n"）;
getchar（）;
scanf（"%d",&a）;
printf（"pleaseinputthesecondnumberwhitchmeans积分上限\n"）;
getchar（）;
scanf（"%d",&b）;
printf（"pleaseinputthethirdnumberwhitchmeans等分数\n\t\t\t\t\t注意:
支持4等分到256等分\n\t\t\t\t\t注意:
一定是2的整数幂\n"）;
scanf（"%d",&k）;
//
if（4==（int）k）
{
_k=2;
}
elseif（8==（int）k）
{
_k=3;
}
elseif（16==（int）k）
{
_k=4;
}
elseif（32==（int）k）
{
_k=5;
}
elseif（64==（int）k）
{
_k=6;
}
elseif（128==（int）k）
{
_k=7;
}
elseif（256==（int）k）
{
_k=8;
}
else
{
printf（"您输入的等分区间不是2的整数幂,请重新输入!
\n\n"）;
getchar（）;
getchar（）;
getchar（）;
continue;
}
aaaaa:
printf（"pleaseselectthefunctionoff（x）.\n"）;
printf（"ifyouwanttochosefunction1（pow（x,2）*sin（x））pleaseinput1;\n"）;
printf（"ifyouwanttochosefunction2（sin（x）*cos（x））pleaseinput2;\n"）;
printf（"ifyouwanttochosefunction3（sin（x）*sin（x）*x）pleaseinput3;\n"）;
printf（"ifyouwanttochosefunction4（sin（x）/x）pleaseinput4.\n"）;
scanf（"%d",&flag）;
if（flag==1||flag==2||flag==3||flag==4）
{
//null
}
else
{
printf（"Youhaveinputedwrongnum!
\nPleasetryagain!
\n\n\n"）;
gotoaaaaa;
}
fuHuaTiXing（_k,a,b）;
simpson（_k）;
cotes（_k）;
romberg（_k）;
printf（"outputfuHuaTiXing（复化梯形）result.\n"）;
for（loop_k=0;loop_k<=_k;loop_k++）
{
printf（"TT%d=%.14f\n",（int）pow（2,loop_k）,TT[loop_k]）;
}
printf（"outputsimpsonresult.\n"）;
for（loop_k=0;loop_k<=_k-1;loop_k++）
{
printf（"S%d=%.14f\n",（int）pow（2,loop_k）,S[loop_k]）;
}
printf（"outputcotesresult.\n"）;
for（loop_k=0;loop_k<=_k-2;loop_k++）
{
printf（"C%d=%.14f\n",（int）pow（2,loop_k）,S[loop_k]）;
}
printf（"outputrombergresult.\n"）;
for（loop_k=0;loop_k<=_k-3;loop_k++）
{
printf（"R%d=%.14f\n",（int）pow（2,loop_k）,R[loop_k]）;
}
printf（"\n\n\n"）;
printf（"继续计算请按1，退出请按0\n"）;
scanf（"%d",&g_quit1）;
if（0==g_quit1）
{
exit（0）;
}
}
}
else
{
exit（0）;
}
return0;
}
doublef（doublex）//原函数sinx/x
{
doubleresult_f;
switch（flag）
{
case1:
result_f=pow（x,2）*sin（x）;
return（result_f）;
break;
case2:
result_f=sin（x）*cos（x）;
return（result_f）;
break;
case3:
result_f=sin（x）*sin（x）*x;
return（result_f）;
break;
case4:
result_f=sin（x）;
if（x!
=0）
{
result_f=result_f/x;
return（result_f）;
}
elsereturn1.0;
break;
default:
break;
}
}
voidfuHuaTiXing（intk,inta,intb）//复化梯形求积里面T1跟T2
{
intloop_i;
doubleresult_b_a;
result_b_a=double（b-a）/2;
for（loop_i=0;loop_i<=k;loop_i++）//k>=3
{
inttemp;
if（loop_i==0）//kmeanszhishu
{
TT[0]=result_b_a*（f（a）+f（b））;//代表T1
}
else
{
temp=pow（2,loop_i）;
TT[loop_i]=TT[loop_i-1]/2+（（double）（b-a）/（temp））*leiJia（a,b,loop_i）;
}
}
for（loop_i=0;loop_i<10;loop_i++）
{
TTT[loop_i]=TT[loop_i];//把局部数组转存到全局数组中方便下面的访问
}
}
doubleleiJia（inta,intb,inti）//复化梯形求积里面的累加
{
intloop_m;
intloop_i;
doubleresult_leijia=0.0;
inttemp;
doubletemp1,temp2;
switch（i）
{
case1:
loop_m=（int）pow（2,（i-1））;
break;
case2:
loop_m=（int）pow（2,（i-1））;
break;
case3:
loop_m=（int）pow（2,（i-1））;
break;
case4:
loop_m=（int）pow（2,（i-1））;
break;
case5:
loop_m=（int）pow（2,（i-1））;
break;
case6:
loop_m=（int）pow（2,（i-1））;
break;
default:
break;
}
for（loop_i=1;loop_i<=loop_m;loop_i++）
{
temp=pow（2,i）;
temp1=（double）（a+（2*loop_i-1）*（（double）（b-a）/temp））;
//temp2=f（temp1）;
result_leijia=result_leijia+f（temp1）;//第一次temp1=0.5
}
return（result_leijia）;
}
//求simpson
voidsimpson（intk）//k代表要计算的ximpson的个数
{
intloop_i=0,loop_j=0;
for（loop_i=0;loop_i<=k-1;loop_i++）
S[loop_i]=（4*TT[loop_i+1]-TT[loop_i]）/3;
}
voidcotes（intk）
{
intloop_i;
for（loop_i=0;loop_i<=k-2;loop_i++）
C[loop_i]=（16*S[loop_i+1]-S[loop_i]）/15;
}
voidromberg（intk）
{
intloop_i;
for（loop_i=0;loop_i<=k-3;loop_i++）
R[loop_i]=（pow（4,3）*C[loop_i+1]-C[loop_i]）/（pow（4,3）-1）;
}
升级会员 