初探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類。

下圖是兩種算法的比較時間：

好了，這篇到此結束，如果有更好用的庫和算法，也歡迎大家留言。