数值分析龙贝格实验报告

实验三 龙贝格方法

【实验类型】 验证性

【实验学时】 2学时

【实验内容】

1.理解龙贝格方法的基本思路

2.用龙贝格方法设计算法，编程求解一个数值积分的问题。

【实验前的预备知识】

1．计算机基础知识2．熟悉编程基本思想3．熟悉常见数学函数；

【实验方法或步骤】

龙贝格方法的基本思路龙贝格方法是在积分区间逐次二分的过程中，通过

对梯形之值进行加速处理，从而获得高精度的积分值。

1． 龙贝格方法的算法

步骤1 准备初值()f a 和()f b ，用梯形计算公式计算出积分近似值

()()12b a T f a f b -=+⎡⎤⎣

⎦ 步骤2 按区间逐次分半计算梯形公式的积分近似值令

2i b a h -=，0,1,2,...i =计算12102122n n n i i h T T f x -+=⎛⎫=+ ⎪⎝⎭

∑，2i n = 步骤3 按下面的公式积分梯形公式：()223n n n n T T S T -=+

辛普生公式：()2215n n n n S S C S -=+

龙贝格公式：()2263n n n n C C R C -=+

步骤4 精度控制 当2n n R R ε-<，（ε为精度）时，终止计算，并取2n R 为近似值否则将步长折

半，转步骤2。

[实验程序]

#include

#include

# define Precision 0.00001//积分精度要求

# define e 2.71828183

#define MAXRepeat 10 //最大允许重复

double function(double x)//被积函数

{

double s;

s=2*pow(e,-x)/sqrt(3.1415926);

return s;

}

double Romberg(double a,double b,double f(double x))

{

int m,n,k;

double y[MAXRepeat],h,ep,p,xk,s,q;

h=b-a;

y[0]=h*(f(a)+f(b))/2.0;//计算T`1`(h)=1/2(b-a)(f(a)+f(b)); m=1;

n=1;

ep=Precision+1;

while((ep>=Precision)&&(m

{

p=0.0;

for(k=0;k

{

xk=a+(k+0.5)*h; // n-1

p=p+f(xk); //计算∑f(xk+h/2),T

} // k=0

p=(y[0]+h*p)/2.0; //T`m`(h/2),变步长梯形求积公式

s=1.0;

for(k=1;k<=m;k++)

{

s=4.0*s;// pow(4,m)

q=(s*p-y[k-1])/(s-1.0);//[pow(4,m)T`m`(h/2)-T`m`(h)]/[pow(4,m)-1 ],2m阶牛顿柯斯特公式，即龙贝格公式

y[k-1]=p;

p=q;

}

ep=fabs(q-y[m-1]);//前后两步计算结果比较求精度

m=m+1;

y[m-1]=q;

n=n+n; // 2 4 8 16

h=h/2.0;//二倍分割区间

}

return q;

}

main()

{

double a,b,Result;

cout<<"请输入积分下限:"<

cin>>a;

cout<<"请输入积分上限:"<

cin>>b;

Result=Romberg( a, b, function);

cout<<"龙贝格积分结果："<

}