Romberg求积法.doc

精品文档安徽中医药大学题目：Romberg求积法c语言编程姓名：杨撞撞学号：13713042班级：13医软（1）班目录1简介2计算公式3算法描述4程序流程图5算法程序表示6算法结果截图l 1.简介龙贝格求积公式也称为逐次分半加速法。它是在梯形公式、辛普森公式和柯特斯公式之间的关系的基础上，构造出一种加速计算积分的方法。 作为一种外推算法, 它在不增加计算量的前提下提高了误差的精度.在等距基点的情况下，用计算机计算积分值通常都采用把区间逐次分半的方法进行。这样，前一次分割得到的函数值在分半以后仍可被利用,且易于编程。l 2.计算公式梯形公式复化辛普森公式复化科特斯公式 龙贝格求积公式其对应的公式为：T2n=1/2（Tn+Hn) (梯形公式）Sn=4/（4-1）T2n-1(4-1)Tn (辛普森公式)Cn=42/(42-1)S2n-1/(42-1)Sn (柯特斯公式)Rn=43/（43-1）C2n-1/(43-1)Cn (龙贝格求积法公式)l 3.算法描述3.1龙贝格算法基本描述先算出T0（0）,从而计算出T0(1),以此类推，直到计算出|T0(0)-Tn-1（0）|e即可利用加速推算公式推算出结果。 3.2龙贝格算法程序包步骤 1.输入积分上限 2输入程序下限3输入区间等分数4输入要求的函数5计算出所求函数的积分，分别是：u 复化梯形求积结果u 辛普森求积结果u 柯特斯求积结果u 龙贝格求积结果 l 4.程序流程图例题：用Romberg方法计算积分I=0sin（x）/xdx的相关算法流程图表示如下图l 5.算法程序表示#include#include#define A(x)(sin(x)/x) /宏定义若干常用函数A,B,C,D,E,G#define B(x)(cos(x*x+2*x+1)#define C(x)(atan(sqrt(x*x+1)#define D(x)(sqrt(exp(x)+sin(2*x)#define E(x)(x*x*x+3*x*x+5)#define G(x)(log10(x)/pow(2,x)double t20,s20,c20,r20;/定义全局数组double dh,fan,a,b,m; /定义全局变量int jj=0;char hs;double F(double x) /用switch调用若干被积函数switch(hs)case A:fan=A(x);break;case B:fan=B(x);break;case C:fan=C(x);break;case D:fan=D(x);break;case E:fan=E(x);break;case G:fan=G(x);break;default :printf(输入错误！);return(fan);/返回被积函数值double H(int i) /求和函数并返回和SUMint j; double zh,SUM=0.0; /定义求和变量SUM并赋初值for(j=1;j=pow(2,i-1);j+)zh=(a+(2*j-1)*(b-a)/pow(2,i); SUM=SUM+F(zh); /调用F(x)函数 SUM=(b-a)*SUM/pow(2,i);return(SUM);double Txing(int k) /梯形公式if(k=0) dh=tjj=(b-a)/2)*(F(a)+F(b); /分半次数为零时T形公式求积else dh=0.5*Txing(k-1)+H(k); /Txing函数递归调用循环输出并返回dh t+jj=dh; m=pow(2,jj);printf(T%0.0lf=%0.7lft,m,tjj);/输出并返回dhreturn(dh);double Simpson(int k) /辛普森公式int i,j; Txing(k); /调用梯形公式for(i=0;i=k-1;i+)si=(4.0*ti+1-ti)/3.0; /递推辛普森公式 printf(n);for(j=0;j=k-1;j+) /循环输出 m=pow(2,j);printf(S%0.0lf=%0.7lf,m,sj); printf(t);return(sk-1); /返回最后一个值sk-1double Cotes(int k) /科特斯公式int i,j; Simpson(k);for(i=0;i=k-2;i+)ci=(16.0*si+1-si)/15.0; /递推科特斯公式 printf(n);for(j=0;j=k-2;j+) /循环输出m=pow(2,j);printf(C%0.0lf=%0.7lft,m,cj);return(ck-2); /返回最后一个值ck-2double Romberg(int k) /隆贝格公式int i,j; /调用科特斯公式 Cotes(k); for(i=0;i=k-3;i+) ri=(64.0*ci+1-ci)/63.0; /递推隆贝格公式printf(n);for(j=0;j=k-3;j+) /循环输出 m=pow(2,j); printf(R%0.0lf=%0.7lft,m,rj); printf(n);return(rk-3); /返回最后一个值rk-3main()int k;char y;printf(请从以下公式中选择积分函数：n);printf(A: sin(x)/(x)tB:cos(x2+2x+1)tC:atan(sqrt(x2+1)nn);printf(D:sqrt(ex+sin(2x)tE:x3+3*x2+5tG:log10(x)/pow(2,x)nn);printf(请选择函数F(x)（大写）：n);scanf(%c,&hs);getchar();printf(请输入您选用的求积公式第一个字母（大写）：n);y=getchar();switch(y) /根据输入的大写字母判断所选求积公式 case T:printf(请输入分半次数：n); scanf(%d,&k); printf(请输入积分区间a,b且（ab）:n); scanf(%lf,%lf,&a,&b); if(a=0) a=0.0000000001; /输入当积分区间为零时默认为一极小数 printf(-n); printf(nTxing=%0.7fn,Txing(k);break; /输出计算结果小数点后保留7位 case S:printf(请输入分半次数大于等于1的正整数:n); scanf(%d,&k); printf(请输入积分区间a,b且（ab）:n); scanf(%lf,%lf,&a,&b); if(a=0) a=0.0000000001; /输入当积分区间为零时默认为一极小数 printf(-n); printf(Simpson=%0.7fn,Simpson(k);break; /输出计算结果小数点后保留7位case C:printf(请输入分半次数大于等于2的正整数:n); scanf(%d,&k); printf(请输入积分区间a,b且（ab）:n); scanf(%lf,%lf,&a,&b); if(a=0) a=0.0000000001; /输入当积分区间为零时默认为一极小数 printf(-n); printf(Cotes=%0.7fn,Cotes(k);break; /输出计算结果小数点后保留7位case R:printf(请输入分半次数大于等于3的正整数:n); scanf(%d,&k); printf(请输入积分区间a,b且（ab）:n); scanf(%lf,%lf,&a,&b); if(a=0) a=0.0000000001; /输入当积分区间为零时默认为一极小数 printf(-n); printf