龙贝格算法.cpp

这是计算方法上的所有算法

#include"iostream.h"
#include"math.h"
#include"iomanip.h"
#define EPS 1e-6
double f(double x)//f为求积函数指针
{
	return 4/(1+x*x);
}

void romberg(double(*f)(double),double a,double b)//a,b分别为积分上下限
 {
	double T,sum,s,m;
   int i,j,k;
  double d[6][6];
  s=0.0;
  for(k=0;k<11;k++)
{
	  if(k==0)
{
		  sum=(b-a)/2*(f(a)+f(b));
d[0][k]=sum;
continue;
}
T=sum;
for(i=1;i<=pow(2,k-1);i++)
{
	m=a+(2*i-1)*(b-a)/pow(2,k);
s=s+f(m);
}
s=s*(b-a)/pow(2,k);
sum=T/2+s;
s=0.0;
if(fabs(sum-s)<=EPS)
break;
d[0][k]=sum;
}
for(i=1;i<=5;i++)
for(j=0;j<6-i;j++)
{
	d[i][j]=(pow(4,i)*d[i-1][j+1]-d[i-1][j])/(pow(4,i)-1);
if(fabs(d[i][j+1]-d[i][j])<=EPS)
break;
}
for(i=1;i<=6;i++)
{ 
	for(j=0;j<i;j++)
cout<<"T("<<j<<","<<i-j-1<<")="<<setiosflags(ios::fixed)<<d[j][i-j-1]<<setw(3);
cout<<endl;
}
cout<<"积分结果 ="<<d[j-1][i-j-1]<<endl;
}
void main()
{
double f(double);
romberg(f,0.0,1.0);
}