来自「一个用C++编写的积分算法,可以用VC编译.」· 代码 · 共 444 行

#include  "stdlib.h"
#include  "math.h"
#include    "iostream.h"
#include  "stdio.h"
#include  "积分.h"
double gauss::derivate (double A[max],double B[max],double x,int n)
{

	int j;
		for (j=0;j<n;j++)
		{
			if(j==0)
				B[j]=A[j];
			else 
			{
				B[j]=A[j]+x*B[j-1];
			}
		}
		return B[n-1];
}
gauss::rowgauss (double A[max][max],double B[max],int n)
{
	int i,j,k;
	double temp;
	int row,line[max];
	double result[max];
	for(i=0;i<n;i++)
		line[i]=i;
	for(i=0;i<n;i++)
	{
		row=0;
		temp=0;
		for(j=0;j<n;j++)
		{
			if(temp<fabs(A[i][j]))
			{
				temp=fabs(A[i][j]);
				row=j;
			}
			else 
				continue;
		}
		if(temp==0)
			exit(0);
		else
		{
			for(k=0;k<n;k++)
			{
				temp=A[k][row];
				A[k][row]=A[k][i];
				A[k][i]=temp;
			}
			k=line[row];
			line[row]=line[i];
			line[i]=k;
			for(j=i+1;j<n;j++)
			{
				temp=A[j][i]/A[i][i];
				for(k=i;k<n;k++)
					A[j][k]=A[j][k]-temp*A[i][k];
				B[j]=B[j]-temp*B[i];
			}
		}
	}
	result[n-1]=B[n-1]/A[n-1][n-1];
	for(i=n-2;i>=0;i--)
	{
		temp=0;
		for(j=i+1;j<n;j++)
			temp+=A[i][j]*result[j];
		result[i]=(B[i]-temp)/A[i][i];
	}
	for(i=0;i<n;i++)
	{
		row=line[i];
		B[row]=result[i];
	}

}
gauss::newton (double A[max],double result[max],int n)
{
	double B[max],C[max];
	double temp,factor,r;
	int i,j;
	for(i=0;i<n-1;i++)
	{
		r=3;
		do
		{
			factor=1;
			j=10;
			temp=derivate(A,B,r,n-i);
			temp=temp/derivate(B,C,r,n-i-1);
			do
			{
				if(j==0)
					break;
				else
				{
					result[i]=r-factor*temp;
					j--;
					factor=factor/2;
				}
			}
			while(fabs(derivate(A,B,result[i],n-i))>fabs(derivate(A,B,r,n-i)));
			r=result[i];			
		}
		while (fabs(temp)>1e-5);
		derivate(A,B,result[i],n-i);

		
		for(j=0;j<n-i;j++)
		{
			A[j]=B[j];
		}
		
	}
	for(i=0;i<n-1;i++)
	{
		if(fabs(result[i])<1e-3)
			result[i]=0;
		else
			continue;
	}

}

