龙贝格求积分.docx

1、龙贝格求积分数值计算方法大作业 注:由朱福利邹素云共同完成龙贝格求积公式【简介】 龙贝格求积公式也称为逐次分半加速法。它是在梯形公式、辛卜生公式和柯特斯公式之间的关系的基础上，构造出一种加速计算积分的方法。 作为一种外推算法, 它在不增加计算量的前提下提高了误差的精度. 在等距基点的情况下，用计算机计算积分值通常都采用把区间逐次分半的方法进行。这样，前一次分割得到的函数值在分半以后仍可被利用,且易于编程 。 【算法】 对区间a， b，令h=b-a构造梯形值序列T2K。 T1=hf(a)+f(b)/2 把区间二等分，每个小区间长度为 h/2=(b-a)/2，于是 T2 =T1/2+h/2f(a+

2、h/2） 把区间四(2)等分，每个小区间长度为h/2 =（b-a）/4，于是 T4 =T2/2+h/2f（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=0 b=1 f(0)=4 f(1)=2 由梯形公式得 T1=1/2f(0)+f(1)=3 计算f(1/2)=16/5 用变步长梯形公式得 T2=1/2T1+f(1/2)=3.1 由加速公式得 S1=1/3(4T2-T1)=3.133333333

3、求出f(1/4) f(3/4) 进而求得 T4=1/2T2+1/2f(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/2T4+1/4f(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 把区间再二分，重复上述步

4、骤算得 T16=3.140941613 S8=3.141592652 C4=3.141592662 R2=3.141592640 程序步骤1. 输入积分上线2. 输入积分下线3. 输入区间等分数4. 输入要求的函数5. 计算出所求函数的积分，分别是 复化梯形求积结果 辛普森求积结果 科特斯求积结果 龙贝格求积结果流程图源程序/ MyTask.cpp : Defines the entry point for the console application./#include stdafx.h#include math.h#include#define IS_INT(x) ( (x) = int

5、( (x) * 10 / 10 ) void fuHuaTiXing(int k, int a, int b);/复化梯形求积里面 T1跟T2double f(double x);/原函数sinx / xdouble fuHuaTiXing(int k, double y,int a, int b);/复化梯形求积里面 T1跟T2double leiJia(int a, int b, int i);/复化梯形求积里面的累加void simpson(int k);/k代表要计算的simpson的个数void cotes(int k);/k =3void romberg(int k ); doub

6、le TTT512 = 0.0;double TT128 = 0.0;double S64 = 0.0;double C32 = 0.0;double R16 = 0.0;int flag;int main(int argc, char* argv) int a,b,k;/a 积分下线 b 积分上限 k 等分数 /k最大支持256等分 int _k; double result; int loop_k; int g_quit = 0; int g_quit1 = 0; /欢迎界面 printf(tt*n); printf(tt*n); printf( 欢迎使用RomBerg方法计算积分 n);

7、 printf( 本程序支持最小4等分最大256等分 n); printf(tt*n); printf(tt*n); printf(想继续运行吗，继续请按1，退出请按0n); scanf(%d,&g_quit); if(1 = g_quit) while(1) printf(please input 3 numbers whitch means three constent.n); printf(please input the first number whitch means 积分下限n); getchar(); scanf(%d,&a); printf(please input the

8、second number whitch means 积分上限n); getchar(); scanf(%d,&b); printf(please input the third number whitch means 等分数 nttttt注意:支持4等分到256等分 n ttttt注意:一定是2的整数幂n); scanf(%d,&k); / if(4 = (int)k) _k = 2; else if(8 = (int)k) _k = 3; else if(16 = (int)k) _k = 4; else if(32 = (int)k) _k = 5; else if(64 = (int)

9、k) _k = 6; else if(128 = (int)k) _k = 7; else if(256 = (int)k) _k = 8; else printf(您输入的等分区间不是2的整数幂,请重新输入!nn); getchar(); getchar(); getchar(); continue; aaaaa: printf(please select the function of f(x). n); printf(if you want to chose function 1 ( pow(x,2)*sin(x) ) please input 1;n); printf(if you w

10、ant to chose function 2 ( sin(x)*cos(x) ) please input 2;n); printf(if you want to chose function 3 ( sin(x)*sin(x)*x ) please input 3;n); printf(if you want to chose function 4 ( sin(x)/x ) please input 4.n); scanf(%d,&flag); if(flag=1 | flag=2 | flag=3 | flag=4 ) /null else printf(You have inputed

11、 wrong num! nPlease try again!nnn); goto aaaaa; fuHuaTiXing(_k,a,b); simpson(_k); cotes(_k); romberg(_k); printf(output fuHuaTiXing(复化梯形) result.n); for(loop_k=0;loop_k=_k;loop_k+) printf(TT%d =%.14fn,(int)pow(2,loop_k),TTloop_k); printf(output simpson result.n); for(loop_k=0;loop_k=_k-1;loop_k+) pr

12、intf(S%d = %.14fn,(int)pow(2,loop_k),Sloop_k); printf(output cotes result.n); for(loop_k=0;loop_k=_k-2;loop_k+) printf(C%d = %.14f n,(int)pow(2,loop_k),Sloop_k); printf(output romberg result.n); for(loop_k=0;loop_k=_k-3;loop_k+) printf(R%d = %.14fn,(int)pow(2,loop_k),Rloop_k); printf(nnn); printf(继续

13、计算请按1，退出请按0n); scanf(%d,&g_quit1); if(0 = g_quit1) exit(0); else exit(0); return 0;double f(double x)/原函数sinx/x double result_f; switch(flag) case 1: result_f = pow(x,2)*sin(x); return (result_f); break; case 2: result_f = sin(x)*cos(x); return (result_f); break; case 3: result_f = sin(x)*sin(x)*x;

14、return (result_f); break; case 4: result_f = sin(x); if(x != 0) result_f = result_f / x; return (result_f); else return 1.0; break; default :break; void fuHuaTiXing(int k, int a, int b)/复化梯形求积里面 T1跟T2 int loop_i; double result_b_a; result_b_a = double(b - a) / 2; for(loop_i = 0; loop_i=3 int temp; i

15、f(loop_i = 0)/k means zhishu TT0 = result_b_a * ( f(a) + f(b) );/代表T1 else temp = pow(2,loop_i); TTloop_i = TTloop_i - 1 / 2 + ( (double)(b - a) / (temp ) ) * leiJia(a,b,loop_i); for(loop_i=0;loop_i10;loop_i+) TTTloop_i = TTloop_i;/把局部数组转存到全局数组中 方便下面的访问 double leiJia(int a, int b, int i)/复化梯形求积里面的累加

16、 int loop_m ; int loop_i; double result_leijia = 0.0; int temp; double temp1,temp2; switch(i) case 1: loop_m = (int)pow(2,(i-1); break; case 2: loop_m = (int)pow(2,(i-1); break; case 3: loop_m = (int)pow(2,(i-1); break; case 4: loop_m = (int)pow(2,(i-1); break; case 5: loop_m = (int)pow(2,(i-1); bre

17、ak; case 6: 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

18、void simpson(int k)/k代表要计算的ximpson的个数 int loop_i = 0, loop_j = 0; for(loop_i =0;loop_i=k-1;loop_i+) Sloop_i = ( 4 * TTloop_i + 1 -TTloop_i ) / 3; void cotes(int k) int loop_i; for(loop_i=0;loop_i=k-2;loop_i+) Cloop_i = ( 16 * Sloop_i+1 - Sloop_i ) / 15;void romberg(int k ) int loop_i; for(loop_i=0;loop_i=k-3;loop_i+) Rloop_i = ( pow(4,3) * Cloop_i +1 - Cloop_i ) / (pow(4,3) - 1) ;