积分算法.cpp

数值分析：数值积分算法

#include "math.h"
#include "stdio.h"
#define aa -1
#define bb 1
#define err 1e-10
#define N 200


//选主元gauss消去法，用来计算ui//
void functiongauss(double a[][N+1],double ui[N+1],double g[N+1],int n)    
{ 
	
	int i,j,k; 
	double temp,s,l; 
	
	for(k=0;k<=n;k++) 
	{ 
		i=k; 		                                                    //选列主元 
		for(j=k+1;j<=n;j++) 
		{ if(fabs(a[j][k])>fabs(a[i][k])) 
		i=j; 
		} 
		
		if(i!=k) 		                                               //换行 
		{	for(j=k;j<=n;j++) 
		{ 
			temp=a[k][j]; 
			a[k][j]=a[i][j]; 
			a[i][j]=temp; 
			
		} 
		temp=g[k];
		g[k]=g[i];
		g[i]=temp;
		}
		for(i=k+1;i<=n;i++) 			                          //消元 
		{ 
			l=1.0*a[i][k]/a[k][k]; 
			for(j=k+1;j<=n;j++) 
			{a[i][j]=a[i][j]-a[k][j]*l;}
			g[i]=g[i]-g[k]*l;
		} 
        
		
	} 
	ui[n]=g[n]/a[n][n]; 	                                         //回代 
	for(k=n-1;k>=0;k--) 
	{ 
		s=0.0; 
		for(j=k+1;j<=n;j++) 
		{ 
			
			s+=a[k][j]*ui[j]; 
			
		} 
		ui[k]=(g[k]-s)/a[k][k]; 
	} 
}

//复化梯形积分法//
void tixing()    
{   FILE *fp;
   	double h,sum;
	double xi[N+1]={0},ui[N+1]={0},uxi[N+1]={0},g[N+1]={0},lenda[N+1]={0};
    double a[N+1][N+1]={0};

	int i,j,n,k;
    fp=fopen("tixing.txt","w");
	n=1;
 	do 
 	{	
		sum=0;
	
 	    n++;
	    for (k=0;k<=n;k++)
		{
		  h=2.0/n;
		  xi[k]=aa+k*h;
		  g[k]=exp(4*xi[k])+(exp(4+xi[k])-exp(-xi[k]-4))/(xi[k]+4);
		  uxi[k]=exp(4*xi[k]);
		}
        for (k=0;k<=n;k++)
        {         
			for (j=0;j<=n;j++)
			{
			 if(j==0||j==n) a[k][j]=0.5*h*exp(xi[k]*xi[j]);
             else a[k][j]=h*exp(xi[k]*xi[j]);
			}
			if(k==0||k==n) 	lenda[k]=0.5*h;
		
			else lenda[k]=h;

		}
		for (k=0;k<=n;k++)
        {
      		for (j=0;j<=n;j++)
			{
             if(k==j) a[k][j]=1+a[k][j];
			}

        }

			
 	   functiongauss(a,ui,g,n);
       for (k=0;k<=n;k++)
	   {
	    sum+=lenda[k]*(uxi[k]-ui[k])*(uxi[k]-ui[k]);
	   }
  	 if(n>=350) break;
 	} while(sum>err);
	printf("%d  %.10e\n",n,sum);
	for (k=0;k<=n;k++)
	{
	   fprintf(fp,"%d  %f   %f   %f\n",k,xi[k],uxi[k],ui[k]);
	}
	fclose(fp);
}


