线性方程组的直接解法

线性方程组的直接解法

一、实验目的

1. 掌握高斯消去法、直接三角分解法求解线性方程组的基本原理以及线性方程组条件数的计算方法；
2. 已知线性方程组 Ax=b，其中
   $$A=\left[\begin{array}{cccc}1& 2& -12& 8\\ 5& 4& 7& -2\\ -3& 7& 9& 5\\ 6& -12& -8& 3\end{array}\right]b=\left[\begin{array}{c}27\\ 4\\ 11\\ 49\end{array}\right]$$

分别编程实现用高斯消去法、直接三角分解法求解以上线性方程组的根。

二、实验设备

1. 编程语言：C；
2. 实验环境：已安装相关编程环境的计算机1台。

三、算法实现及结果分析

#define  _CRT_SECURE_NO_WARNINGS
#include<stdio.h>
double A[10][10];//系数矩阵
double B[10];//常数矩阵
double X[10];//根矩阵
double L[10][10];//下三角矩阵
double U[10][10];//上三角矩阵
double y[10];
double x[10];
int n;
//高斯消去法
void gausselimination() {
	int i, j, k;
	double AB[10][10];
	//获得增广矩阵
	for (i = 0;i < n;i++) {
		for (j = 0;j < n;j++) {
			AB[i][j] = A[i][j];
		}
		AB[i][n] = B[i];
	}
	//求得阶梯矩阵
	double l;
	for (i = 0;i < n - 1;i++) {
		for (j = i;j < n - 1;j++) {
			l = -AB[j + 1][i] / AB[i][i];
			for (k = i;k <= n;k++) {
				AB[j + 1][k] += AB[i][k]*l;
			}
		}
	}
	for (i = 0;i < n;i++) {
		for (j = 0;j <= n;j++) {
			printf("%lf ", AB[i][j]);
		}
		printf("\n");
	}
	//根据阶梯矩阵求取线性方程组的根
	for (i = n - 1;i >= 0;i--) {
		double sum = 0;
		for (j = i + 1;j < n;j++) {
			sum += AB[i][j] * X[j];
		}
		X[i] = (AB[i][n] - sum) / AB[i][i];
	}
	printf("高斯消去法求得线性方程组的根为：");
	for (int i = 0;i < n;i++) {
		printf("X(%d)=%lf  ", i + 1, X[i]);
	}
	printf("\n");
}

