数值分析之线性方程组直接解法（一）

66 阅读 0 评论 44 点赞

我是靠谱客的博主优秀雪碧，最近开发中收集的这篇文章主要介绍数值分析之线性方程组直接解法（一），觉得挺不错的，现在分享给大家，希望可以做个参考。

概述

数学上，高斯消元法（或译：高斯消去法）（Gaussian Elimination），是线性代数中的一个算法，可用来为线性方程组求解，求出矩阵的秩，以及求出可逆方阵的逆矩阵。

高斯消元法的原理

（一）计算过程

高斯消去法就是通过矩阵的行变换达到消元的目的，从而将方程组的系数矩阵由对称矩阵变为三角矩阵，最后获得方程组的解。

线性方程组

假设方程组的系数矩阵A非奇异（大致意思就是方程组有非零解的条件，具体定义可百度），我们要将系数矩阵A化成上三角（或下三角亦可），再回代求解。

1）首先构造系数矩阵的增广矩阵，然后再开始进行消元，如下：

第一次消元

第二次消元

第三次消元

总的计算量

2）其次是回代求解方程组，如下：

回代求解

这样就可以求得线性方程组的解了。

（二）具体例子：

线性方程组为：


2x + y - z = 8
（L1）
-3x - y + 2z = -11
（L2）
-2x + y + 2z = -3
（L3）

1. 首先第一步是构造增广矩阵，如下：


2
1
-1
|
8
-3
-1
2
|
-11
-2
1
2
|
-3

2. 第二步是利用第一行消去第二行与第三行的第一列，如下：

3. 第三步是利用第二行消去第三行的第二列，如下：

4. 第四步是回代求结果，如下：


1） 根据第三个方程可直接求得 z = -1；
2） 将z代入第二个方程可得 y = 3；
3） 再将y和z代入第一个方程可得 x = 2。由此就可以得到该方程组的解为：
（x = 2，y = 3，z = -1）

（三）存在的问题及改进

在以上的高斯消元法中，消元能进行的条件是主元素（即主对角线上的元素）不能等于0，即a(i,i)，（i = 1,2…）不等于0。在数学理论上，只要满足前面的条件，都可以求出线性方程组的解。但是在计算机程序实现的过程中，还需要考虑其他因素，因为计算机内部无法存储无限小数，必然会有舍入误差，因此有时虽然a(i,i)不等于0，但是绝对值|a(i,i)|很小，这时计算过程的舍入误差会导致消去法数值不稳定，以致结果不可靠。例如下面这个例子：

线性方程组

计算结果

为了解决这个问题，这里进入一个改进的计算方法，列主元高斯消元法。列主元高斯消元法的计算原理与高斯消元法的原理大致相同，只是在每一次消元时要先选择以哪一行元素作为标准进行消元，列主元高斯消元法选择的是以绝对值最大的那一行作为标准行来进行消元，这么说可能不太清楚，具体看下面的例子就可以明白了：

1）首先构造增广矩阵为

2）其次进行选择列主元并开始进行消元
第一次消元选择列主元为a(3,1) = 3.996，交换第一行和第三行，再进行消元计算得，

第二次消元选择列主元为a(3,2) = 2.0028，交换第二行与第三行，再进行消元计算得，

消元过程至此结束，回代计算依次得到解为：
（x3 = 0.90043，x2 = -0.69850，x1 = 1.9273）

高斯消元法的代码实现

在这里我给出一份我自己写的高斯消去法和列主元高斯消去法的代码，这两个的实现基本都差不多，因此我把它们写在一起了，里面有一个变量可以控制是使用高斯法还是列主元法。
这代码里面有一个误差范围的变量，这个变量是用来判断当前选择的主元素是否为0的取值范围，当有小值主元素时，如果设置的太小可能会影响计算结果。如上面高斯消元法存在问题的那个例子，如果设置误差范围太小，计算结果可能会完全不正确。当然这也要看方程组是不是会出现这种小值主元素，如果没有小值主元素，这个误差范围应该不会影响计算结果。其他的注释写的应该挺清楚的，这里就不多说了。下面看代码：

