main.cpp

此源码为高斯-拉格朗日积分法

#include <iostream>
#include <iomanip>
#include <fstream>
#include <math.h>

#define  precisions   8			//定义有效数字
using namespace std;
typedef double num_type;

const num_type epciron=1e-10;
const num_type x_0=-1;
const num_type x_n=1;

int Gauss_2(double* *mat_A,double* mat_B,double* xx,int n);
//***************************************************************************
//*复化梯形方法
//***************************************************************************
int Re_Trapezoidal()
{
	int i,j,n=1;
	num_type error;
	num_type h;
	ofstream out_stream;
	out_stream.open("result1.dat");
	printf("复化梯形方法!\n");
	cout<<setprecision(precisions);
	out_stream<<setprecision(precisions);
	while (1)
	{
		n=n*2;
		num_type *x,*g,*u;	
		x=new num_type[n+1];	//创造动态一维数组存储节点xi的值
		g=new num_type[n+1];	//创造动态一维数组存储方程常数项的值g(xi)
		u=new num_type[n+1];	//创造动态一维数组存储方程的解u(xi)
		num_type* *A=new num_type*[n+1];//创造动态二维数组存储系数矩阵
		for (i=0;i<n+1;i++)
		A[i]=new num_type[n+1];

		h=(x_n-x_0)/n;
		for (i=0;i<n+1;i++)
			x[i]=x_0+i*h;
		for (i=0;i<n+1;i++)
		{
			for (j=0;j<n+1;j++)
			{
				if (j==0)
					A[i][j]=h*exp(x[0]*x[i])/2;
				else if (j==n)
					A[i][j]=h*exp(x[n]*x[i])/2;
				else
					A[i][j]=h*exp(x[i]*x[j]);
				if (j==i)
					A[i][j]+=1;
			}
			g[i]=exp(4*x[i])+(exp(x[i]+4)-exp(-x[i]-4))/(x[i]+4);
		}

		Gauss_2(A,g,u,n+1);
		error=0;
		for (i=0;i<n+1;i++)
		{
			if((i==0)||(i==n)) error+=(h/2)*(exp(4*x[i])-u[i])*(exp(4*x[i])-u[i]);
			else error+=h*(exp(4*x[i])-u[i])*(exp(4*x[i])-u[i]);
		}
		cout<<"n+1= "<<n+1<<"error="<<' '<<error<<endl;
		if(error<=epciron)
		{			
		/*	for (i=0;i<n+1;i++)
				{
				cout<<"xxx"<<x[i]<<' ';
				cout<<"TURE"<<exp(4*x[i])<<' ';
				cout<<"SIMI"<<u[i]<<' ';
				cout<<endl;
				}
			cout<<endl;
			for (i=0;i<n+1;i++)
				cout<<u[i]<<' ';
			cout<<endl;*/
			for (i=0;i<n+1;i++)//将细化结果保存到result1.dat中
				{
				out_stream<<x[i]<<' '<<u[i]<<' ';
				}
			out_stream<<endl;
			out_stream.close();
			break;
		}
		delete []x;
		delete []g;
		delete []u;
		for (i=0;i<n+1;i++)
			delete [] A[i];
		delete [] A;		
	}

	

	printf("Hello World!\n");
	return 1;
	
}
//***************************************************************************
//*复化simpson方法
//***************************************************************************
int Re_Simpson()
{
	int i,j,n=1;
	num_type error;
	num_type h;
	ofstream out_stream;
	out_stream.open("result2.dat");
	printf("复化simpson方法!\n");
	cout<<setprecision(precisions);
	out_stream<<setprecision(precisions);
	while (1)
	{
		n=n*2;
		num_type *x,*g,*u;	
		x=new num_type[n+1];	//创造动态一维数组存储节点xi的值
		g=new num_type[n+1];	//创造动态一维数组存储方程常数项的值g(xi)
		u=new num_type[n+1];	//创造动态一维数组存储方程的解u(xi)
		num_type* *A=new num_type*[n+1];//创造动态二维数组存储系数矩阵
		for (i=0;i<n+1;i++)
		A[i]=new num_type[n+1];

		h=(x_n-x_0)/n;
		for (i=0;i<n+1;i++)
			x[i]=x_0+i*h;
		for (i=0;i<n+1;i++)
		{
			for (j=0;j<n+1;j++)
			{
				if (j==0)
					A[i][j]=h*exp(x[0]*x[i])/3;
				else if (j==n)
					A[i][j]=h*exp(x[n]*x[i])/3;
				else if (j%2)
					A[i][j]=h*exp(x[i]*x[j])*4/3;
				else 
					A[i][j]=h*exp(x[i]*x[j])*2/3;
				if (j==i)
					A[i][j]+=1;
			}
			g[i]=exp(4*x[i])+(exp(x[i]+4)-exp(-x[i]-4))/(x[i]+4);
		}
        
		Gauss_2(A,g,u,n+1);
		
	//***************************************************************************
		error=0;
		for (i=0;i<n+1;i++)
		{
			if((i==0)||(i==n)) error+=h/3*(exp(4*x[i])-u[i])*(exp(4*x[i])-u[i]);
			else if (i%2) error+=4*h/3*(exp(4*x[i])-u[i])*(exp(4*x[i])-u[i]);
			else error+=2*h/3*(exp(4*x[i])-u[i])*(exp(4*x[i])-u[i]);
		}
		cout<<"n+1= "<<n+1<<"error="<<' '<<error<<endl;
		if(error<=epciron)
		{			
			for (i=0;i<n+1;i++)
			{
				cout<<"xxx"<<x[i]<<' ';
				cout<<"TURE"<<exp(4*x[i])<<' ';
				cout<<"SIMI"<<u[i]<<' ';
				cout<<endl;
			}
			cout<<endl;
			for (i=0;i<n+1;i++)
				cout<<u[i]<<' ';
			cout<<endl;
			for (i=0;i<n+1;i++)//将细化结果保存到result2.dat中
				{
				out_stream<<x[i]<<' '<<u[i]<<' ';
				}
			out_stream<<endl;
			out_stream.close();
			break;
		}
		delete []x;
		delete []g;
		delete []u;
		for (i=0;i<n+1;i++)
			delete [] A[i];
		delete [] A;		
	}
	printf("OK!\n");
	return 1;
	
}

