GSL----积分部分(翻译+精简) SAN
每个算法都是计算积分表达式的近似值:
其中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>
using namespace std;
//visualsan@yahoo.cn
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>
using namespace std;
//visualsan@yahoo.cn 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>
using namespace std;
//visualsan@yahoo.cn 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
转自:http://www.cnblogs.com/JustHaveFun-SAN/archive/2011/03/21/visualsan.html