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

總結

以上是生活随笔為你收集整理的GSL----积分部分（翻译+精简）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。