高斯消元法的C++简单实现

高斯消元法

  首先,我们导入几个概念。

定义1: 一个矩阵称为阶梯形(行阶梯形),若它有以下三个性质:

  1.每一非零行在每一零行之上;

  2.某一行的先导元素所在的列位于前一行先导元素的后面;

  3.某一行先导元素所在列下方元素都是零。

  比如,高斯消元法的C++简单实现

定义2:若一个阶梯形矩阵还满足以下性质,称它为简化阶梯形(简化行阶梯形):

  1.每一非零行的先导元素是1;

  2.每一先导元素1是该元素所在列的惟一非零元素。

  比如,高斯消元法的C++简单实现

定理1:每个矩阵行等价于惟一的简化阶梯形矩阵。即简化阶梯形矩阵是惟一的。


   下面,我们用一个具体例子来说明高斯消元法的主要步骤。

  原矩阵:高斯消元法的C++简单实现

  第一步,由最左的非零列开始,这是一个主元列。主元位置在该列顶端。

高斯消元法的C++简单实现

  第二步,在主元列中选取一个非零元作为主元。若有必要的话,对换两行使这个元素移到主元位置上。

高斯消元法的C++简单实现

  第三步,用倍加行变换将主元下面的元素变成0.

高斯消元法的C++简单实现

  第四步,暂时不管包含主元位置的行以及它上面的各行,对剩下的子矩阵使用上述的三个步骤直到没有非零行需要处理为止。

  对每一行重复上述步骤。

高斯消元法的C++简单实现

高斯消元法的C++简单实现

高斯消元法的C++简单实现

  第五步,由最右面的主元开始,把每个主元上方的各元素变成0.若某个主元不是1,用倍乘变换将它变成1.

高斯消元法的C++简单实现

高斯消元法的C++简单实现

高斯消元法的C++简单实现

高斯消元法的C++简单实现

  最后,我们就得到了原矩阵的简化阶梯形。

  其中,第1~4步称为行化简算法的向前步骤,产生唯一的简化阶梯形的第5步,称为向后步骤。

 


 

C++实现

  我们尝试用C++来实现以上步骤。这里只是简单的实现,也就是用代码描述了上述步骤,没有考虑过多的问题。欢迎大家在评论里指出问题,提出更好的建议,以便于日后改进。

  大概的实现思路就是先实现向前步骤:

  首先,我们对于每一行找到第一个不为零的元素,并且将这一行置为1 * * * *的形式,用这一行乘上倍数加到之后的每一行。

  再实现向后步骤:

  然后,我们从最后一行开始,选择主元,加到之前的每一行上,使得该列的元素都为零。

  最后,我们就完成了化简,得到了简化阶梯形。

  以上算法只是一个粗略实现,主要体现在:

  1.对于主元的选定不够最优;

  2.会出现精度问题;

  3.对于某些情况无法处理。

  先暂时贴上代码,之后有时间再进行优化。

 1 #include <iostream>
 2 #include <cstdio>
 3 
 4 using namespace std;
 5 
 6 int main()
 7 {
 8     double martix[100][100];
 9     int n, m; // n行m列
10 
11     scanf("%d %d", &n, &m);
12 
13     // 输入
14     for(int i = 0; i < n; i++)
15         for(int j = 0; j < m; j++)
16             scanf("%lf", &martix[i][j]);
17 
18     // 向前步骤
19     for(int i = 0; i < n - 1; i++)
20     {
21         // 找主元
22         int pos = 0;
23         for(int j = 0; j < m; j++)
24             if(martix[i][j])
25             {
26                 pos = j;
27                 break;
28             }
29 
30         if(martix[i][pos] != 1 && martix[i][pos] != 0)
31         {
32             double tmp = martix[i][pos];
33             for(int j = pos; j < m; j++)
34             {
35                 martix[i][j] = martix[i][j] / tmp;
36             }
37         }
38         for(int j = i + 1; j < n; j++)
39         {
40             if(!martix[j][pos])
41                 continue;
42             double tmp = martix[j][pos];
43             for(int k = pos; k < m; k++)
44             {
45                 martix[j][k] = martix[j][k] - martix[i][k] * tmp;
46             }
47         }
48     }
49 
50     // 向后步骤
51     for(int i = n - 1; i > 0; i--)
52     {
53         int pos = 0;
54         for(int j = 0; j < m; j++)
55             if(martix[i][j])
56             {
57                 pos = j;
58                 break;
59             }
60 
61         if(martix[i][pos] != 1 && martix[i][pos] != 0)
62         {
63             double tmp = martix[i][pos];
64             for(int j = pos; j < m; j++)
65             {
66                 martix[i][j] = martix[i][j] / tmp;
67             }
68         }
69 
70         for(int j = 0; j < i; j++)
71         {
72             if(!martix[j][pos])
73                 continue;
74             double tmp = martix[j][pos];
75             for(int k = pos; k < m; k++)
76             {
77                 martix[j][k] = martix[j][k] - martix[i][k] * tmp;
78             }
79         }
80     }
81 
82     // 输出
83     for(int i = 0; i < n; i++)
84     {
85         for(int j = 0; j < m; j++)
86             printf("%-10.2f", martix[i][j]);
87         printf("
");
88     }
89     return 0;
90 }