GSL----积分部分（翻译）

GSL----积分部分（翻译）

TR:SAN EMIAL:[e-mail] NUAA 2011.3.22

---------------------------------------------------------------

每个算法都是计算积分表达式的近似值：

其中w(x)是权函数，一般积分时取w(x)=1.通过提供绝对精度epsabs和相对精度epsrel，达到预期的计算精度满足：

注意这是不是一个同时满足的约束条件，而是当满足绝对误差或者相对误差或者都满足时计算结束。可以令epsabs=0或epsrel=0,这表明你想计算积分的精确值，此时gsl可能会因为条件太苛刻而计算失败，但是通常情况下gsl会给出一个尽可能满意的结果。

积分函数命名中满足下列简称规则：

Q: quadrature routine

N: non-adaptive integrator 非自适应积分器

A: adaptive integrator 自适应积分器

G: general integrator 常规积分

W: 带权重的积分

S: singularities can be readily integrated 能积分带奇异点

P point of special difficulty can be supplied

I: 积分界限为无穷

0: 带周期性权重函数cos或sin

F: 傅里叶积分

C: 柯西准则求解（ Cauchy principal value）

1. QNG no-adaptive gauss-krondrod integration 非自适应gauss-krondrod积分

int gsl_integration_qng (const gsl_function * f,

double a, double b,

double epsabs, double epsrel,

double *result, double *abserr,

size_t * neval);

这个函数将用10点、21点、43点、87点的gauss-krondrod积分来计算直到误差在允许范围之内，函数返回积分结果result、使用的积分点数neval、绝对误差值估计值abserr。a,b是积分上下限，epsabs为绝对误差上限，epsrel为相对误差上限。

f是一个结构体，它指明要积分的表达式函数。

struct gsl_function_struct

{

double (* function) (double x, void * params);

void * params;

};

其中parms是传递进去的额外参数，可以在目标函数中使用。

目标函数为：

double fx(double x, void * params);

{

Return ….

}

下面这个例子是计算

f(x)=sin(x)/x在1-2之间的积分值：

#include <iostream>
#include <gsl/gsl_integration.h>
usingnamespace std;
//[e-mail]
double fx(double x,void*params)
{
return sin(x)/x;
}
void main()
{

gsl_function f;
f.function=fx;

double r,er;
unsigned int n;
gsl_integration_qng(&f,
1,2,1e-10,1e-10,&r,&er,&n);
cout<<"result="<<r<<endl<<"abserr="<<er<<endl<<"neval="<<n<<endl;
}

------------------------------------------------------------------------------

2. QAG 自适应积分

QAG是一般的自适应积分算法。该算法思路是将积分区域分为若干个子区域，在在子区间中计算各个积分的近似误差，其中最大误差所在区间将被细分为若干个子区间，不断细分直到满足误差要求。该算法通过区域划分来子区间难积分点的误差，从而提高了计算精度。这些子区间有gsl_integration_workspace管理，它将负责管理内存和获取计算结果。

gsl_integration_workspace * gsl_integration_workspace_alloc (const size_t n);

该函数申请n个存放double型区间和它们的积分结果、近似误差的内存空间。

Void gsl_integration_workspace_free (gsl_integration_workspace * w);

该函数将释放w的内存空间。

QAG积分函数：

int gsl_integration_qag (const gsl_function * f,

double a, double b,

double epsabs, double epsrel, size_t limit,

int key,

gsl_integration_workspace * workspace,

double *result, double *abserr);

f:目标函数

a,b:积分区间

epsabs：积分绝对误差上限

epsrel：积分相对误差上限

limit： 子区间划分上限,不能超过gsl_integration_workspace申请的区间数目

key： 高斯积分点选项

GSL_INTEG_GAUSS15 = 1, /* 15 point Gauss-Kronrod rule */

GSL_INTEG_GAUSS21 = 2, /* 21 point Gauss-Kronrod rule */

GSL_INTEG_GAUSS31 = 3, /* 31 point Gauss-Kronrod rule */

GSL_INTEG_GAUSS41 = 4, /* 41 point Gauss-Kronrod rule */

GSL_INTEG_GAUSS51 = 5, /* 51 point Gauss-Kronrod rule */

GSL_INTEG_GAUSS61 = 6 /* 61 point Gauss-Kronrod rule */

Result：积分结果

Abserr：近似绝对误差

下面这个例子将计算

exp( sin(x)/x ) 在-1,2之间的积分值：

