ex_1.cpp

计算方法中列主消元,GUASS法,三角法的算法

//求解线性方程组,采用三种方法:Guass消去法、Guass列主元素消去法
// 和直接三角形分解法(Boolittle方法 \ Crout方法)
// 针对课后第五题作个误差分析

#include"stdio.h"
#include"stdlib.h"
#include"math.h"

//求系数矩阵行列式
void det_A(float *A,int n)
{

}

//对增广矩阵(A,b)作一次消元，即(A(row),b(row)) -> (A(row+1),b(row+1))
void CountA_b(float *a,float *b,int row,int n)
{
	int i,j;

	for(i=row+1;i<=n-1;i++)
	{
		for(j=row+1;j<=n-1;j++)
			*(a+i*n+j)=*(a+i*n+j)-*(a+row*n+j)*(*(a+i*n+row))/(*(a+row*n+row));
		*(b+i)=*(b+i)-*(b+row)*(*(a+i*n+row))/(*(a+row*n+row));
	}
}

//使 A(n) 变为上三角矩阵
void Free_A(float *a,int n)
{
	int i,j;

	for(i=1;i<=n-1;i++)
	for(j=0;j<i;j++)
      *(a+i*n+j)=0;
}

//找出从a[row][row]-->a[n-1][row]中最大元素所在行，
//也即第row列中最大元素所在行
int Max_Row(float *a,int row,int n)
{
	int i,max=row;

	for(i=row+1;i<=n-1;i++)
	  if(fabs(*(a+max*n+row))<fabs(*(a+i*n+row))) max=i;

	return max;
}

//交换增广矩阵(A(r),b(r))第 r 行与第 k 行的元素，包括 b 数组
void Exchange_Row(float *a,float *b,int r ,int k,int n)
{
	int j=r;
	float temp;

	for( ;j<=n-1;j++)
	{
	  temp=*(a+r*n+j);*(a+r*n+j)=*(a+k*n+j);*(a+k*n+j)=temp;
	}

	temp=*(b+r);*(b+r)=*(b+k);*(b+k)=temp;
}

//Guass消去法
void Guass(float *a,float *b,int n)
{
	int k;

	for(k=0;k<=n-2;k++)
	{
	  if(*(a+k*n+k)==0) exit(0); //主元素 a[k][k] 不为零
	  CountA_b(a,b,k,n);      //每进行消去一次，判断一次主元素a[k][k]
	}

    //Free_A(a,n);
}

//Guass列主元素消去法
void Guass_Colum(float *a,float *b,int n)
{
	int k,max;

	for(k=0;k<=n-2;k++)
	{
	  if(*(a+k*n+k)==0) exit(0);    //主元素 a[k][k] 不为零
	  max=Max_Row(a,k,n);        //查找最大元素所在的行
	  if(k!=max) Exchange_Row(a,b,k,max,n);
	  CountA_b(a,b,k,n);         //进行消去一次
	}

	//Free_A(a,n);
}

//上三角形回代过程 兼顾 单位上三角形回代
void Return_Top(float *a,float *b,float *x,int n,char p)
{
	int i,j;
	float temp;

	for(i=n-1;i>=0;i--)
	{
	  temp=0;
	  for(j=i+1;j<=n-1;j++)
		temp+=*(a+i*n+j) * (*(x+j));
	  if(p=='Y') *(x+i)=*(b+i)-temp;
	   else *(x+i)=(*(b+i)-temp)/(*(a+i*n+i));
	}
}

//下三角形回代过程 兼顾 单位下三角形回代
void Return_Bottom(float *a,float *b,float *x,int n,char p)
{
	int i,j;
	float temp;

	for(i=0;i<=n-1;i++)
	{
		temp=0;
		for(j=0;j<=i-1;j++)
		  temp+=*(a+i*n+j)*(*(x+j));
		if(p=='Y') *(x+i)=*(b+i)-temp;
		else *(x+i)=(*(b+i)-temp)/(*(a+i*n+i));
	}
}

