高斯全主元消区法.cpp

这是计算方法上的所有算法

#include <math.h>
#include "stdio.h"
#include "conio.h"
#define N 3
double x[N+1];
static double aa[N][N]={{0.01,2,-0.5},{-1,-0.5,2},{5,-4,0.5}};
static double bb[N]={-5,5,9};

void main()
{
	int i,j;
	double a[N+1][N+1],b[N+1],x[N+1],det;
	int t[N+1];//引入列交换保存x向量的各分量的位置顺序
	double gaussl(double a[][N+1],double b[],double x[],int t[]);

	for(i=1;i<=N;i++)
	{
		for(j=1;j<=N;j++) a[i][j]=aa[i-1][j-1];
		b[i]=bb[i-1];
	}	

	for(i=1;i<=N;i++) t[i]=i;
	det=gaussl(a,b,x,t);
	if(det!=0)
	{	printf("\n方程组的解为:");
		for(i=N;i>=1;i--)
		{
			printf(" x[%d]=%f",t[i],x[i]);
			if(i>1)printf(" -->");

		}
		printf("\n\n系数矩阵的行列式值＝%f",det);
	}
	else printf("\n\n系数矩阵奇异阵，高斯方法失败 !:");
	getch();
			
}

double gaussl(double a[][N+1],double b[],double x[],int t[])
//a传入增广矩阵(有效元素为a[1,1]...a[n,n+1]),x负责传入迭代初值并返回解向量
//(有效值为x[1]...x[n]);返回值：系数矩阵行列式值detA
{double det=1.0,c0,F,m,temp;
 int i,j,k,r,s,n;
 void disp(double a[][N+1],double b[],int x[]);
 printf("消元前增广矩阵 :\n");
 disp(a,b,t);
 for(k=1;k<N;k++)
 { temp=a[k][k];
   r=k;
   s=k;
   for(i=k;i<=N;i++)
	   for(j=k;j<=N;j++)
         if(fabs(temp)<fabs(a[i][j]))
		 {temp=a[i][j];
	      r=i;s=j;
		 }//选主元，选取ik,jk
 if(a[r][s]==0)
	  return 0;
 if(r!=k)
    {for(j=k;j<=N;j++)
      {c0=a[k][j];
       a[k][j]=a[r][j];
       a[r][j]=c0;
       }
    c0=b[k];
    b[k]=b[r];
    b[r]=c0;
    det=-det;
    printf("交换第 %d, %d行:\n",k,r);
    disp(a,b,t);
    }
  if(s!=k)
	  {for(i=0;i<=N;i++){     
	       c0=a[i][k];
           a[i][k]=a[i][s];
           a[i][s]=c0;}
      n=t[k];
	  t[k]=t[s];
	  t[s]=n;
	  det=-det;
      printf("交换第 %d, %d列:\n",k,s);
	  //交换
 //如果ik!=k,则交换[A,b]第k行与第ik行元素
 //如果jk!=k,则交换A的第k列和第jk列元素
      disp(a,b,t);
  }
  for(i=k+1;i<=N;i++)
  {m=a[i][k]/a[k][k];
    for(j=1;j<=k;j++)
		a[i][j]=0.0;
    for(j=k+1;j<=N;j++)
       a[i][j]-=m*a[k][j];
	b[i]-=m*b[k];
  }//消元计算
  //mik=aik/akk (i=k+1,...,n)
  //aij<-aij-mikakj (i,j=k+1,...,n)
  //bi<-bi-mikbk (i=k+1,...,n)
    printf("第%d次消元:\n",k);
    disp(a,b,t);
    
  det*=a[k][k]; 
 }
  
  x[N]=b[N]/a[N][N];
  for(i=N-1;i>0;i--)
    {F=0.0;
     for(j=i+1;j<=N;j++){
		                F+=a[i][j]*x[j];}
	 
     x[i]=(b[i]-F)/a[i][i];
     }//回代计算
   //yn-bn/ann
   //yi=(bi-@@aijyj)/aii (i=n-1,...2,1)

  det*=a[N][N];
 
  return det;
	
}

void disp(double a[][N+1],double b[],int x[])
//显示选主元及消元运算中间增广矩阵结果
{int i,j;
 for(i=1;i<=N;i++)
   {for(j=1;j<=N;j++)
      printf("%f    ",a[i][j]);
    printf(": x%d = %f\n",x[i],b[i]);
    }
}