#include <iostream>
#include <gsl/gsl_integration.h>
usingnamespace std;
//[e-mail] SAN NUAA 2011
double fx(double x,void*params)
{
return exp( sin(x)/x );
}
void main()
{

gsl_function f;
f.function=fx;

double r,er;
unsigned int n;

gsl_integration_workspace*w=gsl_integration_workspace_alloc(10);
if(0==gsl_integration_qag
(&f,//积分函数
-1,2,//积分上下限
1e-10,//绝对误差上限
1e-10,//相对误差上限
10,//子区间数目上限
GSL_INTEG_GAUSS41,//gauss积分积分点选取
w,//内存管理
&r,//积分结果
&er//误差
))
cout<<"result="<<r<<endl<<"abserr="<<er<<endl<<endl;

gsl_integration_workspace_free(w);
}
//result=7.10276
//matlab=
matlab=7.1028

3.QAGS 带奇异值的自适应积分

处理带区间带奇异点积分是通过在奇异点附近不断地产生自适应区间来实现的。随着子区间在尺寸上的缩小，积分的精度也将得到很好的近似，可以通过外推法（extrapolation）来加速。QAGS结合了自适应区间划分和Wynn epsilon-algorithm 加速在奇异点的积分。该函数采用gauss21点积分方案。

函数如下：

int gsl_integration_qags (const gsl_function * f,

double a, double b,

double epsabs, double epsrel, size_t limit,

gsl_integration_workspace * workspace,

double *result, double *abserr);

f:目标函数

a,b:积分区间

epsabs：积分绝对误差上限

epsrel：积分相对误差上限

limit： 子区间划分上限,不能超过gsl_integration_workspace申请的区间数目

Result：积分结果

Abserr：近似绝对误差

下面这个例子将计算

exp( sin(x)/x ) 在-1,2之间的积分值：

#include <iostream>
#include <gsl/gsl_integration.h>
usingnamespace std;
//[e-mail] SAN NUAA 2011
double fx(double x,void*params)
{
return exp( sin(x)/x );
}
void main()
{

gsl_function f;
f.function=fx;

double r,er;
unsigned int n;

gsl_integration_workspace*w=gsl_integration_workspace_alloc(10);
if(0==gsl_integration_qags
(&f,//积分函数
-1,2,//积分上下限
1e-10,//绝对误差上限
1e-10,//相对误差上限
10,//子区间数目上限
w,//内存管理
&r,//积分结果
&er//误差
))
cout<<"result="<<r<<endl<<"abserr="<<er<<endl<<endl;

gsl_integration_workspace_free(w);
}
//result=7.10276
//matlab=7.1028

4.QAGP 考虑已知奇异点的自适应积分

这个积分函数是QAGS的延伸，它考已知的奇异积分点，通过参数double *pts来标示奇异点。比如在（a,b）之间有三个奇异积分点a<x1<x2<x3<b.则

Pts[0]=a;Pts[1]=x1;pts[2]=x2;pts[3]=x3; npts=3; Pts[4]=b;

Npts=5;

如果知道积分区域的奇异值点，函数QAGP会比QAGS快很多。

int gsl_integration_qagp (const gsl_function * f,

double *pts, size_t npts,

double epsabs, double epsrel, size_t limit,

gsl_integration_workspace * workspace,

double *result, double *abserr);

f:目标函数

pts:积分区间，以及奇异点，如Pts[0]=a;Pts[1]=x1;pts[2]=x2;pts[3]=x3; npts=3; Pts[4]=b;

npts：pts中的点数

epsabs：积分绝对误差上限

epsrel：积分相对误差上限

limit： 子区间划分上限,不能超过gsl_integration_workspace申请的区间数目

workspace：负责区间维护的gsl_integration_workspace

Result：积分结果

Abserr：近似绝对误差

--------------------------------------------------------------------------------

5.QAGI积分区间为∞的自适应积分

int gsl_integration_qagi (gsl_function * f,

double epsabs, double epsrel, size_t limit,

gsl_integration_workspace * workspace,

double *result, double *abserr);

这个函数通过适当的变换将函数的积分区间有（-∞，∞）转换为区间（0,1】来计算，引入变量x=(1-t)/t, t=0~1,则积分公式为：

然后用15点高斯积分公式QAGS来计算，因为变换会在原点产生积分奇异点，因此低阶高斯积分会更有效果。

下面这个例子将计算无穷积分：

1. exp(-x*x) （ans=pi^0.5)

2. 1/(x*x+1) (ans=pi)

