利用复化梯形积分和龙贝格积分法计算积分

 

#define EPS 0.0000001

#define maxd 20

 

double f(double x)

{

         if(x==0)

                   return 1.0;

         else

                   return sin(x)/x;

}

 

void fuhuat(double a,double b)//a,b积分区间

{

         int n=1,i;

         double t_2n,t_n,h=1.0,xnew,fnew,e=1.0+EPS;

         t_n=(f(a)+f(b))/2.0;

         cout<<"k="<<n<<""t"<<t_n<<endl;

         while (e>3*EPS)

         {

                   fnew=0.0;

                   for(i=1;i<=pow(2,n-1);i++)

                   {

                            xnew=a+(i-0.5)*h;

                            fnew=fnew+f(xnew);

                   }

                   t_2n=(t_n+h*fnew)/2.0;

 

                   cout<<"k="<<n<<""t"<<t_2n<<endl;

                  

                   e=fabs(t_2n-t_n);

                   t_n=t_2n;

                   n=n+1;

                   h=h/2.0;

         }

         cout<<"Total nodes="<<pow(2,n)+1<<endl;

}

 

void romberg(double a,double b)//a,b积分区间

{

         double t[maxd][4]={0},h=1.0,e=1+EPS;

         double fnew;

         int i,j,k=1,m;

         t[0][0]=h*(f(a)+f(b))/2.0;

         while((k<maxd)&&(e>EPS))

         {

                   fnew=0;

                   for(i=1;i<=pow(2,k-1);i++)

                            fnew=fnew+f(a+(i-0.5)*h);

                   t[k][0]=(t[k-1][0]+h*fnew)/2.0;

                   for(m=1;m<=k;m++)

                   {

                            if(m>3)

                                     break;

                            t[k][m]=(pow(4,m)*t[k][m-1]-t[k-1][m-1])/(pow(4,m)-1);

                   }

                   if(k>=4)

                            e=fabs(t[k][3]-t[k-1][3]);

                   k++;

                   h=h/2.0;

         }

         if(k>maxd)

                   cout<<"失败！"<<endl;

         else

         {

                   cout<<""tT"<<""t"t"<<"S"<<""t"t"<<"C"<<""t"t"<<"R"<<endl;

                   for(i=0;i<k;i++)

                   {

                            cout<<"k="<<i<<""t";

                            for(j=0;j<4;j++)

                                     if(i>=j)

                                               cout<<t[i][j]<<""t";

                            cout<<endl;

                   }

         }

                  

}

转载于:https://www.cnblogs.com/rosion/archive/2008/10/30/1323305.html