龙贝格求-积分.docVIP

龙贝格(Romberg）求积法 1.算法理论 Romberg求积方法是以复化梯形公式为基础，应用Richardson外推法导出的数值求积方法。 由复化梯形公式 可以化为 ＝ 一般地，把区间［a,b］逐次分半k－1次，（k＝1，2,……,n)区间长度（步长)为，其中mk＝2k－1。 记＝ 由＝ 从而 ＝－ (1) 按Richardson外推思想，可将（1）看成关于，误差为的一个近似公式，因而，复化梯形公式的误差公式为 －＝＝ （2） 取＝有 －＝ （3） 误差为的误差公式 ＝＋ 2。误差及收敛性分析 （1）误差，对复化梯形公式误差估计时，是估计出每个子区间上的误差,然后将n个子区间上的误差相加作为整个积分区间上的误差。 (2）收敛性,记，由于 ＝ 上面两个累加式都是积分和，由于在区间上可积可知，只要的分划的最大子区间的长度时,也即时,它们的极限都等于积分值。可见，只要在区间上的可积的条件满足,由复化梯形求积公式计算所得复化梯形值序列都收敛于实际积分值。 3。算法： (1）、输入：a，b，epsilon （2）、令＝b－a，计算＝ （3）、令k＝2，＝ （4)、令＝，计算 （5)、for j＝2，3,……,k （6）、if {输出; return;｝ （7）、k＝k＋1; ＝; 跳转（4) (8)、结束 4.实例 用龙贝格算法计算积分值 I= 5.龙贝格求积源程序： //Romberg method for 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 20 double fun（double x） ｛//被积函数设置 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 T[20］［20]； FILE *fp； if（(fp=fopen(”200402157。txt”,”w”）)==NULL）｛//将运算的中间结果和最终结果保存到文件200402157.txt puts（”\nopen file error！\n）; return ; ｝ fflush（stdin); fprintf(fp,”\n龙贝格求积计算的中间结果：); printf（”\n龙贝格求积计算的中间结果:”); T［1］[1］=h*（（*fun）（a)+（*fun)(b))/2; fprintf(fp,”\nT［1］［1]=％f,T[1］［1]）； printf(”\nT[1］［1］=％f”，T[1］［1]）; k=2；h/=2; while（1）｛ temp=0; for（i=1；i=pow(2，(k-2）)；i++)｛ temp+=(＊fun）(a+（2*i-1)*h); ｝ T[k］[1］=temp*h+T[k—1][1］/2; fprintf（fp，”\nT［％d］[1]=％f\t，k，T［k］［1］）； printf(\nT［%d]［1]=％f\t，k，T［k]［1］)； for（j=2;j=k;j++）｛ T［k］［j］=T[k］[j-1]+(T［k]［j-1］-T［k—1］[j—1]）/(pow（4，j-1）-1); fprintf(fp,”T[%d]［%d]=%f\t,k,j,T［k]［j］）; printf（T[％d][％d]=％f\t，k，j，T[k][j］)； } j——； if((fabs(T［k］［j]—T[k］［j—1]))epsilon)｛ fprintf(fp,”\n积分结果\tI(f)=）； //coutendl〈OK!!”〈endl〈〈T”〈k〈j〈”=T[k]［j］； fprintf（fp,％f\t,T［k][j］)；//积分值输出 printf（”\n积分结果\tI(f)=”）; //cout〈〈endl〈OK！!”〈〈endlT”k〈j〈=”〈T[