#include <iostream>
#include <time.h>
#include <gsl/gsl_integration.h>
usingnamespace std;
//[e-mail] SAN NUAA 2011

//exp(-x*x)
double fx(double x,void*params)
{
return exp(-x*x);
}
//1/(x*x+1)
double fx1(double x,void*params)
{
return1/(x*x+1);
}
void main()
{

gsl_function f;
f.function=fx;

double r,er;
unsigned int n;

clock_t t1,t2;
t1=clock();

gsl_integration_workspace*w=gsl_integration_workspace_alloc(1000);

//exp(-x*x) (-∞,+∞)
if(0==gsl_integration_qagi
(&f,//积分函数
1e-10,//绝对误差上限
1e-10,//相对误差上限
1000,//子区间数目上限
w,//内存管理
&r,//积分结果
&er//误差
))
cout<<"[ exp(-x*x) (-∞,+∞) ]="<<r<<endl;
//1/(x*x+1) (-∞,+∞)
f.function=fx1;
if(0==gsl_integration_qagi
(&f,//积分函数
1e-10,//绝对误差上限
1e-10,//相对误差上限
1000,//子区间数目上限
w,//内存管理
&r,//积分结果
&er//误差
))
cout<<"[ 1/(x*x+1) (-∞,+∞) ]="<<r<<endl;

gsl_integration_workspace_free(w);
}
//result=1.77245 3.15159
/*matlab cmd
syms x;
int(exp(-x*x),-inf,inf); -> pi^0.5
evalf(ans) -> 1.772453850905516
int(1/(1+x*x),-inf,inf); -> pi
evalf(ans) -> 3.141592653589793

*/

6.积分区间为【a,∞）的自适应积分

int gsl_integration_qagiu (gsl_function * f,

double a,

double epsabs, double epsrel, size_t limit,

gsl_integration_workspace * workspace,

double *result, double *abserr);

这个函数计算【a,∞】的积分值，通过引入变量x=a+(1-t)/t来实现，t=(0,1]，然后调用QAGS自适应积分来计算。变换后公式为：

函数参数意义参考QAGS。

7.积分区间为（-∞，b）的自适应积分

int gsl_integration_qagil (gsl_function * f,

double b,

double epsabs, double epsrel, size_t limit,

gsl_integration_workspace * workspace,

double *result, double *abserr);

这个函数计算【-∞，b】的积分值，通过引入变量x=b-(1-t)/t来实现，t=(0,1]，然后调用QAGS自适应积分来计算。变换后公式为：

8.QAWC 柯西主值(Cauchy Principal Value)自适应积分

int gsl_integration_qawc (gsl_function *f,

const double a, const double b, const double c,

const double epsabs, const double epsrel, const size_t limit,

gsl_integration_workspace * workspace,

double * result, double * abserr);

这个函数计算函数f(x)在积分区间（a,b）上带有一个奇异点c的柯西主值(Cauchy Principal Value)。

这个算法在自适应算法QAG为基础上进行改进，确保子区间不会包含奇异点c。当子区间包含或者接近奇异点c时，算法将使用专门改进的25点Clenshaw-Curtis准则来控制奇异点。在远离奇异点的区域，积分策略是15点gauss-kronrod积分方案。

9.QAWS 奇异函数自适应积分

QAWS是针对包含在积分区域端点奇异对数代数函数。该算法需要定义结构体gsl_integration_qaws_table来计算Chebyshev moments。

该结构体相关函数如下：

gsl_integration_qaws_table *

gsl_integration_qaws_table_alloc (double alpha, double beta, int mu, int nu);

该函数为gsl_integration_qaws_table分配内存空间，相关的奇异权重函数w(x)有四个参数

（alpha, beta, mu, nu），：

其中apha>-1,beta>-1,mu=0、1，nu=0、1。Mu和nu取不同值，权重函数w(x)总共有四种变化：

可以通过函数gsl_integration_qaws_table_set来设置alpha, beta, mu, nu。设置函数定义如下：

Int gsl_integration_qaws_table_set (gsl_integration_qaws_table * t, double alpha, double beta, int mu, int nu);

通过函数gsl_integration_qaws_table_free释放内存。

QAWS是计算函数f(x)在权重w(x)下，在区间（a,b）之间的积分值，函数表达式如下：

该算法基于QAG算法，当积分子区间包含其中某个端点时，算法将采用改进的25点Chenshaw-Curtis准则来控制奇异点；在不包含奇异点的积分子区间，将采用一般的15点gauss-kronrod积分方案。

10.QAWO周期函数自适应积分

