数值积分之复化求积法

实际问题中，如果遇上积分的求解，是不是就纠结了呢？呵呵，我们知道积分值表现为x=a与x=b之间，y=0与y=F(x)所组成的曲边梯形面积。根据积分中值定理，我们知道，在[a,b]上可以找到一个函数值F(ξ)使得S＝F(ξ)*(b-a)。按照这种理解，人们不断去搜寻如何逼近这个F(ξ)，于是乎我们有了梯形公式，有了中矩形公式，辛甫生公式……但它们的效果并不是很理想，所以我们得另外寻找更加可靠的方法。

从一般情况出发，对于[a,b]区间内，通过寻找若干个节点用加权平均的方法，使得S＝∑λF(Xi)。复化求积就是通过对[a,b]进行n等分，根据精度需要选择相关的系数进行加权求积。这里的系数选择涉及到代数精度与Newton-Cotes公式里的柯特斯系数求解的问题，限于篇幅，此处就不展开讨论了。

下面给出常用的代数精度所用的加权系数：

一阶系数，1/2，1/2；

二阶系数，1/6，2/3，1/6；

四阶系数，7/90，16/45，2/15，16/45，7/90。

说明：n阶的系数至少有n阶的精度，但我们发现二阶与三阶的精度相当，四阶与五阶的精度相当；由于高阶公式的不稳定，故不宜采用。

一阶系数产生的是复化梯形公式，Tn=∑(F(Xk)/2+F(Xk+1)/2)*h。

二阶系数产生的是复化辛甫生公式，Sn=∑(F(Xk)/6+4*F(Xk+1/2)/6+F(Xk+1)/6)*h。

四阶系数产生的是复化柯特斯公式，Cn=∑(7*F(Xk)/90+16*F(Xk+1/4)/45+2*F(Xk+1/2)/15+16*F(Xk+3/4)/45+7*F(Xk+1)/90)*h。

既然已经有了公式，那就开始编程呗：

double ComTra(double a,double b,int N)//复化梯形
{
	double f=F(a)+F(b),x=a,h=(b-a)/N;
	for(int i=0;i<N-1;i++)
	{
		x+=h;
		f+=2*F(x);
	}
	return f*h/2;
}
double ComSim(double a,double b,int N)//复化辛甫生
{
	double f=0,x=a,h=(b-a)/N;
	for(int i=0;i<N;i++)
	{
		f+=F(x);
		x+=h/2;
		f+=4*F(x);
		x+=h/2;
		f+=F(x);
	}
	return f*h/6;
}
double ComNewCot(double a,double b,int N)//复化柯特斯
{
	double f=0,x=a,h=(b-a)/N;
	for(int i=0;i<N;i++)
	{
		f+=7*F(x);
		x+=h/4;
		f+=32*F(x);
		x+=h/4;
		f+=12*F(x);
		x+=h/4;
		f+=32*F(x);
		x+=h/4;
		f+=7*F(x);
	}
	return f*h/90;
}

在VS2010下进行运行，这里我试验的是F(x)=sin(x)/x的情况，积分区间是[0,1]，其中我令F(0)=1，以防出错。