实验三 龙贝格方法
【实验类型】 验证性
【实验学时】 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<iostream.h>
#include<math.h>
# 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<MAXRepeat))
{
p=0.0;
for(k=0;k<n;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<<"请输入积分下限:"<<endl;
cin>>a;
cout<<"请输入积分上限:"<<endl;
cin>>b;
Result=Romberg( a, b, function);
cout<<"龙贝格积分结果:"<<Result<<endl; return 0;
}。