QAWO是计算被积包含周期函数因子sin(wx)或cos(wx)的积分值。计算该积分需要结构体gsl_integration_qawo_table，该结构体的维护函数如下：

gsl_integration_qawo_table_alloc (double omega, double L,

enum gsl_integration_qawo_enum sine,size_t n);

这个函数用于申请内存空间，参数omega为w,参数L为被积分函数区间的长度，即L=b-a。sine为GSL_INTEG_COSIN或者 GSL_INTEG_SINE。gsl_integration_qawo_table是一个存储在计算积分过程中所需要的三角函数函数系数，参数n决定了系数项数，每一项都对应着区间L的一个子区间，所有n的取值应满足L/n^2大于子区间数。在级数n不足以满足计算精度时，函数gsl_integration_qawo将返回GSL_ETABLE。

函数gsl_integration_qawo_table_set勇于设置结构体；

函数gsl_integration_qawo_table_set_length设置参数L;

函数gsl_integration_qawo_table_free用于释放内存；

积分函数gsl_integration_qawo用于计算函数f(x)带周期性权重sin(wx)或cos(wx)在区间(a,b)上的的积分值：

int gsl_integration_qawo (gsl_function * f,

const double a,

const double epsabs, const double epsrel,

const size_t limit,

gsl_integration_workspace * workspace,

gsl_integration_qawo_table * wf,

double *result, double *abserr);

参数意义：

：f ->被积函数

：a->积分区间起点，终点b=a+L

：epsabs->绝对误差

：epsrel->相对误差

：limit->子区间上限

：workspace->子区间管理

：wf->存储三角函数系数

：result->结果

：abserr->近似误差

下面这个例子将计算：

f(x)=exp(- 1/x)*x; f(x)sin(2.0x)在区间【0.5,2.0】上的积分值

#include <iostream>
#include <time.h>
#include <gsl/gsl_integration.h>
usingnamespace std;
//[e-mail] SAN NUAA 2011
double fx(double x,void*params)
{
return exp(-1/x)*x;
}
void main()
{

gsl_function f;
f.function=fx;

double r,er;
unsigned int n;

clock_t t1,t2;
t1=clock();

gsl_integration_workspace*w=gsl_integration_workspace_alloc(1000);
gsl_integration_qawo_table * wf=
gsl_integration_qawo_table_alloc(
2.0,//omega
1.5,//L=b-a
GSL_INTEG_SINE,//sin(wx)
20);//n取值不是很明确~

if(0==gsl_integration_qawo
(&f,//积分函数
0.5,
1e-10,//绝对误差上限
1e-10,//相对误差上限
1000,//子区间数目上限
w,//内存管理
wf,
&r,//积分结果
&er//误差
))
cout<<"exp(-1/x)*x*sin(2x) [0.5,2.0] ="<<r<<endl;


gsl_integration_workspace_free(w);
gsl_integration_qawo_table_free(wf);
}
//result=0.063929
/*matlab=0.063928999123728*/

11QAWF傅里叶积分

函数gsl_integration_qawf试图计算函数f在【0，∞】区间上的的傅里叶积分值。

int gsl_integration_qawf (gsl_function * f,

const double a,

const double epsabs,

const size_t limit,

gsl_integration_workspace * workspace,

gsl_integration_workspace * cycle_workspace,

gsl_integration_qawo_table * wf,

double *result, double *abserr);

int gsl_integration_qawf (gsl_function * f,

const double a,

const double epsabs,

const size_t limit,

gsl_integration_workspace * workspace,

gsl_integration_workspace * cycle_workspace,

gsl_integration_qawo_table * wf,

double *result, double *abserr);

参数w通过表wf确定，该计算的计算将通过函数QAWO计算每个子区间的积分值来实现。子区间如下：

该函数通过绝对误差abserr来控制计算精度。具体策略如下：在每个区间Ck，算法试图满足误差：TOLk=Uk*abserr,其中Uk=(1-p)p^(k-1),p=9/10。每项的误差之和就是整体误差abserr。

------------------------------------------------------------------------------------------------------------

返回值说明：

成功计算返回0；

错误代码：

GSL_EMAXITER：子区间数使用情况超过了申请的数目limit>n

GSL_EROUND：不能达到所需精度

GSL_ESING：在积分区间存在不可积分的奇异点或者

GSL_EDIVERGE：积分区域不收敛或者收敛速度太慢而无法计算数值积分

-------------------------------------------------------------------------------------------------------------

TR:SAN EMIAL:[e-mail] NUAA 2011.3.22