romberg.cpp

计算方法算法中的romberg算法

#include<iostream.h>
#include<math.h>
#define n 100
typedef double * dd;
double mi(double x,int y);
void romberg(double a,double b,double e);
double f(double x);
void main()
{
	double a,b,e;
	cout<<"输入下限：";
	cin>>a;
	cout<<"输入上限：";
	cin>>b;
	cout<<"输入允许误差：";
	cin>>e;
	romberg(a,b,e);
}

double mi(double x,int y)
{
	double xx=1;
	for(int i=0;i<y;i++)xx*=x;
	return xx;
}

void romberg(double a,double b,double e)
{
	cout<<"f(x)=4/(1+x*x)"<<endl;
	int i=0;
	double sum=0;
	dd * t=new dd[n];
	for(i=0;i<n;i++)t[i]=new double[n];
	t[0][0]=(b-a)*(f(a)+f(b))/2;
	int k=0;
	do
	{
		k++;sum=0;
	    for(i=1;i<=mi(2,k-1);i++)sum+=f(a+(2*i-1)*(b-a)/mi(2,k));
	    t[0][k]=t[0][k-1]/2+(b-a)*sum/mi(2,k);
	    for(i=1;i<=k;i++)
		    t[i][k-i]=(mi(4,i)*t[i-1][k-i+1]-t[i-1][k-i])/(mi(4,i)-1);
	}while(fabs(t[k][0]-t[k-1][0])>=e);
	cout<<"结果是："<<t[k][0]<<endl;
}

double f(double x)
{
	return 4/(1+x*x);
}