romberg算法求积分.cpp

计算方法程序,包括高斯列主元法,复化simpson,标准四阶龙格库塔,Seidel迭代,Romberg算法

#include<iostream.h>
#include<math.h>
double f(double x)//------函数y=sin(x)/x
{
	double m;
	if(x==0)m=1;
	else m=sin(x)/x;
	return m;
}


void main()
{
a1:	double a,b,e,T[20][20];
	int k,i,j;
	cout<<"输入a=";
	cin>>a;
	cout<<"输入b=";
	cin>>b;
	cout<<"输入误差e:";
	cin>>e;
	T[0][0]=(b-a)/2*(f(a)+f(b));
	for(k=1;k<20;k++)           //公式计算
	{
		T[k][0]=0;
		for(i=1;i<=pow(2,k-1);i++)T[k][0]+=f(a+(2*i-1)*(b-a)/pow(2,k));
		T[k][0]*=(b-a)/pow(2,k);
		T[k][0]+=T[k-1][0]/2;
		for(j=1;j<=k;j++)
			T[k-j][j]=(pow(4,j)*T[k-j+1][j-1]-T[k-j][j-1])/(pow(4,j)-1);
		if(fabs(T[0][k]-T[0][k-1])<e)                  //判断误差
		{
			cout<<"由Romberg算法得\n  sin(x)/x ("<<a<<","<<b<<")积分值="<<T[0][k]<<'\n';break;
		}
	}
	cout<<"                   二维数表"<<endl;
	cout<<"-----------------------------------------------"<<endl;
	for(int s=0;s<=k;s++)
	{
		for(int c=0;c<=s;c++)
		{
			cout<<T[s][c]<<'\t';
		}
		cout<<endl;
	}
	cout<<"-----------------------------------------------"<<endl;
	char out='y';
	cout<<"继续?y/n"<<endl;
	cin>>out;
	if(out!='n') goto a1;
}