gaosi.cpp

这是一个关于数值方法的求非线性方程的算法

/*****************************************************************************/
/*           解线性方程组的直接方法（Gauss、Gausspp、Doolittle）             */
/*****************************************************************************/
#include"stdio.h"
#include"math.h"
#define N 100
int n,flag;
float A[N][N],C[N][N],y[N],I=1,O=0;
/*****************************************************************************/

/*****************************************************************************/
float Pivoting(int k)     
{
	/* 选主元函数*/
	int i;
	float templ;   
    for(i=k;i<=n;i++)
		if(fabs(A[i][k])>templ)
		{
			templ=fabs(A[i][k]);    /*将系数矩阵第i列中从第i行到第n行中最大值作为新主元*/
			flag=i;                /* 标志主元templ所在行,作为行交换函数L_exchange()的形参*/
		}
	return templ;
}
/*****************************************************************************/

/*****************************************************************************/
void L_exchange(int r,int k)         
{
	/*行交换函数*/
	float B[N];
	int i;
	for(i=1;i<=n+1;i++)    /*将新（大）主元templ所在行的元素存入数组B[i]*/
		B[i]=A[r][i];
	for(i=1;i<=n+1;i++)   /* 将小主元所在行k与主元所在行r进行交换*/
		A[r][i]=A[k][i];
	for(i=1;i<=n+1;i++)   
		A[k][i]=B[i];
}
/*****************************************************************************/

/*****************************************************************************/
void print()               
{
	/*输出函数*/
	int i,j;
	for(i=1;i<=n;i++)
	{
		for(j=1;j<=n+1;j++)
			printf("%9.3f",A[i][j]);
		printf("\n");
	}
	printf("\n");
}
/*****************************************************************************/

/*****************************************************************************/
void Gresult()      
{
	/*回代求解函数*/
	int k,j;
    float x[N];
	float memery=0;
    x[n]=A[n][n+1]/A[n][n];         /* 回代求解过程 */
    for( k=n-1;k>=1;k--)           
    {
		for(j=k+1;j<=n;j++)
	    memery=memery+A[k][j]*x[j];
        x[k]=(A[k][n+1]-memery)/A[k][k];
	    memery=0;                /* 将计数器memery的值置零以便重新计数*/
	}
	printf("\nHere are the roots for this linear equations:\n");
	for(j=n;j>=1;j--)				/*输出方程组的根*/
		printf("\n\tx%d=%f\n",j,x[j]);
	printf("\n");
}
/*****************************************************************************/

/*****************************************************************************/
void Resultant_Gauss()         
{
	/*Gauss顺序消去法的消元函数*/
	int k,i,j,t;
	float templ;
    for(k=1;k<=n-1;k++)                      
	{
		for(t=k;t<=n;t++)      /*检验方程组系数矩阵的对角元是否为0*/
			if(A[t][k]==0)
				printf("\nError!!The coefficient matrix A is singular!\n");/*输出方程组系数矩阵奇异的信息*/
		for(i=k+1;i<=n;i++)             /* 消元过程 */
		{
			templ=A[i][k]/A[k][k];     
	        for(j=k;j<=n+1;j++)
				A[i][j]=A[i][j]-A[k][j]*templ;
		}
	} 
}
/*****************************************************************************/

/*****************************************************************************/
void Resultant_GaussPP()         
{
	int k,i,j;
	float ark;
	float templ;
	for(k=1;k<=n-1;k++)                      
	{
		ark=Pivoting(k);
        if(ark==0)     /*判断方程是否为线性方程，即是否合法*/
			printf("\nError!!The coefficient matrix A is singular!"); /*输出方程组系数矩阵奇异的信息*/  
	    else if(flag!=k)
			L_exchange(flag,k);      
	    for(i=k+1;i<=n;i++)             /*消元过程*/
		{
			templ=A[i][k]/A[k][k];
            for(j=k;j<=n+1;j++)
				A[i][j]=A[i][j]-A[k][j]*templ;
		}
	} 
}
/*****************************************************************************/

/*****************************************************************************/
void H_Doolittle() 
{
	/* Doolittle分解法*/
	int i,r,k;
	float temp=0,templ=0;
	if(A[1][1]==0)
		printf("\nError!!The coefficient matrix is singular!\n");/*输出方程组系数矩阵奇异的信息*/
	else
	{
		for(i=2;i<=n;i++)
			A[i][1]=A[i][1]/A[1][1];
		r=2;
		while(r<=n)
		{
			for(i=r;i<=n;i++)
			{
				for(k=1;k<r;k++)
					templ=templ+(A[r][k]*A[k][i]);      
			    A[r][i]=A[r][i]-templ;
                templ=0;
			}
			for(i=r+1;i<=n;i++)
			{
				for(k=1;k<r;k++)
					templ=templ+(A[i][k]*A[k][r]);            
			    A[i][r]=(A[i][r]-templ)/A[r][r];
                templ=0;
			}
			r++;
		}
	}   
}
/*****************************************************************************/

/*****************************************************************************/

