Posted on 2008-05-07 15:05
山泉彎延 閱讀(1511)
評論(1) 編輯 收藏 引用 所屬分類:
數(shù)值分析
/*
Romberg Algorithm
開發(fā)者:newplan
開發(fā)日期:08.05.07
*/
/*=======================================*/
/*INCLUDES*/
#include <cstdlib>
#include <iostream>
#include <cmath>
/*=======================================*/
/*MACROS USED IN THIS FILE*/
#define MAX 20
#define PRECISION 0.000008
/*=======================================*/
/*DECLARE NAMES IN STL NAMESAPCE */
using std::cout;
using std::endl;
/*=======================================*/
/*CLASS FUNC (FUNCTION OBJECT): THE ORIGINAL FUNCTION WE WANT TO INTEGRAL*/
class func{
public:
func(double x=1.0):exp(x){}
double operator()(const double& dnum)const{return pow(dnum,exp);}
private:
double exp;
};
/*=======================================*/
/*CLASS ECHELONFUNC (FUNCTION OBJECT)梯形法的遞推公式*/
class echelonFunc{
public:
echelonFunc(double begining,double ending,func & myfunc);
double operator()();
private:
double h;
int n;
double T;
func myfunc;
double begining;
double ending ;
};
/*=======================================*/
echelonFunc::echelonFunc(double begining,double ending,func & myfunc)
{
this->begining=begining;
this->ending=ending;
this->h=ending-begining;
this->n=0;
this->T=0;
this->myfunc=myfunc;//FUCNTION
}
/*------------------------------*/
/* INCREASE FUNCTION 遞推函數(shù)*/
double echelonFunc::operator()()
{ if(this->n==0)
{
this->T=h*0.5*(myfunc(this->begining)+myfunc(this->ending));
this->n=1;
return this->T;
}
double len=0.5*h;
double sum=0;
int k=0;
for(k=0;k<this->n;k++)
{
sum+=myfunc(len);
len=len+h;
}
this->T=0.5*this->T+0.5*h*sum;
this->h/=2;
this->n*=2;
return this->T;
}
/*=======================================*/
/*THE MAIN CLASS IN THIS PROGRAM*/
class Romberg{
public:
Romberg(double begining,double ending,double exp);
~Romberg();
private:
void RombergCPU();/*THE MOST IMPORTANT FUNCTION IN THIS PROGRAM*/
echelonFunc *echol;
double T[MAX][MAX];/*STO THE ROMBERG TABLE*/
};
/*------------------------------*/
Romberg::Romberg(double begining ,double ending ,double paraexp)
{
func myfunc(paraexp);
echol = new echelonFunc(begining,ending,myfunc);
RombergCPU();
}
/*------------------------------*/
Romberg::~Romberg()
{
delete echol;
}
/*------------------------------*/
void Romberg::RombergCPU()
{ clock_t Start; //TIME STARAT
clock_t End; //TIME END
double *p[MAX];//WE USE THIS POINTER ARRAY TO ACCELERATE ALGOTITHM
double **q;
int i,j;
Start = clock();//TIME START FROM HERE
for(i = 0,q = p; q < p+MAX; q++,i++)
*q= &T[i][0];
double a,b,pows;
cout<<"-----------------------Romberg Algorithm---------------------"<<endl;
*p[0]=(*echol)();
cout<<" "<<*p[0]<<endl;
p[0]++;
do{
*p[0]=(*echol)();
cout<<" "<<*p[0];
p[0]++;
for(i=1;;i++)
{
pows=pow(4.0,double(i));
a=pows/(pows-1);
b=1/(pows-1);
*p[i]=a*(*(p[i-1]-1))-b*(*(p[i-1]-2));//ROMBERG ALGORITHM
cout<<" "<<*p[i];
if(p[i]==&T[i][0])
{
p[i]++;
break;
}
p[i]++;
}
cout<<endl;//fabs(T[i][0]-T[i-1][0])
}while(fabs(T[i][0]-T[i-1][0])>PRECISION);
End = clock();//TIME END HERE
cout<<"-------------------------------------------------------------"<<
endl<<" TIME SPEND:"<<(double)(End-Start)/CLOCKS_PER_SEC<<endl;
}
/*=======================================*/
/*MAIN FUNCTION*/
int main(int argc, char *argv[])
{
Romberg romberg(0,1,1.5);//ROMBERG API :BEGIN(0) END(1) EXP(1.5)
system("PAUSE");
return EXIT_SUCCESS;
}