龙贝格求积分.docx

龙贝格求积分.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）;

}