【计算方法】#01 高斯消去法和列主元高斯消去法的原理简介及C++实现

1. 高斯消去法
- 1.1 算法的适用条件
- 1.2 算法步骤和公式
- 1.3 算法复杂度分析
- 1.4 算法的C++实现
2 列主元高斯消去法
- 2.1 经典方法的致命问题
- 2.2 按列选主元步骤的算法描述
- 2.3 算法复杂度分析
- 2.4 算法优势
- 2.5 算法的C++实现
References

求解方程组：Ax=bAx=bAx=b

1. 高斯消去法

1.1 算法的适用条件

满足以下条件中的任一即可

系数矩阵AAA的各阶顺序主子式均不等于零(充要)；
系数矩阵AAA是对称正定矩阵；
系数矩阵AAA是严格对角占优矩阵，即对角线元素大于对应行的其余元素之和。

1.2 算法步骤和公式

消元过程(第kkk次消元)：
{aij(k)=aij(k−1)−aik(k−1)akk(k−1)akj(k−1)=aij(k−1)−likakj(k−1),i=k+1,k+2,...,n;j=k+1,k+2,...,nbi(k)=bi(k−1)−likbk(k−1),i=k+1,k+2,...,nlik≜aik(k−1)akk(k−1),i=k+1,k+2,...,n\begin{cases} &a_{ij}^{(k)}=a_{ij}^{(k-1)}-\frac{a_{ik}^{(k-1)}}{a_{kk}^{(k-1)}}a_{kj}^{(k-1)}=a_{ij}^{(k-1)}-l_{ik}a_{kj}^{(k-1)}, \\ &i=k+1,k+2,...,n; j=k+1,k+2,...,n \\ \\ &b_i^{(k)} = b_i^{(k-1)}- l_{ik}b_{k}^{(k-1)}, i=k+1,k+2,...,n \\ \\ &l_{ik} \triangleq \frac{a_{ik}^{(k-1)}}{a_{kk}^{(k-1)}}, i=k+1, k+2, ..., n \end{cases} ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧aij(k)=aij(k−1)−akk(k−1)aik(k−1)akj(k−1)=aij(k−1)−likakj(k−1),i=k+1,k+2,...,n;j=k+1,k+2,...,nbi(k)=bi(k−1)−likbk(k−1),i=k+1,k+2,...,nlik≜akk(k−1)aik(k−1),i=k+1,k+2,...,n

回代过程：
{xn=bn(n−1)ann(n−1),xk=bk(k−1)−∑j=k+1nakj(k−1)xjakk(k−1),k=n−1,n−2,...,2,1\left\{ \begin{aligned} x_n &= \frac{b_n^{(n-1)}}{a_{nn}^{(n-1)}},\\ x_k &= \frac{b_k^{(k-1)}-\sum\limits_{j=k+1}^{n}a_{kj}^{(k-1)}x_j}{a_{kk}^{(k-1)}}, k=n-1,n-2,...,2,1 \end{aligned} \right. ⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧xnxk=ann(n−1)bn(n−1),=akk(k−1)bk(k−1)−j=k+1∑nakj(k−1)xj,k=n−1,n−2,...,2,1

1.3 算法复杂度分析

令运算量为乘除法个数，则

消元过程：N1=∑k=1n−1(n−k)2+2(n−k)=n3/3+n2/2−5n/6N_1=\sum\limits_{k=1}^{n-1}(n-k)^2+2(n-k)=n^3/3+n^2/2-5n/6N1=k=1∑n−1(n−k)2+2(n−k)=n3/3+n2/2−5n/6

回代过程：N2=1+(n−1)+∑k=1n−1(n−k)=n2/2+n/2N_2=1 + (n-1) + \sum\limits_{k=1}^{n-1}(n-k)=n^2/2 + n/2N2=1+(n−1)+k=1∑n−1(n−k)=n2/2+n/2

高斯消去法总运算量：N=N1+N2=n3/3+n2−n/3N=N_1+N_2=n^3/3+n^2-n/3N=N1+N2=n3/3+n2−n/3

时间复杂度：O(n3)O(n^3)O(n3)

空间复杂度：O(n2)O(n^2)O(n2)，开销为系数矩阵AAA和一维数组bbb

1.4 算法的C++实现

typedef float elem;
vector<elem> Gauss_Elimination(vector<vector<elem>>& A, vector<elem>& b, elem epsilon){/* 高斯消去法参数简介：A - n维的系数矩阵b - 方程组等号右边的向量，维度1xnepsilon - 矩阵元素的阈值，当主元过小时高斯消去法将数值不稳定使用前的约定：1. 确保A不为空，且Ax=b中各个参数的维度一致注：执行完后，A和b将会改变*/int n = A.size();vector<elem> x(n);// 消元过程for (int k=0; k<n-1; ++k){ // 第k+1步消元// 主元太小，不能继续if(abs(A[k][k]) < epsilon){cout<<"The pivot element is too small. Process aborted."<<endl;return {};}for (int i=k+1; i<n; ++i){ // l[i][k]A[i][k] /= A[k][k];}// b[i]for (int i=k+1; i<n; ++i){b[i] -= A[i][k]*b[k];}// A[i][j]for (int i=k+1; i<n; ++i){for (int j=k+1; j<n; ++j){A[i][j] -= A[i][k] * A[k][j];}}}// 回代过程  for (int k=n-1; k>=0; --k){for (int j=k+1; j<n; ++j){b[k] -= A[k][j] * x[j];}x[k] = b[k] / A[k][k];}return x;
}

2 列主元高斯消去法

2.1 经典方法的致命问题

(∣A∣≠0)⇏(akk(k−1)≠0)(\vert A \vert \neq 0) \nRightarrow (a_{kk}^{(k-1)} \neq 0)(∣A∣=0)⇏(akk(k−1)=0)，即使akk(k−1)≠0a_{kk}^{(k-1)} \neq 0akk(k−1)=0，但∣akk(k−1)∣\vert a_{kk}^{(k-1)}\vert∣akk(k−1)∣很小，这就会引起其他元素数量级的剧增和舍入误差的增长，导致计算结果不可靠，甚至计算不能进行下去(上溢)。

2.2 按列选主元步骤的算法描述

相比经典的高斯消去法，在消元过程前额外多了按列选主元步骤和交换行步骤，其余流程是一致的。

按列选主元步骤的算法描述： 选取第kkk列的主元，即寻找它所在的行qqq，有：
q=arg max⁡i∈[k,n)∣aik∣q=\argmax\limits_{i\in [k, n)}\vert a_{ik}\vertq=i∈[k,n)argmax∣aik∣

交换行步骤的算法描述： 找到主元所在行qqq后，分别交换AAA和bbb的第kkk行和第qqq行元素，即
{swap(akj,aqj),j=k,k+1,...,nswap(bk,bq)\begin{cases} &\rm{swap}(a_{kj}, a_{qj}),j=k,k+1,...,n \\ & \rm{swap}(b_{k}, b_{q}) \end{cases} {swap(akj,aqj),j=k,k+1,...,nswap(bk,bq)

2.3 算法复杂度分析

时间复杂度：O(n3)O(n^3)O(n3)，增加的两个步骤的复杂度均为O(n)O(n)O(n)

空间复杂度：O(n2)O(n^2)O(n2)

2.4 算法优势

由于∣aik(k−1)∣/∣akk(k−1)∣≤1(i=k+1,...,n)\vert a_{ik}^{(k-1)}\vert/\vert a_{kk}^{(k-1)}\vert \le 1 \ (i=k+1, ..., n)∣aik(k−1)∣/∣akk(k−1)∣≤1 (i=k+1,...,n)，列主元消去法有利于控制误差的传播，故具有较好的数值稳定性；

2.5 算法的C++实现

typedef float elem;
vector<elem> Gauss_Elimination_ColumnPivot(vector<vector<elem>>& A, vector<elem>& b, elem epsilon){/* 列主元的高斯消去法参数简介：A - n维的系数矩阵b - 方程组等号右边的向量，维度1xnepsilon - 矩阵元素的阈值，当主元过小时高斯消去法将数值不稳定使用前的约定：1. 确保A不为空，且Ax=b中各个参数的维度一致注：执行完后，A和b将会改变*/int n = A.size();vector<elem> x(n);// 消元过程for (int k=0; k<n-1; ++k){ // 第k+1步消元elem maxf = 0;int q; // 最大主元所在行// 选主元for (int i=k; i<n; ++i){if(abs(maxf) < abs(A[i][k])){maxf = A[i][k];q = i;}}// 主元太小，不能继续if(abs(maxf) < epsilon){cout<<"The pivot element is too small. Process aborted."<<endl;return {};}// 交换行if(q != k){for (int j=k; j<n; ++j){swap(A[q][j], A[k][j]);}swap(b[q], b[k]);}// 消元过程for (int i=k+1; i<n; ++i){ // l[i][k]A[i][k] /= A[k][k];}// b[i]for (int i=k+1; i<n; ++i){b[i] -= A[i][k]*b[k];}// A[i][j]for (int i=k+1; i<n; ++i){for (int j=k+1; j<n; ++j){A[i][j] -= A[i][k] * A[k][j];}}}// 回代过程  for (int k=n-1; k>=0; --k){for (int j=k+1; j<n; ++j){b[k] -= A[k][j] * x[j];}x[k] = b[k] / A[k][k];}return x;
}

References

[1] 李乃成, 梅立泉. 数值分析[M]. 科学出版社, 2011.

笔者水平有限，如有错误欢迎批评指正~