#include<cstdio>
#include<cstdlib>
#include<cmath>
/* 列主元高斯消去 */
double select_column(int n, int k);
/* 判断是顺序高斯还是列主元高斯，默认是顺序高斯 */
bool is_column = false;
/* 测试输出 */
bool is_test = false;
/* 各矩阵及相关向量 */
double matrix_A[100][100];
double vector_b[100]= {0},vector_x[100]= {0};
int main()
{
int n;
//矩阵的行列数，行列数相同
double ison;
// 误差值
printf("input the size of matrix( 0<size<100 ):n");
scanf("%d", &n);
if (n <= 0 || n > 100) return -1;
int i, j, k;
for (i = 0; i < n; i++)
//初始化数据
{
for (j = 0; j < n; j++)
matrix_A[i][j] = 0;
}
printf("input the param matrix:n");
//输入系数矩阵
for (i = 0; i < n; i++)
for (j = 0; j < n; j++)
scanf("%lf", &matrix_A[i][j]);
printf("input the result vector:n");
//输入结果向量
for (i = 0; i < n; i++)
scanf("%lf", &vector_b[i]);
printf("input the error range:n");
//输入结果的误差范围
scanf("%lf", &ison);
/* 求解高斯矩阵 */
double temp;
for(k=0; k<n-1; k++)
{
/* 列主元消元使用 */
if (is_column)
{
double max_s = select_column(n, k);
printf("nthe current max: %lf", max_s);
}
/* 误差范围判断 */
if(fabs(matrix_A[k][k])<=ison)
{
printf("A[k][k] is equal to 0,return error!n");
return -1;
}
for(i=k+1; i<n; i++)
{
temp=matrix_A[i][k]/matrix_A[k][k];
vector_b[i]-=temp*vector_b[k];
for(j=k; j<n; j++)
{
matrix_A[i][j]-=temp*matrix_A[k][j];
}
}
if(is_test) //测试
{
printf("n%d the change matrix A:n",k+1);
for (i = 0; i < n; i++)
{
for (j = 0; j < n; j++)
printf("%lf ", matrix_A[i][j]);
printf("n");
}
printf("%d the change vector B:n",k+1);
for (i = 0; i < n; i++)
printf("%lf ", vector_b[i]);
printf("n");
}
}
if(fabs(matrix_A[n-1][n-1])<=ison)
{
printf("A[k][k] is equal to 0,return error!n");
return -1;
}
/* 回代求解未知变量 */
vector_x[n-1]=vector_b[n-1]/matrix_A[n-1][n-1];
for(i=n-2; i>=0; i--)
{
temp=0;
for(j=i+1; j<n; j++)
temp+=matrix_A[i][j]*vector_x[j];
vector_x[i]=(vector_b[i]-temp)/matrix_A[i][i];
}
/* 打印结果向量X */
printf("nthe vector X:n");
for (i = 0; i < n; i++)
printf("x%d=%lf ", i+1, vector_x[i]);
printf("n");
return 0;
}
double select_column(int n, int k)
{
int i, j;
double max_sum = 0;
int index = k;
/* 找出矩阵A剩余行列的主元 */
for (i = k; i < n; i++)
{
/* 记录主元的值及其所在行 */
if (fabs(matrix_A[i][k]) > max_sum)
{
max_sum = fabs(matrix_A[i][k]);
index = i;
}
}
/* 交换矩阵的两行，及结果向量的两个元素 */
if (index != k)
{
double temp;
/* 交换A矩阵的两行 */
for (j = 0; j < n; j++)
{
temp = matrix_A[index][j];
matrix_A[index][j] = matrix_A[k][j];
matrix_A[k][j] = temp;
}
/* 交换向量b的两个数值 */
temp = vector_b[index];
vector_b[index] = vector_b[k];
vector_b[k] = temp;
}
return matrix_A[k][k];
}