gauss::legendre (double A[max],int n)
{
	int i,j;
	double B[max][max];
	B[0][0]=1;
	B[1][0]=0;
	B[1][1]=1;
	for(i=1;i<n;i++)
	{
		B[i+1][i+1]=(2*i+1)*B[i][i]/(i+1);
		B[i+1][i]=(2*i+1)*B[i][i-1]/(i+1);
		B[i+1][0]=-i*B[i-1][0]/(i+1);
		for(j=1;j<i;j++)
			B[i+1][j]=((2*i+1)*B[i][j-1]-i*B[i-1][j])/(i+1);
	}
	for(i=0;i<n;i++)
	{	A[i]=B[n-1][n-1-i];
//		cout<<A[i]<<"       ";
	}
	cout<<endl;
}
gauss::qiebixuefu (double A[max],int n)
{
	double B[max][max];
	int i,j;
	B[0][0]=1;
	B[1][0]=0;
	B[1][1]=1;
	for(i=1;i<n;i++)
	{
		B[i+1][i+1]=2*B[i][i];
		B[i+1][i]=2*B[i][i-1];
		B[i+1][0]=-B[i-1][0];
		for(j=1;j<i;j++)
			B[i+1][j]=2*B[i][j-1]-B[i-1][j];
	}
	for(i=0;i<n;i++)
	{	A[i]=B[n-1][n-1-i];
		cout<<A[i]<<"       ";
	}
	cout<<endl;

}
double   newton::coefficient(int k,int n)
{
	int i,j;
	double temp,temp1;
	double A[max];
	double B[max][max];
	temp=0;
	for(i=0;i<=n;i++)
		for(j=0;j<n;j++)
			B[i][j]=0;
	if(k!=0)
	{
		for(i=0;i<=n;i++)
			B[i][i]=1;
		for(i=1;i<=n;i++)
		{
			if(i<k)
			{
				for (j=0;j<i;j++)
				{
					if(j==0)
						B[i][0]=-i*B[i-1][0];
					else
						B[i][j]=B[i-1][j-1]-i*B[i-1][j];
				}
			}
			else if(i>k)
			{
				for (j=0;j<i;j++)
				{
					if(j==0)
						B[i-1][0]=-i*B[i-2][0];
					else
						B[i-1][j]=B[i-2][j-1]-i*B[i-2][j];
				}
			}
		}
		for(i=0;i<n;i++)
		{
			A[i]=B[n-1][i];
		//	cout<<A[i]<<"   ";
		}
		for(i=0;i<n;i++)
		{
			temp1=1;
			for(j=0;j<i+2;j++)
				temp1*=n;
			temp1=A[i]*temp1/(i+2);
			temp+=temp1;
		}		
	}
	else
	{
		B[0][0]=-1;
		for(i=0;i<=n;i++)
			B[i][i+1]=1;
		for(i=1;i<=n;i++)
			{	
				B[i][0]=-(i+1)*B[i-1][0];
				for (j=1;j<=i;j++)
						B[i][j]=B[i-1][j-1]-(i+1)*B[i-1][j];
			}
		for(i=0;i<=n;i++)
		{
			A[i]=B[n-1][i];
		//	cout<<A[i]<<"   ";
		}
		for(i=0;i<=n;i++)
		{
			temp1=1;
			for(j=0;j<i+1;j++)
				temp1*=n;
			temp1=A[i]*temp1/(i+1);
			temp+=temp1;
		}
	}

//	cout<<temp<<endl;
	return temp;

}
double newton::newtoncores (double a,double b,double (*f)(double),int n)
{
	double C[max];
	int i,j;
	double temp1,temp2,temp;
	int sign;
	for(i=0;i<=n;i++)
	{
		temp1=temp2=1;
		for(j=1;j<=i;j++)
			temp1*=j;
		for(j=1;j<=n-i;j++)
			temp2*=j;
		temp=coefficient(i,n);
		sign=(n-i)%2;
		if(sign)
			C[i]=-temp/(temp1*temp2*n);
		else
			C[i]=temp/(temp1*temp2*n);
	}
	temp=b-a;
	temp1=temp/n;
	temp2=0;
	for(i=0;i<=n;i++)
		temp2=temp2+C[i]*f(a+i*temp1);
	temp2=temp2*temp;
	return temp2;
}
double newton::newtoncores (double a,double b,double (*f)(double),int m,int n)
{
	double C[max];
	int i,j;
	double temp1,temp2,temp,result,r;
	int sign;
	for(i=0;i<=n;i++)
	{
		temp1=temp2=1;
		for(j=1;j<=i;j++)
			temp1*=j;
		for(j=1;j<=n-i;j++)
			temp2*=j;
		temp=coefficient(i,n);
		sign=(n-i)%2;
		if(sign)
			C[i]=-temp/(temp1*temp2*n);
		else
			C[i]=temp/(temp1*temp2*n);
	}
	result=0;
	temp=(b-a)/m;
	temp1=a;
	r=(b-a)/(m*n);
	for(i=0;i<m;i++)
	{
		temp2=0;
		for(j=0;j<=n;j++)
			temp2+=C[j]*f(temp1+j*r);
		temp2=temp2*temp;
		temp1=temp1+temp;
		result+=temp2;
	}
	return result;
}
double newton::romberg (double a,double b,double (*f)(double))
{
	double temp,r,m;
	double limit;
	double T[max][max];
	int k,i;
	k=1;
	r=b-a;
	T[0][0]=r*(f(a)+f(b))/2;
	do
	{
		m=1;
		for(i=1;i<k;i++)
			m*=2;
		temp=0;
		for(i=0;i<m;i++)
			temp+=f(a+(i+0.5)*r);
		T[k][0]=T[k-1][0]/2+r*temp/2;
		temp=4;
		for(i=1;i<=k;i++)
		{
			T[k][i]=(temp*T[k][i-1]-T[k-1][i-1])/(temp-1);
			temp*=4;
		}

		r=r/2;
	
		limit=fabs(T[k][k]-T[k-1][k-1]);
		k++;
	}
	while(limit>1e-3);

	return  T[k-1][k-1];	
}
double gauss::integral(double a,double b,double (*f)(double),int n)
{
	double value;
	int i,j;
	double temp1,temp2;
	double A[max];
	double result[max];
	double B[max];
	double C[max][max];
	legendre(A,n+1);
	newton(A,result,n+1);
	temp1=temp2=1;
	for(i=0;i<n;i++)
	{
		temp1*=-1;
		temp2*=1;
		B[i]=(temp2-temp1)/(i+1);			
	}
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
		{
			if(i==0)
				C[i][j]=1;
			else
				C[i][j]=C[i-1][j]*result[j];
		}
	}
	rowgauss(C,B,n);
	value=0;
	for(i=0;i<n;i++)
	{
		result[i]=result[i]*(b-a)/2+(a+b)/2;
		value+=B[i]*f(result[i]);
	}
	return value;
}
double f(double x)
{
	
	return(1/x);
}
void  main()
{
	double a,b;
	double result;
	double f(double);
	int n;
	a=1;
	b=3;
	n=5;
	gauss  aa;
	result=aa.integral (a,b,f,n);
	cout<<result<<endl;
/*	double f(double);
	int m,n;
	double a,b;
	n=2;
	m=4;
	a=1;b=3;
	double result;
	newton aa;
	result=aa.romberg(a,b,f);
	cout<<result<<endl;*/
/*	double A[max];
	int n;
	n=7;
	gauss  aa;
	aa.qiebixuefu (A,n);
//	aa.legendre (A,n);*/
/*	int n,i,j;
	double temp1,temp2,temp;
	double C[max][max];
	int sign;
	for(n=1;n<5;n++)
	{
		for(i=0;i<=n;i++)
		{
			temp1=temp2=1;
			for(j=1;j<=i;j++)
				temp1*=j;
			for(j=1;j<=n-i;j++)
				temp2*=j;
			temp=coefficient(i,n);
			sign=(n-i)%2;
			if(sign)
				C[n][i]=-temp/(temp1*temp2*n);
			else
				C[n][i]=temp/(temp1*temp2*n);
		}
		for(i=0;i<=n;i++)
			cout<<C[n][i]<<"    ";
		cout<<endl;
	}*/
}