龙贝格方法.doc

实验三 龙贝格方法
【实验内容】1. 理解龙贝格方法的基本思路
2. 用龙贝格方法设计算法，编程求解一个数值积分的问题。
【实验方法与步骤】（对于必须编写计算机程序的实验，要附上学生自己编写的程序）
1． 龙贝格方法的基本思路
龙贝格方法是在积分区间逐次二分的过程中，通过对梯形之值进行加速处理，从而获得高精度的积分值。
2． 龙贝格方法的算法
步骤1 准备初值 和，用梯形计算公式计算出积分近似值

步骤2 按区间逐次分半计算梯形公式的积分近似值
令
，
计算
，
步骤3 按下面的公式积分（为便于编程，写下列形式）
梯形公式：
辛普生公式：
龙贝格公式：
步骤4 精度控制
当，（为精度）时，终止计算，并取为近似值，否则，将步长折半，转步骤2。

#include<stdio.h>
#include<math.h>
float f(float )
{
return 4/(1+x*x);
}
float Romberg(float a,float b,float( * f)(float),float eps)
{
int n=1,k;
float h=b-a,y,temp;
float T1,T2,S1,S2,C1,C2,R1,R2;
T1=(b-a)/2*(( * f)(a)+( * f)(b));
while(1)
{
temp=0;
for(k=0;k<=n-1;k++)
{ y=a+k*h+h/2;
temp+=( * f)(y);
}
T2=(T1+temp*h)/2;
if(fabs(T2-T1)<eps) return T2;
S2=T2+(T2-T1)/3;
if(n==1){T1=T2;S1=S2;h/=2;n*=2;continue;}
C2=S2+(S2-S1)/15;
if(n==2){C1=C2;T1=T2;S1=S2;h/=2;n*=2;continue;}
R2=C2+(C2-C1)/63;
if(n==4){R1=R2;C1=C2;T1=T2;S1=S2;h/=2;n*=2;continue;}
if(fabs(R2-R1)<eps) return R2;
R1=R2;C1=C2;T1=T2;S1=S2;h/=2;n*=2;
}
}
void main()
{ double eps=5e-6;
printf("R=%f",Romberg(0.0,1.0,f,eps));
}【实验结果】