void Doolittle_L() 
{
	/*Doottle下三角分解及结果输出*/
	float x[N];  /*定义数组x[N],用于存放方程组的解*/
    float templ=0;
    int i, r;
    x[n]=y[n]/A[n][n];
    printf("\nHere are the roots for this linear equations:\n");
	for(r=n-1;r>=1;r--)
    {
		for(i=r+1;i<=n;i++)
			templ+=A[r][i]*x[i];
		x[r]=(y[r]-templ)/A[r][r];
        templ=0;
	}
	for(i=1;i<=n;i++)
    {
		printf("\tx%d=%9.3f",i,x[i]);
		printf("\n");
	}
} 
/*****************************************************************************/

/*****************************************************************************/
void Doolittle_U()
{
	
	/*Doottle上三角阵U的求解函数*/
	float templ=0;
    int i,r;
    y[1]=A[1][n+1];
    for(r=2;r<=n;r++)
    {
		for(i=1;i<=r-1;i++)
        templ=templ+A[r][i]*y[i];
        y[r]=A[r][n+1]-templ;
        templ=0;
	}
	
}
/*****************************************************************************/

/*****************************************************************************/
void Print_D_L()
{
	/*Doottle分解法,下三角矩阵L的输出函数*/
	int i,j; 
	int x=1,y=1;
	printf("\n Here is the Doolittle_L:\n");
	for(i=1;i<=n;i++)
    {
		for(j=1;j<=n;j++)
		{
			if(i>j)
				printf("%9.3f",A[i][j]);
			else if(i==j)
				printf("%9.3f",I);
			else
				printf("%9.3f",O);
		}
		printf("\n");
	}
}
/*****************************************************************************/

/*****************************************************************************/
void Print_D_U()
{
	/*Doottle分解法,上三角矩阵U的输出函数*/
	int i,j;
	float f=1.0;
	printf("\nHere is the Doolittle_U:\n");
    for(i=1;i<=n;i++)
    {
		for(j=1;j<=n;j++)
		{
			if(i<=j)
			{printf("%9.3f",A[i][j]);
				f=f*A[i][j];	/*计算系数矩阵A的行列式值*/
			}
			else
				printf("%9.3f",O);
		}
		printf("\n");
	}
	if(f==0)								/* 判断系数矩阵A的行列式值是否为0 */
		printf("Error!!! The coefficient matrix is singular!");	
											/* 若是，则输出系数矩阵A奇异的信息 */
	else									/* 若不是，则输出系数矩阵A的行列式值 */								
	{
		printf("The value of determinant A is:");		/*打印输出系数矩阵A的行列式值*/
		printf("%9.3f",f);
	}
}
/*****************************************************************************/


/*****************************************************************************/
void Gauss()
{
	Resultant_Gauss();
	Gresult();
}
/*****************************************************************************/

/*****************************************************************************/
GaussPP()
{
	Resultant_GaussPP();
	Gresult();
}
/*****************************************************************************/

/*****************************************************************************/
void Doolittle()
{
	H_Doolittle();
	Print_D_L();
	Print_D_U();
	Doolittle_U();
	Doolittle_L();
}
/*****************************************************************************/

/*****************************************************************************/
void Choose()
{
	/*菜单函数，负责完成对解题方法的选择以及是否退出程序*/
	int s=0;
	while(s!=4)
	{
		printf("If you want to run the program to solve this linear equations,\n");
		printf("please make sure that 's' is from 1 to 3; if you want to break \n");
		printf("out of program,please input 4 !\n");
		printf("\n*******************************************\n");
		printf("***\t 1 to choose Gauss().\t\t***\n");
		printf("***\t 2 to choose Gausspp().\t\t***\n");
		printf("***\t 3 to choose Doottle().\t\t***\n");
		printf("***\t 4 to break programer!\t\t***\n");
		printf("*******************************************\n");

		printf("Please enter your choice: ");		/*提示用户选择解题方法*/
		scanf("%d",&s);
		switch(s)
		{
		case 1:            /*输入1,则选择高斯顺序消去法解线性方程组*/
			Gauss();
			break;
		case 2:             /*输入2,则选择高斯列主元素消去法解线性方程组*/
			GaussPP();
			break;
		case 3:            /*输入3,则选择Doottle分解法解线性方程组*/
			Doolittle();
			break;
		case 4:         /*若输入4,则退出程序*/
			break;
		}
	}
}
/**********************************************************************************/

/**********************************************************************************/
void Input()        /*输入函数*/
{
	int i,j;
    printf("\nPlease input dimension of the linear equations:\n");
	printf("* n=");
    scanf("%d",&n);				/*输入线性方程组的维数*/
   	printf("\nPlease input augmented matrix (A,b) of the linear equations:\n");
	printf("(A,b)=\n\n");
    for(i=1;i<=n;i++)
	{
		for(j=1;j<=n+1;j++)        /*输入线性方程组的增广矩阵(A,b)*/
			scanf("%f",&A[i][j]);
	}
	printf("\n");   
}
/*****************************************************************************/

/*****************************************************************************/
void main()
{
	Input();
	Choose();	
}
/*****************************************************************************/
/*                                 THE PROGRAMER END                         */
/*****************************************************************************/