Romberg求積法

1、安徽中醫(yī)藥大學題目：Romberg求積法c語言編程姓名：楊撞撞學號：13713042班級：13醫(yī)軟（1）班目錄1簡介2計算公式3算法描述4程序流程圖5算法程序表示6算法結果截圖l 1.簡介龍貝格求積公式也稱為逐次分半加速法。它是在梯形公式、辛普森公式和柯特斯公式之間的關系的基礎上，構造出一種加速計算積分的方法。 作為一種外推算法, 它在不增加計算量的前提下提高了誤差的精度.在等距基點的情況下，用計算機計算積分值通常都采用把區(qū)間逐次分半的方法進行。這樣，前一次分割得到的函數(shù)值在分半以后仍可被利用,且易于編程。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輸入?yún)^(qū)間等分數(shù)4輸入要求的函數(shù)5計算出所求函數(shù)的積分，分別是：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) /宏定義若干常用函數(shù)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;/定義全局數(shù)組double dh,fan,a,b,m; /定義全局變量int jj=0;char hs;double F(double x) /用switch調用若干被積函數(shù)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);/返回被積函數(shù)值double H(int i) /求和函數(shù)并返回和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)函數(shù) 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); /分半次數(shù)為零時T形公式求積else dh=0.5*Txing(k-1)+H(k); /Txing函數(shù)遞歸調用循環(huán)輸出并返回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+) /循環(huán)輸出 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+) /循環(huán)輸出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+) /循環(huán)輸出 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("請從以下公式中選擇積分函數(shù)：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("請選擇函數(shù)F(x)（大寫）：n");scanf("%c",&hs);getchar();printf

10、("請輸入您選用的求積公式第一個字母（大寫）：n");y=getchar();switch(y) /根據(jù)輸入的大寫字母判斷所選求積公式 case 'T':printf("請輸入分半次數(shù)：n"); scanf("%d",&k); printf("請輸入積分區(qū)間a,b且（a<b）:n"); scanf("%lf,%lf",&a,&b); if(a=0) a=0.0000000001; /輸入當積分區(qū)間為零時默認為一極小數(shù) printf("-n&

11、quot;); printf("nTxing=%0.7fn",Txing(k);break; /輸出計算結果小數(shù)點后保留7位 case 'S':printf("請輸入分半次數(shù)大于等于1的正整數(shù):n"); scanf("%d",&k); printf("請輸入積分區(qū)間a,b且（a<b）:n"); scanf("%lf,%lf",&a,&b); if(a=0) a=0.0000000001; /輸入當積分區(qū)間為零時默認為一極小數(shù) printf("

12、;-n"); printf("Simpson=%0.7fn",Simpson(k);break; /輸出計算結果小數(shù)點后保留7位case 'C':printf("請輸入分半次數(shù)大于等于2的正整數(shù):n"); scanf("%d",&k); printf("請輸入積分區(qū)間a,b且（a<b）:n"); scanf("%lf,%lf",&a,&b); if(a=0) a=0.0000000001; /輸入當積分區(qū)間為零時默認為一極小數(shù) printf("-n"); printf("Cotes=%0.7fn",Cotes(k);break; /輸出計算結果小數(shù)點后保留7位case 'R':printf("請輸入分半次數(shù)大于等于3的正整數(shù):n"); scanf("%d",&k); printf("請輸入積分區(qū)間a,b且（a<b）:n"); scanf("%lf,%lf",&a,&b); if(a=0) a=0.0000000001; /輸入