Simpson 积分算法

#define PR double

typedef PR (*FUN)(PR);                                //定义指向被积函数的指针

class Integral
{
private:
 PR lower,upper;                                            //lower——积分下限；upper——积分上限
 FUN integrand;                                              //integrand——被积函数
public:
 Integral(PR l=0,PR u=0,FUN f=0):lower(l),upper(u),integrand(f){}
 void PutChangeBd(PR l,PR u);                  //设置积分上限、下限
 void Put_Integrand(FUN f);                         //设置积分函数
 PR Solution(long n);                                    //用Simpson 3/8法求解积分
};

void Integral::PutChangeBd(PR l,PR u)    //设置积分下限、上限
{
 lower=l;
 upper=u;
}

void Integral::Put_Integrand(FUN f)            //设置积分函数
{
 integrand=f;
}

PR Integral::Solution(long n)                       //用Simpson 3/8法求解积分
{
 PR h,result=0,temp=0;
 h=(upper-lower)/3/n;
 result=integrand(lower)+integrand(upper);
 for (int i=1;i<=n;i++)
  temp+=(integrand(lower+(3*i-1)*h)+integrand(lower+(3*i-2)*h));
 result+=3*temp;
 temp=0;
 for (int i=1;i<n;i++)
  temp+=integrand(lower+3*i*h);
 result+=2*temp;
 result*=(upper-lower)/8/n;
 return result;
}