//复化simpson积分法//
void simpson()
{
	FILE *fp;
   	double h,sum;
	double xi[N+1]={0},ui[N+1]={0},uxi[N+1]={0},g[N+1]={0},lenda[N+1]={0};
	double a[N+1][N+1]={0};
	int i,j,n,k;
    fp=fopen("simpson.txt","w");
	n=1;
 	do 
 	{	
       sum=0;
	
 	   n++;
	for (k=0;k<=n;k++)
	{
		h=2.0/n;
		xi[k]=aa+k*h;
		g[k]=exp(4*xi[k])+(exp(4+xi[k])-exp(-xi[k]-4))/(xi[k]+4);
		uxi[k]=exp(4*xi[k]);
	}
	for (k=0;k<=n;k++)
	{         
		for (j=1;j<n;j++)
		{
		    if(j%2==0) a[k][j]=(2.0/3.0)*h*exp(xi[k]*xi[j]);
			else a[k][j]=(4.0/3.0)*h*exp(xi[k]*xi[j]);
		}
	    a[k][0]=(h/3.0)*exp(xi[k]*xi[0]);
		a[k][n]=(h/3.0)*exp(xi[k]*xi[n]);

	}
	for (k=0;k<=n;k++)
	{
      		for (j=0;j<=n;j++)
			{
				if(k==j) a[k][j]=1+a[k][j];
			}
			
	}
	lenda[0]=h/3.0;
	lenda[n]=h/3.0;

	for (k=1;j<k;j++)
	{
		if(k%2==0) lenda[k]=(2.0/3.0)*h;
		else lenda[k]=(4.0/3.0)*h;
	}

	
	
	functiongauss(a,ui,g,n);
	for (k=0;k<=n;k++)
	   {
		sum+=lenda[k]*(uxi[k]-ui[k])*(uxi[k]-ui[k]);
	   }
// 	 if(n>=200) break;
	} while(sum>err);
	printf("%d  %.10e\n",n,sum);
	for (k=0;k<=n;k++)
	{
		fprintf(fp,"%d  %f   %f   %f\n",k,xi[k],uxi[k],ui[k]);
	}
	fclose(fp);
}


//gauss积分法//
void gauss()
{
	FILE *fp;
   	double h,sum;
	double ui[N+1]={0},uxi[N+1]={0},g[N+1]={0};
	double xi[N+1]={-0.9681602395,-0.8360311073,-0.6133714327,-0.3242534234,0,0.3242534234,0.6133714327,0.8360311073,0.9681602395};
	double Ai[9]={0.0812743884,0.1806481607,0.2606106964,0.3123470770,0.3302393550,0.3123470770,0.2606106964,0.1806481607,0.0812743884};
	double a[N+1][N+1]={0};
	int i,j,n,k;
    fp=fopen("gauss.txt","w");
	n=8;
	sum=0;
	
	for (k=0;k<=n;k++)
	{
		g[k]=exp(4*xi[k])+(exp(4+xi[k])-exp(-xi[k]-4))/(xi[k]+4);
		uxi[k]=exp(4*xi[k]);
	}
	for (k=0;k<=n;k++)
	{         
		for (j=0;j<=n;j++)
		{
		  a[k][j]=Ai[j]*exp(xi[k]*xi[j]);
		}
		
	}
	for (k=0;k<=n;k++)
	{
      		for (j=0;j<=n;j++)
			{
				if(k==j) a[k][j]=1+a[k][j];
			}
			
	}
	
	
	functiongauss(a,ui,g,n);
	for (k=0;k<=n;k++)
	   {
		sum+=Ai[k]*(uxi[k]-ui[k])*(uxi[k]-ui[k]);
	   }
	

	printf("%d  %.10e\n",n,sum);
	for (k=0;k<=n;k++)
	{
		fprintf(fp,"%d  %.10f   %.10f   %.10f\n",k,xi[k],uxi[k],ui[k]);
	}
	fclose(fp);
}



//主函数//
main()
{	
	// 	double xi[N+1]={0},ui[N+1]={0},uxi[N+1]={0},g[N+1]={0};
	//     double a[N+1][N+1]={0};

	
	printf("which method to operate:\n");
	printf("please input: tixing: 1  simpson: 2  gauss: 3\n");
	 int x;
	scanf("%d",&x);
	switch(x)
	{
	case 1:    tixing();
 		break;
	case 2: 	simpson();
 		break;
	case 3:	   gauss();
		break;
	default:
		break;
	}
}