初探C++连续无奇点函数的积分算法

定积分算法概述

本文主要介绍一下积分算法。定积分可以理解为函数与x轴、指定的左右界围成的有符号的面积代数和。那么最简单的算法求a到b上，f(x)函数的积分，就是把[a,b]区间平均分割为无数多个divx模块，每个模块都看成长度为f(a+divx)，宽度为divx的矩形，相乘后的有符号面积进行代数和后求出的即为定积分值，这种算法这里就不予实现了，网络上很多。今天将两种算法，起始还是属于积分中比较浅显的，给需要的同胞们解决一下燃眉之急。

自定义算法

我介绍的这种算法比较简单，是用梯形分割法进行反复迭代计算的。话不多说，po一下代码：

template<class C>
void my_Integral_trapzd(C& f1, const DB begin, const DB end, const int  num, const DB eps,DB& result)
/*************************************************
Function:		 my_Integral_trapzd
Description:	采用梯形迭代法的基本一重积分
Input:
*f1:			需要求积分的原始表达式
*begin:			定积分起点
*end:			定积分终点
*multiple:		对称分割的次数，为2^(num-1)次
*result:		存放积分结果
Output:
Return:
*ret:			void
*************************************************/
{
	DB x, tnm, sum, del;
	const DB EPS{ eps };
	static DB s;
	DB olds = 0.0;
	int it, j;
	for (int multiple = 0; multiple < num; multiple++)
	{
		if (multiple == 1)//仅分割轴向区间1次
		{
			s = 0.5*(end - begin)*(f1(begin) + f1(end));//赋值静态变量，准备迭代，并且返回值
			result = s;
		}
		else
		{
			for (it = 1, j = 1; j < multiple - 1; j++)  it <<= 1;//快速分裂分割点，1，2，4，8；
			tnm = it;//令tnm = 分裂迭代次数；
			del = (end - begin) / tnm;//所增加点的间隔；
			x = begin + 0.5*del;
			for (sum = 0.0, j = 0; j < it; j++, x += del)
				sum += f1(x);
			s = 0.5*(s + (end - begin)*sum / tnm);
		}
		if (multiple>5)//防止早期的假性收敛
		{
			if (fabs(s - olds) < EPS*fabs(olds) || (s == 0.0&&olds == 0.0))
			{
				result = s;
				return;
			}
			else  olds = s;
		}
	}
	result = s;
	return;
}
           

其中形参C&c，可以是被积函数对象，该函数对象是类中重载()运算符后直接使用的，重载方式如下：

class M_F
{
public:
	DB operator()(DB x)
	{
		DB Se = pow((1 + pow(abs(paras_alpha * x), paras_n)), 1.0 / paras_n - 1.0);//公式9
		DB temp = 1.0;
		temp *= paras_ks;
		temp *= pow(Se, paras_cauchy);
		temp *= pow(1.0 - pow(1.0 - pow(Se, paras_n / (paras_n - 1)), 1.0 - 1.0 / paras_n), 2);
		return temp;
	}
};
           

至此，算法需要的成分就编写完毕，接下来直接在main函数饮用就可以了。但是因为是自己编写的，没有加什么内存管理方式或者算法优化，所以导致算法运行时间较长。下面介绍一下GSL库的gsl_integration_qag算法。

GSL库gsl_integration_qag算法

对于该库，官方说明文档的解释是：

此函数适应性地施加积分规则，直到积分的估计F超过（a，b）在所需绝对和相对误差范围内实现的，epsabs和 epsrel。该函数返回最终近似值 result和绝对误差的估计值abserr。积分规则由的值确定key，应从以下符号名称中选择：

GSL_INTEG_GAUSS15（键= 1 ）；

GSL_INTEG_GAUSS21（键= 2 ）；

GSL_INTEG_GAUSS31（键= 3 ）；

GSL_INTEG_GAUSS41（键= 4 ）；

GSL_INTEG_GAUSS51（键= 5 ）；

GSL_INTEG_GAUSS61（键= 6 ）；

对应于15、21、31、41、51和61点高斯-克朗罗德规则。高阶规则可为平滑函数提供更好的准确性，而低阶规则可在函数包含局部困难（例如不连续性）时节省时间。

在每次迭代中，自适应积分策略均将具有最大误差估计值的区间一分为二。子间隔及其结果存储在所提供的内存中workspace。子间隔的最大数量由给出limit，不能超过工作空间的分配大小。

话不多说，下面给出该算法的调用方法：

//定义函数系数结构体
struct my_f_params { double alpha; double n; double ks; double cauchy; };
//定义被积原函数
double my_f(double x, void * p)
{
	struct my_f_params * params = (struct my_f_params *)p;
	double paras_alpha = params->alpha;
	double paras_n = params->n;
	double paras_ks = params->ks;
	double paras_cauchy = params->cauchy;

	DB Se = pow((1 + pow(abs(paras_alpha * x), paras_n)), 1.0 / paras_n - 1.0);//公式9
	DB temp = 1.0;
	temp *= paras_ks;
	temp *= pow(Se, paras_cauchy);
	temp *= pow(1.0 - pow(1.0 - pow(Se, paras_n / (paras_n - 1)), 1.0 - 1.0 / paras_n), 2);
	return temp;
}

int main()
{
	//定义变量
	double result, error;
	gsl_function F;
	my_f_params params = { paras_alpha, paras_n, paras_ks,paras_cauchy }; 
	gsl_integration_workspace * w= gsl_integration_workspace_alloc(1000);
	F.function = &my_f;
	F.params = &params;
	tgsl1 = GetTickCount();//获取GSL库计算的起始时间节点
	for (int i = 0; i < up.size(); i++)
	{
		gsl_integration_qag(&F,  down, up[i],0,1e-7, 1000, GSL_INTEG_GAUSS61,w, &result,&error);
		rst_gsl.push_back(result);
		//ntrail_gsl.push_back(w->size);
	}
	tgsl2 = GetTickCount();//获取GSL库计算的终止时间节点
	gsl_integration_workspace_free(w);
}
           

使用该方法，必须要调用GSL编译好的动态链接库，因为dll文件是看不到内部函数结构的，所以，这里没办法给大家po出来。本来想把该方法整理成执勤啊自定义函数的模式，这样更加符合C++的语言风格和特性，但是因为GSL库已经定义好gsl_function 格式，我试过static_cast强制转换，也不能将class类成员方法转变过去，网上说可以把需要调用的成员函数转变为static静态函数，但是operator属于符号重载方法，必须是成员函数，而且，符号重载方法不能设置为静态函数，所以幻想破灭，老老实实的用原装方法。然后我测试了一下100个数的积分时间，用的是#include<windows.h>库中的DWORD类。

下图是两种算法的比较时间：

好了，这篇到此结束，如果有更好用的库和算法，也欢迎大家留言。