加减法(中学所学)是我们平常用的解法之一。 例如，现有如下所示的二元一次方程组。

将等式两边同乘以一个实数，上下系数合并，消去其中一元未知数的方法便是熟知的加减法。

之后，把带入式1，解得。

把上式用行列式表示如下，

之后，第2行乘以，上下相减，得到

的形式。随后，从最下面的式子中解出未知数，代入，得到：

前面的操作是指，把系数行列的左下部分(不包括对角元素)全部变成0，求解如下所示的三角行列(upper triangle matrix)。

这样的处理称作前进消去(forward elimination)。 根据前进消去，未知数可通过解出， 把带入行可以解出，之后把 代入行可以解出 。然后反复执行类似的操作直到解出，这时我们便得到了所有解。这样的后半部分的处理称作后退代入(back substitution)。通过前进消去和后退代入求解n元1次方程组的方法便是高斯消去法。

高斯消去法的代码实现

首先，把常数向量添加到的最后一列组成的行列。

其中第i行j列的元素表示为 

由于要在C,C++上实现，这里的元素序号从0开始计数。

前进消去可以表示成下式

这里的

根据上式，把用更新得到上三角行列。这一步的C++代码如下。

 // 前进消去(forward elimination)//  - 将左下角元素变为0
for(int k = 0; k < n-1; ++k){double akk = A[k][k];for(int i = k+1; i < n; ++i){double aik = A[i][k];for(int j = k+1; j < n+1; ++j){ A[i][j] = A[i][j]-aik*(A[k][j]/akk);}}
}

这里的二维数组A[n][n+1]中存储着系数行列和常数向量。

后退代入的公式如下。

C++代码：

// 后代入(back substitution)
//  - x_n的解为b_n/a_nn，x_n，代入n-1行求出x_(n-1)A[n-1][n] = A[n-1][n]/A[n-1][n-1];
for(int i = n-2; i >= 0; --i){double ax = 0.0;for(int j = i+1; j < n; ++j){ax += A[i][j]*A[j][n];}A[i][n] = (A[i][n]-ax)/A[i][i];
}

计算结果存储于A[i][n]中。

高斯消去法的完整代码如下。

/*!* 高斯消去法(无主元交换)* @param[inout] A 为n×n的系数项与n×1的常数项(b)合并后n×(n+1)的行列* @param[in] n n元一次方程组* @return 1:成功*/
int GaussElimination(vector< vector<double> > &A, int n)
{// 前进消去(forward elimination)//  - 左下角元素为0for(int k = 0; k < n-1; ++k){double akk = A[k][k];for(int i = k+1; i < n; ++i){double aik = A[i][k];for(int j = k+1; j < n+1; ++j){ A[i][j] = A[i][j]-aik*(A[k][j]/akk);}}}// 后退代入(back substitution)A[n-1][n] = A[n-1][n]/A[n-1][n-1];for(int i = n-2; i >= 0; --i){double ax = 0.0;for(int j = i+1; j < n; ++j){ax += A[i][j]*A[j][n];}A[i][n] = (A[i][n]-ax)/A[i][i];}return 1;
}

上面使用STL的vector代替2维数组。

主元选择

我们先重新看一遍高斯消去法中的前进消去公式。

会发现右侧的第二项中包含除以的计算。当前进消去过程中对角元素中有0时，则会停止计算。即使不为0，当有绝对值很小的数出现时也会有误差出现。

针对上述问题，这里我们对主元消元法(pivotal elimination method)进行说明。首先，我们知道联立方程的公式的顺序是可以交换的。也就是说，处理第行时，把其对角元素与同列还未处理的元素 的值进行比较，把绝对值最大的那一行与行互换后，应用到前进消去公式中。上述操作便是主元消元。 像这样的只交换行的操作称为部分主元消元(partial pivoting)， 列方向也交换的操作称为完全主元消元(complete pivoting)。值得注意的是，部分主元消元中未知数的顺序不发生变化，完全主元消元中，未知数的顺序会发生变化。

包含主元消元的高斯消去法代码如下

/*!* 主元选择(Pivoting)*  - 只交换行的部分主元消去法* @param[inout] A * @param[in] n n元一次方程组* @param[in] k 对象行* @return 1:成功*/
inline void Pivoting(vector< vector<double> > &A, int n, int k)
{// 检索k行后k列的绝对值最大元素所在行int p = k; // 绝对值最大的行double am = fabs(A[k][k]);// 最大値for(int i = k+1; i < n; ++i){if(fabs(A[i][k]) > am){p = i;am = fabs(A[i][k]);}}// 若k != p，行交换(主元选择)if(k != p) swap(A[k], A[p]);
}/*!* 高斯消去法(pivoting)* @param[inout] A * @param[in] n n元一次方程组* @return 1:成功*/
int GaussEliminationWithPivoting(vector< vector<double> > &A, int n)
{// 前进消去(forward elimination)for(int k = 0; k < n-1; ++k){// 主元消去Pivoting(A, n, k);double akk = A[k][k];for(int i = k+1; i < n; ++i){double aik = A[i][k];for(int j = k; j < n+1; ++j){A[i][j] = A[i][j]-aik*(A[k][j]/akk);}}}// 后退代入(back substitution)A[n-1][n] = A[n-1][n]/A[n-1][n-1];for(int i = n-2; i >= 0; --i){double ax = 0.0;for(int j = i+1; j < n; ++j){ax += A[i][j]*A[j][n];}A[i][n] = (A[i][n]-ax)/A[i][i];}return 1;
}

高斯･约当法

首先我们假设 是一个 的正方行列。 若 ，则是一个正则行列(regular matrix)(反之为奇异行列(singular matrix))。如果把 的单位行列表示成 ，且 非奇异，则只存在一个 行列 满足 。  这个可以用逆矩阵表示。

对于n元1次方程组，用替换，用替换。增广行列可以表示为

由于，两边同乘，得到

也就是说，增广行列的左侧一半变成了单位行列，右侧一半可以得到逆行列。

此时，逆行列为，

像这样的，扩增前进消去，执行上述处理求得逆行列的方法，便是高斯・约旦法。

高斯・约当法的实现

高斯・约当法公式如下。

求得逆行列代码如下

/*!* 高斯・约当法(无Pivoting)* @param[inout] A n×2n的增广行列* @param[in] n n元一次方程* @return 1:成功*/
int GaussJordan(vector< vector<double> > &A, int n)
{// 将增广行列的右半部分变为单位行列for(int i = 0; i < n; ++i){for(int j = 0; j < n; ++j){A[i][j+n] = (i == j ? 1.0 : 0.0);}}// 根据Gauss-Jordan method计算逆行列for(int k = 0; k < n; ++k){double akk = A[k][k];// 为使得对角元素为1，k行所有元素除以a_kkfor(int j = 0; j < 2*n; ++j){A[k][j] /= akk;}// k列的非对角元素变为0for(int i = 0; i < n; ++i){if(i == k) continue;double aik = A[i][k];for(int j = 0; j < 2*n; ++j){A[i][j] -= A[k][j]*aik;}}}return 1;
}

参考文献

• 佐藤次男, 中村理一郎, “よくわかる数値計算 アルゴリズムと誤差解析の実際”, 日刊工業新聞社, 2001.
• 川上一郎, “数値計算の基礎”, http://www7.ocn.ne.jp/~kawa1/