GSL科学计算库——计算高斯-勒让德积分

假设计算下列积分：
$\int_0^\pi e^xcos(x)dx$
根据高斯勒让德积分，先将积分区间变换为 $[-1,1]$ ，使用线性变换：
$x=\frac{\pi-0}{2}t+\frac{\pi+0}{2}$ ，将原积分变为
$\int_0^\pi e^xcos(x)dx=-\frac{\pi}{2}e^{\frac{\pi}{2}t}\int_{-1}^1 e^{\frac{\pi}{2}t}cos(\frac{\pi}{2}t)dx$
最后再根据 $n$ 点求积系数和求积节点计算高斯-勒让德积分，高斯点和求积系数的选区如下图：
GSL科学计算库——计算高斯-勒让德积分

在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 $\sum_{i=1}^{n}\omega_if(x_i)$
参数解释：
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;
}

GSL科学计算库——计算高斯-勒让德积分

参考链接：
[1] https://www.gnu.org/s/gsl/manual/gsl-ref.pdf
[2] https://zh.wikipedia.org/wiki/高斯求积

秒客网

GSL科学计算库——计算高斯-勒让德积分

相关文章