Halo

A magic place for coding

0%

前言

在数值计算课程中学到了高斯消元法的详细算法过程，自己就用C++和matlab实现了一下。首先我会讲解算法的推导过程，然后我将会把代码分块讲解，想直接看完整代码的可以直接到最后看。

高斯消元法讲解

算法详解

$$\begin{cases} a_{11}x_1 + a_{12}x_2 + \cdots + a_{1n}x_n = b_1\ a_{21}x_1 + a_{22}x_2 + \cdots + a_{2n}x_n = b_2\ \vdots \ a_{n1}x_1 + a_{n2}x_2 + \cdots + a_{nn}x_n = b_n\ \end{cases} \tag{1.1}$$
，或写成矩阵形式
$$\begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n} \ a_{21} & a_{22} & \cdots & a_{2n} \ \vdots & \vdots & & \vdots & \ a_{n1} & a_{n2} & \cdots & a_{nn}\ \end{bmatrix} \begin{bmatrix} x_1\ x_2\ \vdots\ x_n \ \end{bmatrix} = \begin{bmatrix} b_1\ b_2\ \vdots\ b_n \ \end{bmatrix}$$
，简记为**Ax = b**.

(1) 第一步（k=1

$$m_{i1} = a_{i1}^{(1)} / a_{11}^{(1)}, i = 2,3,…,n.$$

$$\begin{bmatrix} a_{11}^{(1)} & a_{12}^{(1)} & \cdots & a_{1n}^{(1)} \ 0 & a_{22}^{(2)} & \cdots & a_{2n}^{(2)} \ \vdots & \vdots & & \vdots & \ 0 & a_{n2}^{(2)} & \cdots & a_{nn}^{(2)}\ \end{bmatrix} \begin{bmatrix} x_1\ x_2\ \vdots\ x_n \ \end{bmatrix} = \begin{bmatrix} b_1^{(1)}\ b_2^{(2)}\ \vdots\ b_n^{(2)} \ \end{bmatrix}$$

$$A^{(2)}x = b^{(2)}，$$

$$\begin{cases} a_{ij}^{(2)} = a_{ij}^{(1)} - m_{i1}a_{1j}^{(1)},i,j = 2,3,…,n.\ b_i^{(2)} = b_i ^{(1)} - m_{i1}b_1^{(1)}, i = 2,3,…,n.\ \end{cases}$$

(2) 第$k$次消元（$k = 1,2,…,n-1$）

$$\begin{bmatrix} a_{11}^{(1)} & a_{12}^{(1)} & \cdots & a_{1k}^{(1)} & \cdots &a_{1n}^{(1)} \ & a_{22}^{(2)} & \cdots & a_{2k}^{(2)} & \cdots & a_{2n}^{(2)} \ & & \ddots & \vdots & & \vdots & \ & & & a_{kk}^{(k)} & \cdots & a_{kn}^{(k)}\ & & & a_{nk}^{(k)} & \cdots & a_{nn}^{(k)} \end{bmatrix} \begin{bmatrix} x_1\ x_2\ \vdots\ x_k \ \vdots \ x_n \ \end{bmatrix} = \begin{bmatrix} b_1^{(1)}\ b_2^{(2)}\ \vdots\ b_k^{(k)}\ \vdots \ b_n^{(k)} \ \end{bmatrix},$$

$$m_{ik} = a_{ik}^{(k)} / a_{kk}^{(k)}, i = k+1,k+2,…,n.$$

$$\begin{cases} a_{ij}^{(k+1)} = a_{ij}^{(k)} - m_{ik}a_{kj}^{(k)},i,j = k+1,…,n.\ b_i^{(k+1)} = b_i ^{(k)} - m_{ik}b_k^{(k)}, i = k+1,…,n.\ \end{cases}，$$

（3）经过$n-1$步消元计算后，我们得到$A^{(n)}x = b^{(n)}$：
$$\begin{bmatrix} a_{11}^{(1)} & a_{12}^{(1)} & \cdots & a_{1n}^{(1)} \ & a_{22}^{(2)} & \cdots & a_{2n}^{(2)} \ & & \ddots & \vdots \ & & & a_{nn}^{(n)} \end{bmatrix} \begin{bmatrix} x_1\ x_2\ \vdots \ x_n \ \end{bmatrix} = \begin{bmatrix} b_1^{(1)}\ b_2^{(2)}\ \vdots \ b_n^{(k)} \ \end{bmatrix},$$

A是非奇异矩阵，且$a_{kk}^{(k)} ≠ 0（k=1,2,…,n-1）$，求解矩阵的回代公式为：
$$\begin{cases} x_n = b_n / a_{nn}^{(n)},\ x_k = b_k^{(k)} - \sum_{j=k+1}^n a_{ki}^{(k)}x_j / a_{kk}, k = n-1,n-2,…,1.\ \end{cases}，$$

代码分析

C++ Version

变量定义

高斯消元法主要是矩阵的运算，这里我们采用二维数组的形式存储矩阵，采用一维数组的形式存储向量。

顺序消元法

顺序消元法分为两个步骤：第一是消元，把原始矩阵变成一个上三角矩阵；第二是回代，通过最后一个方程一直回代解出整个矩阵。

列主元消元法

首先我要解释一下为什么会有这个方法。在顺序消元法中，我们看到了这样一个语句temp[i][k] = A[i][k] / A[k][k];，这是用来求出第i行与主元行的比值。注意到，如果某一个主元A[k][k]特别小的话，被一个相对大的数除，就会存在精度损失的问题。为了克服这个问题，我们希望在一列中，每一次都尽可能地挑选绝对值最大的一个元素作为主元，这就是列主元消元法
该方法的消元回代两个主要步骤和顺序消元法是一致的，不同的地方在于多了选主元换行这两步。

Matlab Version

Matlab实现矩阵运算非常简便，这里直接给出高斯消元法选主元的高斯消元法代码。

大致讨论一下运算复杂度，高斯消元法的计算复杂度大约是O($\frac{1}{3} n^3$)
高斯消元法的代码实现讲解到这里结束了，谢谢！

参考资料：

1.数值分析（第5版） 李庆扬，王能超，易大义

Welcome to my other publishing channels