GSL科學計算庫——計算高斯-勒讓德積分

相關文章：

Windows系統Qt5配置GSL科學計算庫

Windows系統下Cygwin+Dev C ++ 配置GSL科學計算庫

假設計算下列積分：

∫ 0 π e x c o s ( x ) d x \int_0^\pi e^xcos(x)dx ∫0π​excos(x)dx

根據高斯勒讓德積分，先将積分區間變換為 [ − 1 , 1 ] [-1,1] [−1,1]，使用線性變換：

x = π − 0 2 t + π + 0 2 x=\frac{\pi-0}{2}t+\frac{\pi+0}{2} x=2π−0​t+2π+0​，将原積分變為

∫ 0 π e x c o s ( x ) d x = − π 2 e π 2 ∫ − 1 1 e π 2 t c o s ( π 2 t ) d x \int_0^\pi e^xcos(x)dx=-\frac{\pi}{2}e^{\frac{\pi}{2}}\int_{-1}^1 e^{\frac{\pi}{2}t}cos(\frac{\pi}{2}t)dx ∫0π​excos(x)dx=−2π​e2π​∫−11​e2π​tcos(2π​t)dx

最後再根據 n n n點求積系數和求積節點計算高斯-勒讓德積分，高斯點和求積系數的選取如下圖：

在GSL科學計算庫的數值積分子產品實作了多種資料積分方法，這裡簡單介紹一下其兩種方法實作的高斯-勒讓德積分，

• double gsl_integration_glfixed(const gsl_function * f, double a, double b,const gsl_integration_glfixed_table * t)

參數解釋：

const gsl_function * f：積分函數；

double a, double b：積分區間；

const gsl_integration_glfixed_table * t：高斯-勒讓德積分的節點數

• int gsl_integration_fixed(const gsl_function * func, double * result, const gsl_integration_fixed_workspace * w)

  This function integrates the function f(x) provided in func using previously computed fixed quadrature rules. The integral is approximated as ∑ i = 1 n ω i f ( x i ) \sum_{i=1}^{n}\omega_if(x_i) i=1∑n​ωi​f(xi​)

  參數解釋：

  const gsl_function * f：積分函數；

  double * result：存放計算結果；

  const gsl_integration_fixed_workspace * w：針對固定點方法積分設定的工作區間，包含積分類型，區間等。

  GSL固定點積分支援的積分類型有：

  

  詳細内容看代碼注釋。

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_sf_gamma.h>

// 需要求積分的函數 f(x)=m*exp(x)*cos(x), 0 <= x <= pi
double f(double x, void * params)
{
	int m = *(int *) params;
	double f = m*exp(x)*cos(x);
	return f;
}
// GSL提供了單獨的高斯-勒讓德積分方法
double gsl_gl(int n)
{
	//	gsl_integration_fixed_workspace * w;
	gsl_integration_glfixed_table * t;

	// 固定的格式類型，選擇勒讓德積分gsl_integration_fixed_legendre

	int m = 1;
//	int n = 2;
	double expected, result;
	gsl_function F;

	t = gsl_integration_glfixed_table_alloc(n);

	F.function = &f;
	F.params = &m;

	result = gsl_integration_glfixed(&F, 0.0,M_PI, t);

	gsl_integration_glfixed_table_free(t);
	return result;
}
/*-----------------------------------------------
 GSL求解積分步驟：
 1）确定積分方式（插值型積分固定點還是其他）以及workshop；
 	該步驟包括了确定積分類型。插值點數，積分區間 等；例：
 	gsl_integration_fixed_workspace * w;
	const gsl_integration_fixed_type * T = gsl_integration_fixed_legendre;
	w = gsl_integration_fixed_alloc(T, n, 0.0, M_PI, 0.0, 0.0);
 2）設定積分類型（ 勒讓德、切比雪夫、埃米爾特等 ）
 3）設定點數（插值型積分，如高斯-勒讓德積分）
 4） 計算積分
 -----------------------------------------------*/
 
 /* --------------------------------------------- 
 Fixed point quadratures支援的積分類型gsl_integration_fixed_type有：
 1.gsl_integration_fixed_legendre
 2.gsl_integration_fixed_chebyshev
 3.gsl_integration_fixed_gegenbauer
 4.gsl_integration_fixed_jacobi 
 5.gsl_integration_fixed_laguerre
 6.gsl_integration_fixed_hermite
 7.gsl_integration_fixed_exponential
 8.gsl_integration_fixed_rational
 9.gsl_integration_fixed_chebyshev2
 -----------------------------------------------*/
 
 
// 下面的例子是計算高斯勒讓德積分
int main (int argc, char *argv[])
{
	
	gsl_integration_fixed_workspace * w;
	const gsl_integration_fixed_type * T = gsl_integration_fixed_legendre;

	int m = 1;
	int n = 2;
	double expected, result;
	gsl_function F;

	w = gsl_integration_fixed_alloc(T, n, 0.0, M_PI, 0.0, 0.0);

	F.function = &f;
	F.params = &m;

	gsl_integration_fixed(&F, &result, w);

	printf ("-------------使用Fixed point quadratures計算積分-------------\n\n");
	printf ("積分類型：gsl_integration_fixed_legendre\n\n");
	printf ("積分類固定點數：%d\n\n",n);
	printf ("計算結果 = % .8f\n\n\n", result);
	printf ("-------------使用高斯-勒讓德計算積分-------------\n\n");
	printf ("高斯點數：%d\n\n",n);
	printf ("計算結果：%.8f\n\n",gsl_gl(n));
	
	gsl_integration_fixed_free (w);
	return 0;
}
           

參考連結：

[1] https://www.gnu.org/s/gsl/manual/gsl-ref.pdf

[2] https://zh.wikipedia.org/wiki/高斯求積