Romberg求积法

1、安徽中医药大学题目：Romberg求积法c语言编程姓名：杨撞撞学号：13713042班级：13医软（1）班目录1简介2计算公式3算法描述4程序流程图5算法程序表示6算法结果截图l 1.简介龙贝格求积公式也称为逐次分半加速法。它是在梯形公式、辛普森公式和柯特斯公式之间的关系的基础上，构造出一种加速计算积分的方法。 作为一种外推算法, 它在不增加计算量的前提下提高了误差的精度.在等距基点的情况下，用计算机计算积分值通常都采用把区间逐次分半的方法进行。这样，前一次分割得到的函数值在分半以后仍可被利用,且易于编程。l 2.计算公式梯形公式复化辛普森公式复化科特斯公式 龙贝格求积公式其对应的公式为：T2

2、n=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 龙

3、贝格求积结果 l 4.程序流程图例题：用Romberg方法计算积分I=ò0¹sin（x）/xdx的相关算法流程图表示如下图l 5.算法程序表示#include<stdio.h>#include<math.h>#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)(l

4、og10(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

5、 '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

6、) 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&qu

7、ot;);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(

8、"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");retu

9、rn(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

10、("请输入您选用的求积公式第一个字母（大写）：n");y=getchar();switch(y) /根据输入的大写字母判断所选求积公式 case 'T':printf("请输入分半次数：n"); scanf("%d",&k); printf("请输入积分区间a,b且（a<b）:n"); scanf("%lf,%lf",&a,&b); if(a=0) a=0.0000000001; /输入当积分区间为零时默认为一极小数 printf("-n&

11、quot;); printf("nTxing=%0.7fn",Txing(k);break; /输出计算结果小数点后保留7位 case 'S':printf("请输入分半次数大于等于1的正整数:n"); scanf("%d",&k); printf("请输入积分区间a,b且（a<b）:n"); scanf("%lf,%lf",&a,&b); if(a=0) a=0.0000000001; /输入当积分区间为零时默认为一极小数 printf("

12、;-n"); printf("Simpson=%0.7fn",Simpson(k);break; /输出计算结果小数点后保留7位case 'C':printf("请输入分半次数大于等于2的正整数:n"); scanf("%d",&k); printf("请输入积分区间a,b且（a<b）: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且（a<b）:n"); scanf("%lf,%lf",&a,&b); if(a=0) a=0.0000000001; /输入