//直接三角分解法(L的对角元为1)
void LU1(){
	int i,j,k,r;
	for (i = 0;i < n;i++) {
		for (j = 0;j < n;j++) {
			L[i][j] = 0;
			U[i][j] = 0;
		}
		x[i] = 0;
		y[i] = 0;
	}
	//求U的第一行元素和L的第一列元素
	for(i=0;i<n;i++){
		U[0][i]=A[0][i];
		L[i][0]=A[i][0]/U[0][0];
	}
	//求U的第二行元素
	for(j=1;j<n;j++){
		U[1][j]=A[1][j]-L[1][0]*U[0][j];
	}
	//求L的第二列元素
	for(i=1;i<n;i++){
		L[i][1]=(A[i][1]-L[i][0]*U[0][1])/U[1][1];
	}
	//下三角矩阵的对角元为1
	for(k=0;k<n;k++){
		L[k][k]=1;
	}
	for(k=2;k<n;k++){
		for(j=k;j<n;j++){
			U[k][j]=A[k][j];
			for(r=0;r<k;r++){
				U[k][j]-=L[k][r]*U[r][j];
			}
		}
		for(i=k+1;i<n;i++){
			L[i][k]=A[i][k];
			for(r=0;r<k;r++){
				L[i][k]-=L[i][r]*U[r][k];
			}
			L[i][k]/=U[k][k];
		}
	}
	printf("矩阵L为：\n");
	for(i=0;i<n;i++){
		for(j=0;j<n;j++){
			printf("%lf ",L[i][j]);
		}
		printf("\n");
	} 
	printf("矩阵U为：\n");
	for(i=0;i<n;i++){
		for(j=0;j<n;j++){
			printf("%lf ",U[i][j]);
		}
		printf("\n");
	} 
	for(k=0;k<n;k++){
		y[k]=B[k];
		for(j=0;j<k;j++){
			y[k]-=L[k][j]*y[j];
		}
	}
	x[n-1]=y[n-1]/U[n-1][n-1];
	for(k=n-2;k>=0;k--){
		x[k]=y[k];
		for(j=k+1;j<n;j++){
			x[k]-=U[k][j]*x[j];
		}
		x[k]/=U[k][k];
	}
	printf("直接三角分解法(L对角元素为1)求得线性方程组的根为："); 
	for(i=0;i<n;i++){
		printf("X(%d)=%f  ", i + 1, x[i]);
	}
	printf("\n");
}
//直接三角分解法(U的对角元为1)
void LU2() {
	int i, j, k, r;
	for (i = 0;i < n;i++) {
		for (j = 0;j < n;j++) {
			L[i][j] = 0;
			U[i][j] = 0;
		}
		x[i] = 0;
		y[i] = 0;
	}
	//求L的第一行元素和U的第一列元素
	for (i = 0;i < n;i++) {
		L[i][0] = A[i][0];
		U[0][i] = A[0][i] / L[0][0];
	}
	//求L的第二列元素
	for (i = 1;i < n;i++) {
		L[i][1] = A[i][1] - U[0][1] * L[i][0];
	}
	//求U的第二行元素
	for (j = 1;j < n;j++) {
		U[1][j] = (A[1][j] - U[0][j] * L[1][0]) / L[1][1];
	}
	//上三角矩阵的对角元为1
	for (k = 0;k < n;k++) {
		U[k][k] = 1;
	}
	for (k = 2;k < n;k++) {
		for (j = k;j < n;j++) {
			L[j][k] = A[j][k];
			for (r = 0;r < k;r++) {
				L[j][k] -= L[j][r] * U[r][k];
			}
		}
		for (i = k + 1;i < n;i++) {
			U[k][i] = A[k][i];
			for (r = 0;r < k;r++) {
				U[k][i] -= U[r][i] * L[k][r];
			}
			U[k][i] /= L[k][k];
		}
	}
	printf("矩阵L为：\n");
	for (i = 0;i < n;i++) {
		for (j = 0;j < n;j++) {
			printf("%lf ", L[i][j]);
		}
		printf("\n");
	}
	printf("矩阵U为：\n");
	for (i = 0;i < n;i++) {
		for (j = 0;j < n;j++) {
			printf("%lf ", U[i][j]);
		}
		printf("\n");
	}
	for (k = 0;k < n;k++) {
		y[k] = B[k];
		for (j = 0;j < k;j++) {
			y[k] -= L[k][j] * y[j];
		}
		y[k] /= L[k][k];
	}
	x[n - 1] = y[n - 1] / U[n - 1][n - 1];
	for (k = n - 2;k >= 0;k--) {
		x[k] = y[k];
		for (j = k+1;j < n;j++) {
			x[k] -= U[k][j] * x[j];
		}
	}
	printf("直接三角分解法(L对角元素为1)求得线性方程组的根为：");
	for (i = 0;i < n;i++) {
		printf("X(%d)=%f  ", i + 1, x[i]);
	}
}
int main() {
	int i,j;
	printf("请输入矩阵维度n=\n");
	scanf("%d", &n);
	printf("请输入系数矩阵\n");
	for (i = 0;i < n;i++) {
		for (j = 0;j < n;j++) {
			scanf("%lf", &A[i][j]);
		}
	}
	printf("请输入常数矩阵\n");
	for (i = 0;i < n;i++) {
		scanf("%lf", &B[i]);
	}
	printf("高斯消去法:");
	gausselimination();
	printf("直接三角分解法:");
	LU1();
	LU2();
	return 0;
}

经过运行，我们可以看到成功求得的线性方程组解及L，U矩阵。在直接三角分解法中，书中给的例子是L矩阵的对角元素为1，我们考虑将U矩阵的对角元素设为1，最终也是成功求得L，U矩阵及线性方程组的解。高斯消元法是将增广矩阵变换成同解的上三角矩阵，我们或许可以考虑将进一步将矩阵转化成单位阶梯矩阵进行求解。在直接三角矩阵分解法中，我们需要注意L，U矩阵元素的求解顺序。