//直接三角形分解法中Boolittle方法
void Doolittle(float *a,float *LU,int n)
{
	int i,k,p;
	float temp;

	for(k=0;k<=n-1;k++)
	{
		for(i=k;i<=n-1;i++)
		{
			temp=0;
			for(p=0;p<=k-1;p++)
				temp+=*(LU+k*n+p)*(*(LU+p*n+i));
			*(LU+k*n+i)=*(a+k*n+i)-temp;
		}

		for(i=k+1;i<=n-1 && k!=n-1;i++)
		{
			temp=0;
			for(p=0;p<=k-1;p++)
				temp+=*(LU+i*n+p)*(*(LU+p*n+k));
			*(LU+i*n+k)=(*(a+i*n+k)-temp)/(*(LU+k*n+k));
		}
	}
}

//直接三角形分解法中Crout方法
void Crout(float *a,float *LU,int n)
{
	int i,k,p;
	float temp;

	for(k=0;k<=n-1;k++)
	{
	    for(i=k;i<=n-1;i++)
		{
			temp=0;
			for(p=0;p<=k-1;p++)
			temp+=*(LU+i*n+p)*(*(LU+p*n+k));
			*(LU+i*n+k)=*(a+i*n+k)-temp;
		}

		for(i=k+1;i<=n-1 && k!=n-1;i++)
		{
			temp=0;
			for(p=0;p<=k-1;p++)
			temp+=*(LU+k*n+p)*(*(LU+p*n+i));
			*(LU+k*n+i)=(*(a+k*n+i)-temp)/(*(LU+k*n+k)); 
		}
	}
}

//为系数矩阵 A 申请空间
float* MallocSpace_A(int n)
{
	float *base;
	base=(float *)malloc(n*n*sizeof(float));

	return base;
}

//为系数矩阵 b 申请空间
float* MallocSpace_b(int n)
{
	float *base;
	base=(float *)malloc(n*sizeof(float));

	return base;
}

//为分解的上、下三角形矩阵 LU 申请空间
float* MallocSpace_LU(int n)
{
	float *base;
	base=(float *)malloc(n*n*sizeof(float));

	return base;
}

//为解 X 申请空间
float* MallocSpace_X(int n)
{
	float *base;
	base=(float *)malloc(n*sizeof(float));

	return base;
}

//为解 Std_X 申请空间
float* MallocSpace_Std_X(int n)
{
	float *base;
	base=(float *)malloc(n*sizeof(float));

	return base;
}

//为暂存量 Y 申请空间
float* MallocSpace_Y(int n)
{
    float *base;
	base=(float *)malloc(n*sizeof(float));

	return base;
}

//输入系数矩阵和常量矩阵
void InputElement(float *a1,float *b1,float *std_x,int n)
{
	int i,j;
	float temp;

	printf("输入系数矩阵:\n");
	for(i=0;i<=n-1;i++)
	for(j=0;j<=n-1;j++)
    {
		scanf("%f",&temp);
		*(a1+i*n+j)=temp;
	}

	printf("输入常量矩阵:\n");
	for(i=0;i<=n-1;i++)
	{
		scanf("%f",&temp);
		*(b1+i)=temp;
	}

	printf("输入标准的解:\n");
	for(i=0;i<=n-1;i++)
	{
		scanf("%f",&temp);
		*(std_x+i)=temp;
	}
}

//为系数矩阵作个备份
float* Copy_A(float *a1,int n)
{
  int i,j;
  float *base;
  
  base=(float *)malloc(n*n*sizeof(float));

  for(i=0;i<=n-1;i++)
  for(j=0;j<=n-1;j++)
  *(base+i*n+j)=*(a1+i*n+j);

  return base;
}

//为常量矩阵作个备份
float* Copy_b(float *b1,int n)
{
  int i;
  float *base;
  
  base=(float *)malloc(n*sizeof(float));

  for(i=0;i<=n-1;i++)
  *(base+i)=*(b1+i);

  return base;

}

//打印求解线性方程组的系数矩阵
void OutputElement_A(float *A,int n)
{
	int i,j;

	printf("\n此方程组的系数矩阵为:\n");
	for(i=0;i<=n-1;i++)
	{
	  for(j=0;j<=n-1;j++)
	  printf("%f   ",*(A+i*n+j));

	  printf("\n");
	}
}

//打印求解线性方程组的常量矩阵
void OutputElement_b(float *b,int n)
{
	int i;

	printf("\n此方程组的常量矩阵为:\n");
	for(i=0;i<=n-1;i++)
	printf("%f   ",*(b+i));

	printf("\n");

}

