gauss&doolittle.cpp

用Gauss消去法和Doolittle三角分解法求解解线性方程组的程序源代码．

#include<stdio.h>
#include<math.h>
#define N 100

 float a[N][N],l[N][N],u[N][N],b[N],x[N],y[N];
 int   i,j,k,n,r; 
 
 void shu_ru()
 {
	printf("输入矩阵的维数:n=");		
	scanf("%d",&n);
 printf("输入增广矩阵(A,b):\n");
	for(i=1;i<=n;i++)
	{for(j=1;j<=n+1;j++)
	 scanf("%f",&a[i][j]);
	}printf("\n");
}
/*********************************************/ 
 void gauss()//消元过程
 {  float M;
	for(j=1;j<=n;j++)
		for(i=j+1;i<=n;i++)
		{
		M=a[i][j]/a[j][j];
		for(k=1;k<=n+1;k++)
		a[i][k]=a[i][k]-M*a[j][k];				
		b[i]=b[i]-M*b[i];
		}
    printf("消元之后的矩阵A：\n");
		for(i=1;i<=n;i++)
			for(j=1;j<=n+1;j++){
				printf("%10.4f",a[i][j]);
				if(j==n+1)printf("\n");}
}
/************************************************/
void huidai() /*回代求解过程*/
{
	x[n]=a[n][n+1]/a[n][n];          
   for( k=n-1;k>=1;k--)
   {
     float H=0;
     for(j=k+1;j<=n;j++)
	  {
       H=H+a[k][j]*x[j];
      }
       x[k]=(a[k][n+1]-H)/a[k][k];
   }
}
/*******************************************************/
 //列主元素法
void  xuan_xiao()
{   int p;	
	double max;
    float t,M;
  for(k=1;k<n;k++)
		{
           max=fabs(a[k][k]);//将最大值赋给主对角线上
         for(p=k+1;p<=n;p++) 
			{
				if(max<fabs(a[p][k]))
				{max=fabs(a[p][k]);r=p;}
				else max=fabs(a[p][k]);
			}			    
						//交换行并打印交换后的矩阵
	for(j=1;j<=n+1;j++)
		{t=a[r][j];	a[r][j]=a[k][j];a[k][j]=t;}	
								
printf("\n第%d次交换的结果为：\n",k);
  for(i=1;i<=n;i++)
	for(j=1;j<=n+1;j++) {
		 printf("%10.4f",a[i][j]);
		            if(j%(n+1)==0)printf("\n");	 
				  }  
			 
	  //计算并打印消去的矩阵
		 for(i=k+1;i<=n;i++)
		{ 	
			M=a[i][k]/a[k][k];
				for(j=k;j<=n+1;j++)
				 a[i][j]=a[i][j]-M*a[k][j];}
						
printf("消第%d次得到的矩阵 :\n",k);
        for(i=1;i<=n;i++)
		      for(j=1;j<=n+1;j++)
		{	printf("%10.4f",a[i][j]);
		      if(j%(n+1)==0)printf("\n");}  
		}	  
}	  
/*************************************************/
void san_jiaoshuru()
{
    printf("输入矩阵的维数:n=");		
	scanf("%d",&n);
  printf("输入矩阵A:\n");
	for(i=1;i<=n;i++)
	{for(j=1;j<=n;j++)
	 scanf("%f",&a[i][j]);
	}printf("\n");
}
/************************************************/
void fen_jie()
{ float d,e;
      for(i=1;i<=n;i++)
		for(j=1;j<=n;j++)
	if(i==j) l[i][j]=1.0;//将矩阵L的主对角线赋值为1	
   //将数组A分解成数组U和数组L
	//r代表的是行数,i代表列数
	for(r=1;r<=n;r++)	//当行数为第1行时,分别求得u的第1行和l第1列
		if(r==1){
			for(i=1;i<=n;i++)u[r][i]=a[r][i];//公式4.4
			for(i=2;i<=n;i++)l[i][r]=a[i][r]/u[r][r];//公式4.5
		}
        //当行数为第2到第n时,分别求得u的第2到第n行和l第2到第n列
		else {
			for(i=r;i<=n;i++){
				d=0;
				for(k=1;k<=r-1;k++) d=d+l[r][k]*u[k][i];
				u[r][i]=a[r][i]-d;}//公式4.6			
			   for(i=r+1;i<=n;i++)if(r!=n){
				e=0;
				for(k=1;k<=r-1;k++)e=e+l[i][k]*u[k][r];
				l[i][r]=(a[i][r]-e)/u[r][r];}//公式4.7
         	}
	printf("输出矩阵L:\n");
		for(i=1;i<=n;i++)
			for(j=1;j<=n;j++){
				printf("%10.4f",l[i][j]);
				if(j==n)printf("\n");}
			printf("\n");
         
	printf("输出矩阵U:\n");
           for(i=1;i<=n;i++)
			 for(j=1;j<=n;j++){
				printf("%10.4f",u[i][j]);
				if(j==n)printf("\n");}
			 printf("\n");
}
/***************************************************/
void shuru__b()
{
	printf("输入矩阵b:\n");
	for(i=1;i<=n;i++)
	scanf("%f",&b[i]);
	printf("\n");
}
/**********************************************************/
void  qiuzhi()
{  float f,g;
   for(r=1;r<=n;r++)//LY=B求得数组Y
		if(r==1)y[r]=b[r];//公式4.10
		else {
			f=0;
			for(i=1;i<=r-1;i++) f=f+l[r][i]*y[i];
		    y[r]=b[r]-f;}//公式4.10
  
	for(r=n;r>=1;r--)//UX=Y求得数组X
		if(r==n) x[r]=y[r]/u[n][n];//4.11
		else {
			g=0;
			for(i=r+1;i<=n;i++)g=g+u[r][i]*x[i];
			x[r]=(y[r]-g)/u[r][r];}//4.11
}
/****************************************/
 shu_chu() //输出解
 {    
   printf("该线性方程组的解为:\n");
   for(i=1;i<=n;i++)
   {
     printf(" \nx%d=%f\n",i,x[i]); 
   }
   printf("\n");
    return i;
}
/************************************************************/
 //主函数
void main()
{  int m;
 loop:
    printf("     1,Gauss消去法\n");
	printf("     2,列主元素消去法\n");
	printf("     3,Doolittle三角分解法\n");
	printf("选择解线性方程组的方法：");
	scanf("%d",&m);
switch(m){
case 1:printf("\nGauss消去法:\n");
	   shu_ru();     //输入矩阵
	   gauss();        //Gauss消去法
       huidai();      //回代
	   shu_chu();     //输出结果x1,x2,....,xn	
	 goto loop;

case 2:printf("\n列主元素消去法:\n");
       shu_ru();    //输入矩阵
	   xuan_xiao();   //选主元,换行,消元
	   huidai();      //回代
	   shu_chu();    	 //输出结果x1,x2,....,xn	
	 goto loop;

case 3:printf("\n三角分解法:\n");
	   san_jiaoshuru(); //输入方阵
	   fen_jie();      //分解矩阵
	   shuru__b();
	   qiuzhi();          
	   shu_chu(); //输出结果x1,x2,....,xn
	goto loop;
	default: break;
}
}