龙贝格求积分

1、精品龙贝格（ Romberg）求积法1. 算法理论Romberg求积方法是以复化梯形公式为基础，应用Richardson外推法导出的数值求积方法。h 2由复化梯形公式T2 f ( a)2 f (ah2 )f (b)可以化为2T21 h1 f (a) f (b) h1 f (a h2 )221 T1h1 f ( ah2 ) 2一般地，把区间 a,b 逐次分半 k 1 次，（ k 1,2,n ）区间长度（步ba，其中 mk 2k 1。记TkTk(1)长）为 hkmkmk1(1)2由Tk(1)2Tk 1hk 1j 1f (a ( 2 j1)hk )从而bf ( x)dxTk(1)bah2f '

2、; ' ( k )a 12（1）按 Richardson外推思想，可将（ 1 ）看成关于 hk ，误差为 O (hk2 ) 的一个近似公式，因而，复化梯形公式的误差公式为bf ( x)dxTk(1)a24. 2i（2 ） K 1 hkK 2 hki1K i hk1bf (x)dx1a(1)2 i（ 3）取 hk 1 hk 有 Tk 1 2iK i hk 122i 1误差为 O(h2 j ) 的误差公式Tk( j ) Tk( j1) Tk( j 1)Tk( j11)4 j 112. 误差及收敛性分析感谢下载载精品（ 1）误差，对复化梯形公式误差估计时，是估计出每个子区间上的误差，然后将

3、n 个子区间上的误差相加作为整个积分区间上的误差。（ 2）收敛性，记 xih ，由于 Tn ( f )h n f (xi ) f ( xi 1 )2 i01n1n(f ( xi ) xif ( xi ) xi )2i0i 1上面两个累加式都是积分和，由于f ( x) 在区间 a,b 上可积可知，只要 a,b的分划的最大子区间的长度0 时，也即 n时，它们的极限都等于积分值I ( f ) 。可见，只要 f (x) 在区间 a,b 上的 Rieman可积的条件满足， 由复化梯形求积公式计算所得复化梯形值序列Tn ( f ) 都收敛于实际积分值I ( f ) 。3. 算法：（ 1）、输入： a，b

4、，epsilonh1 f (a)f (b)（ 2）、令 h1 b a，计算 T11 21h1（ 3）、令 k2 ， hk 2mk（ 4）、令 mk 2k1，计算 T1T2fajhk1(k 1 ),1h(21)2k 1k )j1（ 5）、forj 2,3,kTkjTk , jTk , j 1Tk 1, j 114 j11（ 6）、ifTk , jTk, j 1输出 Tk , j ;return;（ 7）、k k 1;hk 1 hk 1 ;跳转（ 4）2（ 8）、结束4. 实例感谢下载载精品用龙贝格算法计算积分值I18=2dx0 1 x5. 龙贝格求积源程序：/Romberg method for

5、 Integral/Projectname hu1.cpp/Executable file hu1.exe/Date:2004.12/By Hu Bentao#include "stdafx.h"#include "stdio.h"#include "iostream.h"#include "math.h"#include "conio.h"#include "stdlib.h"#define N1 20#define N2 20double fun(double x)/ 被

6、积函数设置return (8/(1+x*x);感谢下载载精品void Romberg(double a,double b,double(*fun)(double),double epsilon)int i,j,k;double h=b-a,temp;double T2020;FILE *fp;if(fp=fopen("200402157.txt","w")=NULL)/将运算的中间结果和最终结果保存到文件200402157.txtputs("nopen file error!n");return ;fflush(stdin);fpri

7、ntf(fp,"n龙贝格求积计算的中间结果：");printf("n龙贝格求积计算的中间结果：");T11=h*(*fun)(a)+(*fun)(b)/2;fprintf(fp,"nT11=%f",T11);printf("nT11=%f",T11);k=2;h/=2;while(1)temp=0;for(i=1;i<=pow(2,(k-2);i+)感谢下载载精品temp+=(*fun)(a+(2*i-1)*h);Tk1=temp*h+Tk-11/2;fprintf(fp,"nT%d1=%ft&qu

8、ot;,k,Tk1);printf("nT%d1=%ft",k,Tk1);for(j=2;j<=k;j+)Tkj=Tkj-1+(Tkj-1-Tk-1j-1)/(pow(4,j-1)-1);fprintf(fp,"T%d%d=%ft",k,j,Tkj);printf("T%d%d=%ft",k,j,Tkj);j-;if(fabs(Tkj-Tkj-1)<epsilon)fprintf(fp,"n积分结果 tI(f)=");/cout<<endl<<"OK!"<

9、;<endl<<"T"<<k<<j<<"="<<Tkj;fprintf(fp,"%ft",Tkj);/积分值输出printf("n积分结果 tI(f)=");/cout<<endl<<"OK!"<<endl<<"T"<<k<<j<<"="<<Tkj;printf("%ft",T

10、kj);/积分值输出break;k+=1;h/=2;感谢下载载精品fclose(fp);/ hu1.cpp : Defines the entry point for the console application. #include "stdafx.h"#include "stdio.h" #include "iostream.h" #include "math.h" #include "conio.h" #include "stdlib.h"externvoidRomberg(doublea,doubleb,double(*fun)(double),doubleepsilon);extern double fun(double x);int main(int argc, char* argv)double a,b,epsilon;cout<<endl<<"Please input:n积分下限 a="cin>>a;/ 输入积分下限cout<<"积分上限 b="cin>&