romberg.cpp

来自「一些数值计算方面的算法」· C++ 代码 · 共 54 行

//利用龙贝格积分法重新计算
#include <iostream.h>
#include <math.h>
#define EPS 0.0000000001
#define max 20
double f(double x)//定义被积分函数
{
	if(x==0)return (1);
	else return(sin(x)/x);
}
void romberg(double a,double b)//计算函数
{
	double T[max][4]={0},h=1.0,e=1.0+EPS;
	double fnew;
	int i,j,k=1,m;
	T[0][0]=h*(f(a)+f(b))/2.0;
	while((k<max)&&(e>EPS))
	{
		fnew=0;              //计算梯形序列
		for(i=1;i<=(int)(pow(2,k-1));i++)fnew=fnew+f(a+(i-0.5)*h);
		T[k][0]=(T[k-1][0]+h*fnew)/2.0;
		for(m=1;m<=k;m++)  //外推S C R序列
		{
			if(m>3)break;//外推到R序列为止
			T[k][m]=(pow(4,m)*T[k][m-1]-T[k-1][m-1])/(pow(4,m)-1);
		}
		if(k>=4)e=fabs(T[k][3]-T[k-1][3]);
		k=k+1;
		h=h/2.0;
	}
		if(k>max)cout<<"方法失败!"<<endl;
        else
		{
			cout<<"        "<<"T\t\tS\t\tC\t\tR"<<endl;
			for(i=0;i<k;i++)
			{
				cout<<"k="<<i<<"      ";
				for(j=0;j<4;j++)
				{
					if(i>=j)cout<<T[i][j]<<"      ";
				}
				cout<<endl;
			}
		}
}
void main()
{
	double a=0.0,b=1.0;
	romberg(a,b);//调用函数
}