integral.cpp

C++实现的数值分析源程序

//梁霄 11/24/2005
#include"Integral.h"
#define MAXSIZE 100
//被积函数
double CommonIntegral::Function_T(double x)//梯形法被积函数
{
	if(x==0)
		x=0.00001;
	return sin(x)/x;
}
double CommonIntegral::Function_S(double x)//辛普森法被积函数
{
	return 1+exp(-x)*sin(4*x);
}
double CommonIntegral::Function_NC(double x)//牛顿-柯特斯法被积函数
{
	return 1+exp(-x)*sin(4*x);
}
//一般机械求积
//三种求积公式：梯形法，牛顿-科特斯法，辛普森法
double CommonIntegral::Trapezoid(double lower,double upper)
{
	return (upper-lower)*(Function_T(lower)+Function_T(upper))/2;
}
double CommonIntegral::Simpson(double lower,double upper)
{
	return (upper-lower)*(Function_S(lower)
							+4*Function_S((upper+lower)/2)
							+Function_S(upper))/6;
}
double CommonIntegral::Newton_Cotes(double lower,double upper)
{
	int k;
	double step=(upper-lower)/4;
	double x[5];

	x[0]=0.000000001;
	for(k=1;k<=3;k++)
		x[k]=lower+step*k;
	x[4]=upper;

	return (upper-lower)*(7*Function_NC(x[0])
							+32*Function_NC(x[1])
							+12*Function_NC(x[2])
							+32*Function_NC(x[3])
							+7*Function_NC(x[4]))/90;
}
//复化机械求积
//复化梯形法，复化牛顿-科特斯法，复化辛普森法
double ComplexIntegral::ComplexTrapezoid(double lower,double upper,int segmentnum=1)
{
	int k;
	double step=(upper-lower)/segmentnum;
	double x[MAXSIZE];

	x[0]=0.000000001;
	for(k=1;k<=segmentnum-1;k++)
	{
		x[k]=lower+step*k;
	}
	x[segmentnum]=upper;

	double temp=0.0;
	for(k=1;k<=segmentnum-1;k++)
		temp+=2*Function_T(x[k]);

	return (upper-lower)*(Function_T(x[0])
							+temp
							+Function_T(x[segmentnum]))/(2*segmentnum);
}
double ComplexIntegral::ComplexSimpson(double lower,double upper,int segmentnum=1)
{
	int k;
	double step=(upper-lower)/segmentnum;
	double x[MAXSIZE];

	for(k=0;k<=segmentnum-1;k++)
	{
		x[k]=lower+step*k;
	}
	x[segmentnum]=upper;

	double y[MAXSIZE];
	for(k=0;k<=segmentnum-1;k++)
		y[k]=(x[k]+x[k+1])/2;

	double tempx=0.0;
	for(k=1;k<=segmentnum-1;k++)
		tempx+=2*Function_S(x[k]);

	double tempy=0.0;
	for(k=0;k<=segmentnum-1;k++)
		tempy+=4*Function_S(y[k]);

	return (upper-lower)*(Function_S(x[0])
							+tempx
							+tempy
							+Function_S(x[segmentnum]))/(6*segmentnum);
}
double ComplexIntegral::ComplexNewton_Cotes(double lower,double upper,int segmentnum=1)
{
	int k;
	double step=(upper-lower)/segmentnum;
	double x[MAXSIZE];

	for(k=0;k<=segmentnum-1;k++)
	{
		x[k]=lower+step*k;
	}
	x[segmentnum]=upper;

	double y[MAXSIZE];
	for(k=0;k<=segmentnum-1;k++)
		y[k]=(x[k]+x[k+1])/2;

	double z[MAXSIZE];
	for(k=0;k<=segmentnum-1;k++)
	{
		z[k*2]=(x[k]+y[k])/2;
		z[k*2+1]=(y[k]+x[k+1])/2;
	}

	double tempx=0.0;
	for(k=1;k<=segmentnum-1;k++)
		tempx+=14*Function_NC(x[k]);

	double tempy=0.0;
	for(k=0;k<=segmentnum-1;k++)
		tempy+=12*Function_NC(y[k]);

	double tempz=0.0;
	for(k=0;k<=2*segmentnum-1;k++)
		tempz+=32*Function_NC(z[k]);

	return (upper-lower)*(7*Function_NC(x[0])
							+tempx
							+tempy
							+tempz
							+7*Function_NC(x[segmentnum]))/(90*segmentnum);
}
//变步长梯形机械求积
double AltIncrementIntegral::AltIncrementTrapezoid(double lower,double upper,double e)
{
	int n=1;

	while(pow((ComplexTrapezoid(lower,upper,n)-ComplexTrapezoid(lower,upper,n*2)),2)>=e)
		n*=2;
	return ComplexTrapezoid(lower,upper,n*2);
}