//***************************************************************************
//*Gauss_Legendre方法
//***************************************************************************
int Gauss_Legendre()
{
	int i,j,n=0;
	num_type error;
	ifstream in1_stream;
	ifstream in2_stream;
	ofstream out_stream;

	in1_stream.open("table1.dat");
	in2_stream.open("table2.dat");
	out_stream.open("result3.dat");
	printf("Gauss_Legendre方法!\n");
	cout<<setprecision(precisions);
	out_stream<<setprecision(precisions);
	while (1)
	{
		n=n+1;
		if (n==10) 
		{cout<<endl<<"can not do it!!"<<endl;return 0;}
		num_type *x,*g,*u,*AA;	
		x=new num_type[n];	//创造动态一维数组存储节点xi的值
		AA=new num_type[n];	//创造动态一维数组存储Ai的值
		g=new num_type[n];	//创造动态一维数组存储方程常数项的值g(xi)
		u=new num_type[n];	//创造动态一维数组存储方程的解u(xi)
		num_type* *A=new num_type*[n];//创造动态二维数组存储系数矩阵
		for (i=0;i<n;i++)
		A[i]=new num_type[n];
		for (i=0;i<n;i++)
		{
		 in1_stream>>x[i];
		 in2_stream>>AA[i];
		}

		for (i=0;i<n;i++)
		{
			for (j=0;j<n;j++)
			{
				A[i][j]=AA[j]*exp(x[i]*x[j]);
				if(j==i)
				A[i][j]+=1;
			}
			g[i]=exp(4*x[i])+(exp(x[i]+4)-exp(-x[i]-4))/(x[i]+4);
		}
		Gauss_2(A,g,u,n);
		
	//***************************************************************************
		error=0;
		for (i=0;i<n;i++)
		{
			error+=AA[i]*(exp(4*x[i])-u[i])*(exp(4*x[i])-u[i]);
		}
		cout<<"n= "<<n<<"error="<<' '<<error<<endl;
		if(error<=epciron)
		{			
			for (i=0;i<n;i++)
			{
			cout<<"xxx"<<x[i]<<' ';
			cout<<"TURE"<<exp(4*x[i])<<' ';
			cout<<"SIMI"<<u[i]<<' ';
			cout<<endl;
			}
			cout<<endl;
			for (i=0;i<n;i++)
				cout<<u[i]<<' ';
			cout<<endl;
			for (i=0;i<n;i++)//将细化结果保存到result3.dat中
			{
			out_stream<<x[i]<<' '<<u[i]<<' ';
			}
			out_stream<<endl;
			in1_stream.close();
			in2_stream.close();
			out_stream.close();
			break;
		}
		delete []x;
		delete []g;
		delete []u;
		for (i=0;i<n;i++)
			delete [] A[i];
		delete [] A;		
	}
	


	printf("ok!\n");
	return 1;
	
}
//***************************************************************************
//*列主元素Gauss消去法求解线性方程组，输入A，B矩阵
//返回值0：方法错误；返回值1：正确
//***************************************************************************
int Gauss_2(double* *mat_A,double* mat_B,double* xx,int n)
{
	int k,i,j;
	int find_k;
	double m,max_num,temp;

	
	for (k=0;k<n-1;k++)   
	{
		find_k=k;
		max_num=mat_A[k][k];
		for (i=k;i<n-1;i++)
		{
			if (max_num<mat_A[i+1][k])    //找第k列k行之后最大元素  
			{
				find_k=i+1;
                max_num=mat_A[i+1][k];
			}
		}
		if (find_k!=k)        //找第k列k行之后最大元素后，矩阵交换行
		{
			for (i=0;i<n;i++)
			{
				temp=mat_A[k][i];
				mat_A[k][i]=mat_A[find_k][i];
                mat_A[find_k][i]=temp;
			}
			temp=mat_B[k];
			mat_B[k]=mat_B[find_k];
            mat_B[find_k]=temp;
		}

        //化为上三角阵
		if (mat_A[k][k]==0)
		return 0;
		for (i=k+1;i<n;i++)
		{
			m=mat_A[i][k]/mat_A[k][k];  
			mat_A[i][k]=0;

			for (j=k+1;j<n;j++)
			{
			    mat_A[i][j]=mat_A[i][j]-m*mat_A[k][j];
			}
            mat_B[i]=mat_B[i]-m*mat_B[k];
		}
		
	}
	//回代过程
	xx[n-1]=mat_B[n-1]/mat_A[n-1][n-1];
	for (k=n-2;k>=0;k--)
	{
		m=0;
		for (j=k+1;j<n;j++)
		{
			m=m+mat_A[k][j]*xx[j];
		}
		xx[k]=(mat_B[k]-m)/mat_A[k][k];
	}
	return 1;
}
void main()
{
	int i;
	printf("Hello World!\n");

    num_type temp[8194];
	num_type temp2[258];
	ifstream in_stream;
	ifstream in2_stream;
	ofstream out_stream;
	ofstream out2_stream;
	ofstream out3_stream;
	in_stream.open("result1.dat");
	in2_stream.open("result1.dat");
	out_stream.open("New_result1.dat");
	out2_stream.open("ture_result.dat");
	out3_stream.open("打印.dat");
	cout<<setprecision(precisions);
	out_stream<<setprecision(precisions);
	out2_stream<<setprecision(precisions);
	out3_stream<<setprecision(precisions);

		for (i=0;i<8194;i++)
		{
		 in_stream>>temp[i];
		}
		for (i=0;i<258;i++)
		{
		 in2_stream>>temp2[i];
		}
		for (i=0;i<129;i++)//将细化结果保存到result3.dat中
			{

			out_stream<<temp[i*64]<<' '<<temp[i*64+1]<<' ';
			out2_stream<<temp[i*64]<<' '<<exp(4*temp[i*64])<<' ';
			out3_stream<<setw(14)<<temp[i*64]<<' '<<"【tu】"<<setw(14)<<exp(4*temp[i*64])<<' '<<"【s1】"<<setw(14)<<temp[i*64+1]<<' ';
			out3_stream<<"【s2】"<<setw(14)<<temp2[i*2+1]<<' ';
			 out3_stream<<endl;

			}
	out_stream<<endl;
	in_stream.close();
	in2_stream.close();
	out_stream.close();
	out2_stream.close();
	out3_stream.close();
	//Re_Trapezoidal();
    //Re_Simpson();
	//Gauss_Legendre();
	//printf("Hello World!\n");
}