高斯积分

最近看书，总是提到高斯积分。自己也知道有这么个方法，但是一直不是很清楚。所以这次尝试写了一下，的确明白了一些。但是，写的不好，代码没有优化，也没有扩展，只能算到三次积分，且最多能取5点。以后有机会希望能拓展，希望我那些做数值算法的同学能帮忙就好了。

如果那位还有更好的，能不能铁给我，万分感谢。

/***************************************************************************
  Copyright (C), ALL, begtostudy. NERC. ZZU.
  File name:   
  Author:begtostudy Version:   Date: 2006.11.15
  Description:   
  Others:        
  Function List: 
  History:       
    1. Date:
       Author:
       Modification:
***************************************************************************/

#include <assert.h>
typedef double (*pFun)(double);
double GaussIntegral(pFun fun,int n,double l=-1,double u=1);
typedef double (*pFun2)(double,double);
double GaussIntegral(pFun2 fun,int n,double l=-1,double u=1,double l2=-1,double u2=1);
typedef double (*pFun3)(double,double,double);
double GaussIntegral
(pFun3 fun,int n,double l=-1,double u=1,double l2=-1,double u2=1,double l3=-1,double u3=1);

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

 /***************************************************************************
  Copyright (C), ALL, begtostudy. NERC. ZZU.
  File name:   
  Author:begtostudy Version:   Date: 2006.11.15
  Description:   
  Others:        
  Function List: 
  History:       
    1. Date:
       Author:
       Modification:
***************************************************************************/
double (*GaussLegendreConffient(int n,double (*p)[2]=NULL))[2]
{
 if (p==NULL)
 {
  p = new double[n][2];//第一列是积分权系数，第二列是积分点坐标
 }

 switch(n) {
 case 1:
  p[0][0] = 2;

  p[0][1] = 0;
 case 2:
  p[0][0] = 1;
  p[1][0] = 1;
  
  p[0][1] = -0.577350269189626;
  p[1][1] = -p[0][1];
  break;
 case 3:
  p[0][0] = 0.555555555555556;
  p[1][0] = 0.888888888888889;
  p[2][0] = p[0][0];
  
  p[0][1] = -0.774596669241483;
  p[1][1] = 0;
  p[2][1] = -p[0][1];

  break;
 case 4:
  p[0][0] = 0.347854845137454;
  p[1][0] = 0.652145154862546;
  p[2][0] = p[1][0];
  p[3][0] = p[0][0];
  
  p[0][1] = -0.861136311594053;
  p[1][1] = -0.339981043584856;
  p[2][1] = -p[1][1];
  p[3][1] = -p[0][1];
  break;
 case 5:
  p[0][0] = 0.236926885056189;
  p[1][0] = 0.478628670499366;
  p[2][0] = 0.568888888888889;
  p[3][0] = p[1][0];
  p[4][0] = p[0][0];
  
  p[0][1] = -0.906179845938664;
  p[1][1] = -0.538469310105683;
  p[2][1] = 0;
  p[3][1] = -p[1][1];
  p[4][1] = -p[0][1];
  break;
 case 6:
  p[0][0] = 0.171324492379170;
  p[1][0] = 0.360761573048139;
  p[2][0] = 0.467913934572691;
  p[3][0] = p[2][0];
  p[4][0] = p[1][0];
  p[5][0] = p[0][0];
  
  p[0][1] = -0.932469514203152;
  p[1][1] = -0.661209386466265;
  p[2][1] = -0.238619186083197;
  p[3][1] = -p[2][1];
  p[4][1] = -p[1][1];
  p[5][1] = -p[0][1];

  break;
 default:
  assert(0);
  delete []p;
  p = NULL;
 }
 
 return p;
}

double GaussIntegral(pFun fun,int n,double l/*=-1*/,double u/*=1*/)
{
 assert(l!=u);

 double GI=0;
 double (*p)[2] = new double[n][2];
 GaussLegendreConffient(n,p);
 int i;
 if (l==-1 && u==1)
 {
  for (i=0;i<n;i++)
  {
   GI += p[i][0]*(*fun)(p[i][1]);
  }
 }
 else
 {
  for (i=0;i<n;i++)
  {
   double ti = (l+u)/2+(u-l)/2*p[i][1];
   GI += p[i][0]*(*fun)(ti);
  }
  GI *= (u-l)/2.0;
 }

 return GI;
}
//----------------------End 1
double GaussIntegral
(pFun2 fun,int n,double l/*=-1*/,double u/*=1*/,double l2/*=-1*/,double u2/*=1*/)
{
 assert(l!=u);

 double GI=0;
 double (*p)[2] = new double[n][2];
 GaussLegendreConffient(n,p);
 int i,j;
 if (l==-1 && u==1 && l2==-1 && u2==1)
 {
  for(j=0;j<n;j++)
  for (i=0;i<n;i++)
   GI += p[i][0]*p[j][0]*(*fun)(p[i][1],p[j][1]);
 }
 else
 {
  for(j=0;j<n;j++)
  {double tj = (l2+u2)/2+(u2-l2)/2*p[j][1];
  for (i=0;i<n;i++)
  {
   double ti = (l+u)/2+(u-l)/2*p[i][1];
   GI += p[i][0]*p[j][0]*(*fun)(ti,tj);
  }
  }
  GI *= (u-l)/2.0;
  GI *= (u2-l2)/2.0;
 }

 return GI;
}
//----------------------End 2
double GaussIntegral
(pFun3 fun,int n,double l/*=-1*/,double u/*=1*/,double l2/*=-1*/,double u2/*=1*/,
 double l3/*=-1*/,double u3/*=1*/)
{
 assert(l!=u);

 double GI=0;
 double (*p)[2] = new double[n][2];
 GaussLegendreConffient(n,p);
 int i,j,k;
 if (l==-1 && u==1 && l2==-1 && u2==1 && l3==-1 && u3==1)
 {
  for(k=0;k<n;k++)
  for(j=0;j<n;j++)
  for (i=0;i<n;i++)
   GI += p[i][0]*p[j][0]*p[k][0]*(*fun)(p[i][1],p[j][1],p[k][1]);
 }
 else
 {
  for(k=0;k<n;k++)
  {double tk = (l3+u3)/2+(u3-l3)/2*p[k][1];
  for(j=0;j<n;j++)
  {double tj = (l2+u2)/2+(u2-l2)/2*p[j][1];
  for (i=0;i<n;i++)
  {
   double ti = (l+u)/2+(u-l)/2*p[i][1];
   GI += p[i][0]*p[j][0]*p[k][0]*(*fun)(ti,tj,tk);
  }
  }
  }
  GI *= (u-l)/2.0;
  GI *= (u2-l2)/2.0;
  GI *= (u3-l3)/2.0;
 }

 return GI;
}
//----------------------End 3

//以下是算例函数

#include <math.h>
double testfun(double x)
{
 return x*x*cos(x);
}

//书上结果0.558608
double testfun2(double x)
{
 return 1/(1+x*x);
}

书上结果0.786885

声明：以上积分公式均已验证，如有问题请您留言，谢谢。

传闻中的begtostudy之编程技术文集

 2007.1.26更新：可以做6点积分了，并且将原来的积分精度提的更高，呵呵。