//打印求解线性方程组的解
void OutputElement(float *x,int n)
{
	int i;

	printf("\n此方程组的解为:\n");
	for(i=0;i<=n-1;i++)
	printf("%f\n",*(x+i));

}

//为系数矩阵重新赋值
void Give_A(float *a1,float *a2,int n)
{
	int i,j;

	for(i=0;i<=n-1;i++)
	for(j=0;j<=n-1;j++)
	  *(a2+i*n+j)=*(a1+i*n+j);

}

//为常量矩阵重新赋值
void Give_b(float *b1,float *b2,int n)
{
	int i;

	for(i=0;i<=n-1;i++)
	  *(b2+i)=*(b1+i);

}

// 针对课后第五题作个误差分析
void Count_Error(float *x,float *std_x,int n)
{
	int i;
    float temp;

	printf("课后第五题误差分析情况:\n");
    for(i=0;i<n;i++)
	{
      temp=(float)(fabs(*(x+i)-*(std_x+i))/fabs(*(std_x+i)));
	  printf("%.4f",temp*100);
	  putchar('%');        //不能用printf("%");语句
	  printf("   ");
	}
	printf("\n");
}

//选择方法
void Select_Methom(float *a,float *b,float *lu,float *x,float *std_x,float *y,int n)
{
	int select1,select2;

	printf("\n请选择解线性方程组的方法,提示信息如下:\n");
	printf("\n选择 Guass 消去法,          请输入 1\n");
	printf("选择 Guass 列主元素消去法,  请输入 2\n");
	printf("选择 直接三角形分解法,      请输入 3\n");
	printf("选择 退出,                  请输入 4\n");
	printf("\n请选择你所要采用的方法:  ");
	scanf("%d",&select1);

	switch(select1)
	{
	case 1:printf("\n--------------------采用Guass消去法--------------------\n");
		Guass(a,b,n);
		Return_Top(a,b,x,n,'N');
		OutputElement(x,n);
		break;
	case 2:printf("\n--------------------采用Guass列主元素消去法--------------------\n");
		Guass_Colum(a,b,n);
		Return_Top(a,b,x,n,'N');
		OutputElement(x,n);
		break;
    case 3:printf("\n--------------------采用直接三角形分解法--------------------\n");
		printf("请选择分解的方法 :\n");
		printf("选择 Boolittle 方法,请输入 1\n");
		printf("选择 Crout 方法,请输入 2\n");
		printf("采用哪种方法: ");
		scanf("%d",&select2);

		switch(select2)
		{
		case 1:printf("--------------------Bolittle--------------------\n");
			Doolittle(a,lu,n);
			Return_Bottom(lu,b,y,n,'Y');//单位下三角形回代
			Return_Top(lu,y,x,n,'N');
            OutputElement(x,n);
			break;
		case 2:printf("--------------------Crout--------------------\n");
			Crout(a,lu,n);
			Return_Bottom(lu,b,y,n,'N');
			Return_Top(lu,y,x,n,'Y');//单位上三角形回代
			OutputElement(x,n);
			break;
		}
		break;
	case 4:exit(0);
	}
	// 针对课后第五题作个误差分析
    Count_Error(x,std_x,n);
}
  
//主函数
void main()
{
	int n;  
	char flag='Y';   //flag 用来表示求解过程是否要继续，开始为 'Y'
	float *a1,*a2,*b1,*b2;
	float *lu,*x,*std_x,*y;

	printf("----------本程序是用来求解线性方程组的解的!!----------\n");
	printf("请输入求解线性方程组的阶数(n): ");
	scanf("%d",&n);

	a1=MallocSpace_A(n);
	b1=MallocSpace_b(n);
	lu=MallocSpace_LU(n);
	x=MallocSpace_X(n);
	std_x=MallocSpace_Std_X(n);
	y=MallocSpace_Y(n);

	InputElement(a1,b1,std_x,n);

	a2=Copy_A(a1,n);
	b2=Copy_b(b1,n);

	while(flag=='Y')
	{
      Give_A(a1,a2,n);
	  Give_b(b1,b2,n);
	  Select_Methom(a2,b2,lu,x,std_x,y,n);
	  getchar();     //去消输入 n 时,产生的回车符'\n',使scanf()语句不再继续执行
	  printf("\n是否还需须尝试用其它的方法进行求解(Y:继续;任意:退出): ");
      scanf("%c",&flag);
	}
}