求线性方程组解的直接方法

线性方程组具有一般形式：

$$ \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} $$

其中，$a_{ij}, b_i$为常数，$x_1, x_2, \cdots, x_n$为未知量。上式也可写成矩阵形式：

$$ \boldsymbol{A}x=b $$

如果线性方程组的系数矩阵$\boldsymbol{A}$的行列式不为$0$，即$\det \boldsymbol{A}\neq0$，线性方程组存在唯一的解。当$n$较小时，Cramer法则给出了线性方程组的公式解：

$$ x_i=\frac{\Delta_i}{\Delta},i=1,2\cdots $$

其中$\Delta$是$\boldsymbol{A}$的行列式，$\Delta_i$是指把$\boldsymbol{A}$的第$i$列换成$b$之后所得方阵的行列式，当$n$较大时，计算量巨大。线性方程组的解法可以分为直接解法和迭代解法。下面考虑使用直接解法，使用有限次四则运算，得到线性方程组的解。

一、消去法

1、高斯（Gauss）消去法

(1) 高斯消去法原理

为统一起见，将上面方程组改写成：

$$ \begin{cases} a_{11}^{(1)} x_1 + a_{12}^{(1)} x_2 + \cdots + a_{1n}^{(1)} x_n = b_1^{(1)} \\ a_{21}^{(1)} x_1 + a_{22}^{(1)} x_2 + \cdots + a_{2n}^{(1)} x_n = b_2^{(1)} \\ \quad\quad\quad\quad\quad\cdots\cdots \\ a_{n1}^{(1)} x_1 + a_{n2}^{(1)} x_2 + \cdots + a_{nn}^{(1)} x_n = b_n^{(1)} \end{cases} $$

右上角的数字用于区分消元前后的系数变换，简记为：

$$ \boldsymbol{A}^{(1)} \boldsymbol{x} = \boldsymbol{b}^{(1)} $$

其中：

$$ \boldsymbol{A}^{(1)} = \boldsymbol{A},\quad \boldsymbol{b}^{(1)} = \boldsymbol{b} $$

一般地，求解$n$阶方程组的高斯消去法步骤如下：

第一步：设 $a_{11}^{(1)}\neq0$，记 $l_{i1}=\frac{a_{i1}^{(1)}}{a_{11}^{(1)}} (i=2,3,\cdots,n)$，将方程组中第 $i$ 个方程减去第 $1$ 个方程乘以 $l_{i1}\ (i=2,3,\cdots,n)$，完成第一次消元，得同解方程组

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

其中：

$$ a_{ij}^{(2)}=a_{ij}^{(1)}-l_{i1}a_{1j}^{(1)},\qquad b_i^{(2)}=b_i^{(1)}-l_{i1}b_1^{(1)}\quad (i,j=2,3,\cdots,n) $$

方程组简记为

$$ \boldsymbol{A}^{(2)}\boldsymbol{x}=\boldsymbol{b}^{(2)} $$

第二步：设 $a_{22}^{(2)}\neq0$, 记 $l_{i2}=a_{i2}^{(2)}/a_{22}^{(2)}\ (i=3,\cdots,n)$。将方程组(2-4)中第 $i$ 个方程减去第 $2$ 个方程乘以 $l_{i2}\ (i=3,\cdots,n)$，完成第二次消元。

如此循环往复，直到完成$n-1$次消元过程，将原方程组化为上三角形方程组：

$$ \begin{cases} \begin{aligned} a_{11}^{(1)}x_1 + a_{12}^{(1)}x_2 + a_{13}^{(1)}x_3 + \cdots + a_{1n}^{(1)}x_n = b_1^{(1)}&\\ a_{22}^{(2)}x_2 + a_{23}^{(2)}x_3 + \cdots + a_{2n}^{(2)}x_n = b_2^{(2)}&\\ a_{33}^{(3)}x_3 + \cdots + a_{3n}^{(3)}x_n = b_3^{(3)}&\\ \vdots\qquad&\\ a_{nn}^{(n)}x_n = b_n^{(n)}& \end{aligned} \end{cases} $$

简记为：

$$ \boldsymbol{A}^{(n)}\boldsymbol{x}=\boldsymbol{b}^{(n)} $$

然后按照变量的逆序，逐步回代得到原来方程组的解。

(2) 高斯消去法代码

下面是这套算法的MATLAB代码：

function x = gauss_solve(A, b)
    n = length(b);
    for k = 1:n-1
        for i = k+1:n
            if abs(A(i,k)) > abs(A(k,k))
                A([k i], :) = A([i k], :);
                b([k i]) = b([i k]);
            end
        end
        if abs(A(k,k)) < eps
            x = [];
            return;
        end
        m = A(k+1:n, k) / A(k,k);
        A(k+1:n, k) = m;
        b(k+1:n) = b(k+1:n) - m * b(k);
        A(k+1:n, k+1:n) = A(k+1:n, k+1:n) - m * A(k, k+1:n);
    end
    if abs(A(n,n)) < eps
        x = [];
        return;
    end
    x = zeros(n,1);
    x(n) = b(n) / A(n,n);
    for i = n-1:-1:1
        x(i) = (b(i) - A(i, i+1:n) * x(i+1:n)) / A(i,i);
    end
end

(3) 高斯消去法计算量

由于计算机作乘除运算所需时间远大于作加减运算所需时间，故我们只讨论乘除运算量。对于三元一次线性方程组的情况，由消去法步骤知，在进行第 $k$ 次消元时，需作除法 $n-k$ 次，乘法 $(n-k)(n-k+1)$ 次，故消元过程中乘除运算总量为：

乘法次数：

$$ \sum_{k=1}^{n-1} (n-k)(n-k+1) = \frac{n}{3} (n^2-1) $$

除法次数：

$$ \sum_{k=1}^{n-1} (n-k) = \frac{n}{2} (n-1) $$

在回代过程中，计算 $x_k$ 需要 $n-k+1$ 次乘除法，整个回代过程需要乘除运算的总量为：

$$ \sum_{k=1}^{n} (n-k+1) = \frac{n}{2} (n+1) $$

所以，高斯消去法的乘除总运算量为：

$$ N = \frac{n}{3}(n^2-1) + \frac{n}{2}(n-1) + \frac{n}{2}(n+1)= \frac{n^3}{3} + n^2 - \frac{n}{3} $$

对于$n$阶线性方程组，运算次数的数量级为$n^3$。

2、主元素消去法

(1) 主元素消去法原理

高斯消去法的计算过程中要求$a_{kk}^{(k)}$（称为主元素）均不为$0$，否则高斯消去法不能进行。另一方面，如果主元素的绝对值非常小，将其作为除数用来计算$l_{ij}$，会使得舍入误差严重扩大，从而导致计算结果不正确。为了抑制舍入误差的增长，应该尽量避免小主元的出现，口蹄疫通过交换方程的次序，选取绝对值大的元素作为主元，基于这种想法导出主元素消去法。一般情况下，为了使程序设计简单，采用列主元素法。

为简便起见，使用增广矩阵

$$ [\boldsymbol{A}, \boldsymbol{b}] = \begin{bmatrix} a_{11} & \cdots & a_{1n} & b_1 \\ a_{21} & \cdots & a_{2n} & b_2 \\ \vdots & & \vdots & \vdots \\ a_{n1} & \cdots & a_{nn} & b_n \\ \end{bmatrix} $$

表示方程组，并直接在增广矩阵上运算。列主元素法的基本思想是在每次消元前，在要消去未知数的系数中找绝对值最大的系数作为主元，通过方程对换将其换到对角线上，然后进行消元。其具体步骤如下：

第一步：首先在矩阵的第$1$列中选取绝对值最大的元，比如 $a_{i1}$，则 $|a_{i1}| = \max\limits_{1 \leq k \leq n} |a_{k1}|$。将增广矩阵中第$i$行与第$1$行互换。为方便起见，记行互换后的增广矩阵为 $[A^{(1)}, b^{(1)}]$，然后进行第一次消元，得矩阵：

$$ [A^{(2)}, b^{(2)}] = \begin{bmatrix} a_{11}^{(2)} & a_{12}^{(2)} & \cdots & a_{1n}^{(2)} & b_1^{(2)} \\ 0 & a_{22}^{(2)} & \cdots & a_{2n}^{(2)} & b_2^{(2)} \\ \vdots & \vdots & & \vdots & \vdots \\ 0 & a_{n2}^{(2)} & \cdots & a_{nn}^{(2)} & b_n^{(2)} \\ \end{bmatrix} $$

第二步：在矩阵 $[A^{(2)}, b^{(2)}]$ 的第$2$列中选主元，比如 $a_{j2}^{(2)}$，使：

$$ |a_{j2}^{(2)}| = \max\limits_{2 \leq k \leq n} |a_{k2}^{(2)}| $$

将矩阵 $[A^{(2)}, b^{(2)}]$ 的第$j$行与第$2$行互换，再进行第二次消元，得矩阵 $[A^{(3)}, b^{(3)}]$。

第三步：在矩阵 $[A^{(k)}, b^{(k)}]$ 的第$k$列中选主元，如 $a_{ik}^{(k)}$ 使：

$$ |a_{ik}^{(k)}| = \max\limits_{k \leq s \leq n} |a_{sk}^{(k)}| $$

将 $[A^{(k)}, b^{(k)}]$ 的第$i$行与第$k$行互换，进行第$k$次消元。如此经过 $n-1$ 步，增广矩阵被化成上三角形矩阵，最后由回代过程求解。

(2) 列主元素消去法代码

下面是这套算法的MATLAB代码：

function x = gauss_solve(A, b)
    n = length(b);
    for k = 1:n-1
        [~, p] = max(abs(A(k:n, k)));        % 寻找最大主元所在行数
        p = p + k - 1;
        if p ~= k                            % 交换
            A([k p], :) = A([p k], :);
            b([k p]) = b([p k]);
        end
        tol = eps * norm(A, inf);            % 防止因为数值稳定性问题对是否有解的判定造成影响
        if abs(A(k,k)) < tol
            x = [];
            return;
        end
        m = A(k+1:n, k) / A(k,k);
        A(k+1:n, k) = m;
        b(k+1:n) = b(k+1:n) - m * b(k);
        A(k+1:n, k+1:n) = A(k+1:n, k+1:n) - m * A(k, k+1:n);
    end
    if abs(A(n,n)) < tol
        x = [];
        return;
    end
    x = zeros(n,1);
    x(n) = b(n) / A(n,n);
    for i = n-1:-1:1
        x(i) = (b(i) - A(i, i+1:n) * x(i+1:n)) / A(i,i);
    end
end

相比于原始的高斯消去法，这里使用用矩阵向量运算代替 for 循环，可以显著提高程序运行速度。

二、直接三角分解法

1、高斯消去法的矩阵形式（LU分解）

(1) LU分解原理

高斯消去法的消元过程可以使用一串初等矩阵左乘增广矩阵表示，例如将第一行同乘$-l_{31}$后与第三行相加，以消去第三行第一列位置的系数，可以使用下面的初等行变换矩阵左乘来表达：

$$ \begin{bmatrix} 1 & 0 & 0 & \cdots & 0 \\ 0 & 1 & 0 & \cdots & 0 \\ -l_{31} & 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & 1 \end{bmatrix} $$

那么对于同一列的消元，可以将上面类型的初等行变换相继左乘得到，例如第一次消元等价于使用矩阵$\boldsymbol{L}_1$左乘增广矩阵：

$$ \boldsymbol{L}_1 = \begin{bmatrix} 1 & 0 & 0 & \cdots & 0 \\ -l_{21} & 1 & 0 & \cdots & 0 \\ -l_{31} & 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ -l_{n1} & 0 & 0 & \cdots & 1 \end{bmatrix} $$

那么第$k$次消元可以使用下面增广矩阵$\boldsymbol{L}_k$左乘得到：

$$ \boldsymbol{L}_k = \begin{bmatrix} 1 & & & & & \\ & \ddots & & & & \\ & & 1 & & & \\ & & -l_{k+1,\,k} & 1 & & \\ & & \vdots & & \ddots & \\ & & -l_{n k} & & & 1 \end{bmatrix} $$

经过$n-1$次消元后，得到：

$$ \begin{aligned} \left[ \boldsymbol{A}^{(n)},\, \boldsymbol{b}^{(n)} \right] &= \begin{bmatrix} a_{11}^{(1)} & a_{12}^{(1)} & \cdots & a_{1n}^{(1)} & b_1^{(1)} \\ & a_{22}^{(2)} & \cdots & a_{2n}^{(2)} & b_2^{(2)} \\ & & \ddots & \vdots & \vdots \\ & & & a_{nn}^{(n)} & b_n^{(n)} \end{bmatrix} \\ &= \boldsymbol{L}_{n-1} \left[ \begin{aligned} \boldsymbol{A}^{(n-1)},\, \boldsymbol{b}^{(n-1)} \end{aligned} \right] \\ &= \boldsymbol{L}_{n-1} \boldsymbol{L}_{n-2} \cdots \boldsymbol{L}_1 \left[ \begin{aligned} \boldsymbol{A}^{(1)},\, \boldsymbol{b}^{(1)} \end{aligned} \right] \end{aligned} $$

因为$\boldsymbol{L}_k$为非奇异矩阵，存在逆矩阵，可以得到：

$$ \boldsymbol{L}^{-1}_k = \begin{bmatrix} 1 & & & & & \\ & \ddots & & & & \\ & & 1 & & & \\ & & l_{k+1,\,k} & 1 & & \\ & & \vdots & & \ddots & \\ & & l_{n k} & & & 1 \end{bmatrix} $$

于是有：

$$ \begin{aligned} \left[ \boldsymbol{A}, \boldsymbol{b} \right] &= \left[ \boldsymbol{A}^{(1)}, \boldsymbol{b}^{(1)} \right] \\ &= \boldsymbol{L}_1^{-1} \boldsymbol{L}_2^{-1} \cdots \boldsymbol{L}_{n-1}^{-1} \left[ \boldsymbol{A}^{(n)},\; \boldsymbol{b}^{(n)} \right] \end{aligned} $$

令：

$$ \boldsymbol{L} = \boldsymbol{L}_1^{-1} \boldsymbol{L}_2^{-1} \cdots \boldsymbol{L}_{n-1}^{-1} $$

$$ \boldsymbol{U}=\boldsymbol{A}^{(n)} $$

可以得到：

$$ \left[\boldsymbol{A},\, \boldsymbol{b}\right] = \left[ \boldsymbol{L} \boldsymbol{U},\, \boldsymbol{L} \boldsymbol{b}^{(n)} \right] \tag{2-14} $$

则有：

$$ \boldsymbol{A}=\boldsymbol{LU}\Rightarrow\boldsymbol{U}=\boldsymbol{L}^{-1}\boldsymbol{A} $$

这说明，消元过程实际上是把系数矩阵 $\boldsymbol{A}$ 分解成单位下三角形矩阵与上三角形矩阵的乘积的过程，其中 $\boldsymbol{L}$ 为单位下三角形矩阵，$\boldsymbol{U}$ 为上三角形矩阵：

$$ \boldsymbol{L} = \begin{bmatrix} 1 & & & & \\ l_{21} & 1 & & & \\ l_{31} & l_{32} & 1 & & \\ \vdots & \vdots & \vdots & \ddots & \\ l_{n1} & l_{n2} & \cdots & l_{n,n-1} & 1 \end{bmatrix} \qquad \boldsymbol{U} = \begin{bmatrix} u_{11} & u_{12} & u_{13} & \cdots & u_{1n} \\ & u_{22} & u_{23} & \cdots & u_{2n} \\ & & \ddots & \vdots & \vdots \\ & & & u_{n-1,n-1} & u_{n-1,n} \\ & & & & u_{nn} \end{bmatrix} $$

上述分解称为杜利特尔（Doolittle）分解，也称为LU分解。当系数矩阵进行三角分解后，求解方程组 $\boldsymbol{A} \boldsymbol{x} = \boldsymbol{b}$ 的问题就变得十分容易（尤其是存在多组$\boldsymbol{b}$的情况），它等价于求解两个三角形方程组 $\boldsymbol{L} \boldsymbol{y} = \boldsymbol{b}$ 和 $\boldsymbol{U} \boldsymbol{x} = \boldsymbol{y}$。因此，线性方程组问题可转化为矩阵的三角分解问题。

高斯消去法要求主元元素 $a_{kk}^{(k)} \neq 0\ (k=1, \cdots, n-1)$，将这一条件与矩阵 $\boldsymbol{A}$ 自身的性质联系起来，可以得到定理如下：

设 $\boldsymbol{A}$ 为 $n$ 阶方阵，若 $\boldsymbol{A}$ 的顺序主子式 $A_i (i=1,2,\cdots, n-1)$ 均不为零，则矩阵 $\boldsymbol{A}$ 存在唯一的LU分解。

(2) LU分解代码

下面是计算得出矩阵$\boldsymbol{L},\boldsymbol{U}$的MATLAB代码：

function [L, U] = lu_decomp(A)
    n = size(A,2);
    U = A;
    L = eye(n);
    for k = 1:n-1
        L(k+1:n, k) = U(k+1:n, k) / U(k, k);
        U(k+1:n, k+1:n) = U(k+1:n, k+1:n) - L(k+1:n, k) * U(k, k+1:n);
        U(k+1:n, k) = 0;
    end
end

下面是通过矩阵$\boldsymbol{L},\boldsymbol{U}$求解方程组的MATLAB代码：

function x = LU_solve(L, U, b)
    n = size(b, 1);
    y = zeros(n, 1);
    x = zeros(n, 1);
    for i = 1:n
        y(i) = b(i) - L(i, 1:i-1) * y(1:i-1);
    end
    for i = n:-1:1
        x(i) = (y(i) - U(i, i+1:n) * x(i+1:n)) / U(i, i);
    end
end

2、列主元素的三角分解法

从直接三角分解公式知，当 $u_{jj}=0\ (j=1,\cdots,n-1)$ 时，计算将中断；或当 $|u_{jj}|\, (j=1,\cdots,n-1)$ 很小时，按公式计算会使舍入误差扩散。但若矩阵 $\boldsymbol{A}$ 非奇异，可增加列选主元过程，即先进行一系列的行交换，保证LU分解能够实现。从而有如下定理：

设矩阵 $\boldsymbol{A}$ 非奇异, 则存在置换矩阵 $\boldsymbol{P}$，使得 $\boldsymbol{PA}$ 有唯一的 Doolittle 分解，其中 $\boldsymbol{L}$ 是单位下三角形矩阵，$\boldsymbol{U}$ 是上三角形矩阵，且 $|l_{ij}| \leq 1$。即存在置换矩阵$\boldsymbol{P}$、单位下三角矩阵$\boldsymbol{L}$和上三角矩阵$\boldsymbol{U}$，使得$\boldsymbol{PA}=\boldsymbol{LU}$。

我们一般不使用这种方法，但是下面给出一份代码，在已知$\boldsymbol{P},\boldsymbol{L},\boldsymbol{U}$的情况下求解线性方程组：

function x = solveWithPLU(P, L, U, b)
    n = size(b, 1);
    y = zeros(n, 1);
    x = zeros(n, 1);
    Pb = P * b;
    for i = 1:n
        y(i) = Pb(i) - L(i, 1:i-1) * y(1:i-1);
    end
    for i = n:-1:1
        x(i) = (y(i) - U(i, i+1:n) * x(i+1:n)) / U(i, i);
    end
end

三、特殊矩阵的三角分解法

1、解三对角方程组的追赶法

(1) 原理

方程组：

$$ \boldsymbol{A}= \begin{bmatrix} b_1 & c_1 & & & & \\ a_2 & b_2 & c_2 & & & \\ & \ddots & \ddots & \ddots & & \\ & & a_{n-1} & b_{n-1} & c_{n-1} & \\ & & & a_n & b_n & \end{bmatrix},\qquad \boldsymbol{x}= \begin{bmatrix} x_1\\[4pt] x_2\\[4pt] \vdots\\[4pt] x_{n-1}\\[4pt] x_n \end{bmatrix},\qquad \boldsymbol{d}= \begin{bmatrix} d_1\\[4pt] d_2\\[4pt] \vdots\\[4pt] d_{n-1}\\[4pt] d_n \end{bmatrix} $$

$$ \boldsymbol{A}\,\boldsymbol{x}=\boldsymbol{d} $$

称为三对角方程组，其中系数矩阵$\boldsymbol{A}$是一种特殊矩阵，其非零元素集中分布在主对角线及相邻的两条次对角线上，为三对角矩阵。设方程组的系数矩阵满足：

$$ \begin{cases} |b_1| > |c_1| > 0\\ |b_i| \ge |a_i| + |c_i| \\[4pt] |b_n| > |a_n| > 0 \end{cases} \qquad (a_i,c_i \ne 0,\ \text{且}\ i=2,3,\cdots,n-1) $$

则 $\boldsymbol{A}$ 可唯一分解为：

$$ \boldsymbol{A}=\boldsymbol{L}\boldsymbol{U}= \begin{bmatrix} 1 & & & & \\ l_2 & 1 & & & \\ & l_3 & 1 & & \\ & & \ddots & \ddots & \\ & & & l_n & 1 \end{bmatrix} \begin{bmatrix} u_1 & c_1 & & & & \\ & u_2 & c_2 & & & \\ & & u_3 & c_3 & & \\ & & & \ddots & \ddots & \\ & & & & \ddots & c_{n-1} \\ & & && &u_n \end{bmatrix} $$

其中 $c_i\ (i=1,2,\cdots,n-1)$ 由原矩阵给出。

[!NOTE]

条件的引出： 这三个条件在数学上被称为严格对角占优（首尾严格，中间非严格但整体不可约）。 在后面的计算中，我们需要不断地除以主对角线相关的元素（即公式中的$u_i$）。如果满足上述对角占优条件，可以从理论上绝对保证在计算过程中所有的主元$u_i\neq0$。这意味着算法不会因为除以零而崩溃，保证了分解的唯一性和计算过程的数值稳定性。

利用待定系数法计算，可以得到：

$$ \begin{cases} u_1 = b_1,\\ l_i = \dfrac{a_i }{ u_{i-1}} \\ u_i = b_i - c_{i-1} l_i \end{cases} \qquad (i=2,3,\cdots,n) $$

当矩阵 $\boldsymbol{A}$ 三角分解后, 解线性方程组化为求解方程组 $\boldsymbol{Ly}=\boldsymbol{d}$ 和 $\boldsymbol{Ux}=\boldsymbol{y}$。

解 $\boldsymbol{Ly}=\boldsymbol{d}$，得：

$$ \begin{cases} y_1 = d_1,\\ y_k = d_k - l_k y_{k-1} \end{cases} \qquad (k=2,3,\cdots,n) $$

再解 $\boldsymbol{Ux}=\boldsymbol{y}$ ，得三对角方程组的解：

$$ \begin{cases} x_n = \dfrac{y_n}{u_n}\\ x_k = \dfrac{y_k - c_k x_{k+1}}{u_k} \end{cases} \qquad (k=n-1,n-2,\cdots,1) $$

称上述求解三对角方程组的方法为追赶法。由于矩阵中出现了大量零元素，从而使得计算公式简化，计算量大为减少，乘除法的运算量为$5n-4$。

(2) 代码

由于矩阵的大部分位置都是 $0$，因此系统不会给你整个矩阵，而是只给您传入一个长为$n$的向量$\boldsymbol{a}$ ，以及长为$n-1$的向量$\boldsymbol{p},\boldsymbol{q}$。分别代表上文中的数组$\boldsymbol{b},\boldsymbol{a},\boldsymbol{c}$，即：

$$ A = \begin{bmatrix} a_1 & p_1 & 0 & \cdots & 0 & 0 \\ q_1 & a_2 & p_2 & \cdots & 0 & 0 \\ 0 & q_2 & a_3 & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & a_{n-1} & p_{n-1} \\ 0 & 0 & 0 & \cdots & q_{n-1} & a_n \end{bmatrix} $$

实现求解的MATLAB代码如下：

function x = solve_tridiag(a, p, q, b)
    n = length(a);
    l = zeros(n, 1);
    u = zeros(n, 1);
    u(1) = a(1);
    for i = 1:n
        if i > 1
            l(i) = q(i-1) / u(i-1);
            u(i) = a(i) - l(i) * p(i-1);
        end
    end
    y = zeros(n, 1);
    y(1) = b(1);
    for i = 2:n
        y(i) = b(i) - l(i) * y(i-1);
    end
    x = zeros(n, 1);
    x(n) = y(n) / u(n);
    for i = n-1:-1:1
        x(i) = (y(i) - p(i) * x(i+1)) / u(i);
    end
end

2、平方根法

(1) 朴素平方根法（Cholesky分解法）

[!NOTE]

前置定理：设$\boldsymbol{A}$是对称正定矩阵，则存在唯一的非奇异下三角形矩阵$\boldsymbol{L}$，使得：

$$ \boldsymbol{A}=\boldsymbol{LL}^\mathrm{T} $$

且$\boldsymbol{L}$的对角元素皆为正数。矩阵的这种分解叫做Cholesky分解。

根据前置定理，使用待定系数法可以导出$\boldsymbol{L}$的计算公式，设：

$$ \boldsymbol{L}= \begin{bmatrix} l_{11} & & & \\ l_{21} & l_{22} & & \\ \vdots & \vdots & \ddots & \\ l_{n1} & l_{n2} & \cdots & l_{nn} \end{bmatrix} $$

比较 $\boldsymbol{A}$ 与 $\boldsymbol{L}\boldsymbol{L}^{\mathrm{T}}$ 的相应元素，可得：

$$ \begin{cases} l_{kk} = \left(a_{kk} - \displaystyle\sum_{j=1}^{k-1} l_{kj}^2\right)^{1/2} & k=1,2,\dots,n\\ l_{ik} = \dfrac{a_{ik} - \displaystyle\sum_{j=1}^{k-1} l_{ij}l_{kj}}{l_{kk}} & i=k+1,\dots,n \end{cases} $$

计算顺序是按列进行，即：

$$ l_{11}\to l_{i1}\ (i=2,3,\dots,n)\to l_{22}\to l_{i2}\ (i=3,\dots,n)\to\cdots $$

当矩阵 $\boldsymbol{A}$ 完成 Cholesky 分解后，求解方程组 $\boldsymbol{Ax}=\boldsymbol{b}$ 就转化为依次求解方程组：

$$ \boldsymbol{Ly}=\boldsymbol{b},\qquad \boldsymbol{L}^{\mathrm{T}}\boldsymbol{x}=\boldsymbol{y} $$

它们的解分别为：

$$ y_k=\dfrac{b_k-\displaystyle\sum_{j=1}^{k-1} l_{kj}y_j}{l_{kk}}\qquad k=1,2,\dots,n $$

$$ x_k=\dfrac{y_k-\displaystyle\sum_{j=k+1}^{n} l_{jk}x_j}{l_{kk}}\qquad k=n,n-1,\dots,1 $$

求解线性方程组的上述方法称为平方根法，也称为 Cholesky 分解法。这种方法无需选主元，计算过程也是稳定的。由于 $\boldsymbol{A}$ 的对称性，平方根法的乘除运算量为约 $n^3/6$ 数量级，约是 Gauss 消去法的一半。但这种方法在求 $\boldsymbol{L}$ 时需作 $n$ 次开方运算，这样又增加了计算量。

(2) 改进平方根法（$\boldsymbol{LDL}^\mathrm{T}$法）

上面朴素方法中的前置定理的证明过程表明，对称正定矩阵$\boldsymbol{A}$又可以作出如下分解：

$$ \boldsymbol{A}=\boldsymbol{LDL}^\mathrm{T} $$

其中$\boldsymbol{L}$为单位下三角形矩阵，$\boldsymbol{D}$为对角矩阵。记：

$$ \boldsymbol{L}= \begin{bmatrix} 1 & & & \\ l_{21} & 1 & & \\ \vdots & \vdots & \ddots & \\ l_{n1} & l_{n2} & \cdots & 1 \end{bmatrix} $$

$$ \boldsymbol{D}=\mathrm{diag}(d_1,\cdots,d_n) $$

然后利用待定系数法，导出$\boldsymbol{LDL}^\mathrm{T}$分解的计算公式，对$k=1,2,\cdots,n$：

$$ \begin{cases} \begin{aligned} d_k &= a_{kk} - \sum_{j=1}^{k-1} l_{kj}^2\,d_j\\ l_{ik} &= \frac{a_{ik} - \sum_{j=1}^{k-1} l_{ij}\,d_j\,l_{kj}}{d_k}\quad (i=k+1,\dots,n) \end{aligned} \end{cases} $$

计算顺序如下：

$$ d_1\to l_{i1}\ (i=2,3,\dots,n)\to d_2\to l_{i2}\ (i=3,\dots,n)\to d_3\to\cdots $$

为避免重复计算，引进辅助量：

$$ u_{ik}=l_{ik}d_k\quad (k=1,2,\dots,n,\; i=k+1,\dots,n) $$

则上式可改写成：

$$ \begin{cases} \begin{aligned} d_k &= a_{kk} - \sum_{j=1}^{k-1} u_{kj}\,l_{kj}\quad (k=1,2,\dots,n)\\ u_{ik} &= a_{ik} - \sum_{j=1}^{k-1} u_{ij}\,l_{kj}\quad (i=k+1,\dots,n)\\ l_{ik} &= \frac{u_{ik}}{d_k} \end{aligned} \end{cases} $$

按上式进行 $\boldsymbol{L}\boldsymbol{D}\boldsymbol{L}^{\mathrm{T}}$ 分解，乘除运算量与 Cholesky 分解相当，且避免了开方运算。

矩阵 $\boldsymbol{A}$ 作 $\boldsymbol{L}\boldsymbol{D}\boldsymbol{L}^{\mathrm{T}}$ 分解后，解方程组 $\boldsymbol{Ax}=\boldsymbol{b}$ 可分两步进行：先解方程组 $\boldsymbol{Ly}=\boldsymbol{b}$，再由 $\boldsymbol{L}^\mathrm{T}\boldsymbol{x}=\boldsymbol{D}^{-1}\boldsymbol{y}$ 求 $\boldsymbol{x}$。具体计算公式为

$$ \begin{cases} \begin{aligned} y_1 &= b_1\\ y_k &= b_k - \sum_{j=1}^{k-1} l_{kj}\,y_j\quad (k=2,\dots,n)\\ x_n &= \frac{y_n}{d_n}\\ x_k &= \frac{y_k}{d_k} - \sum_{j=k+1}^{n} l_{jk}\,x_j\quad (k=n-1,\dots,1) \end{aligned} \end{cases} $$

求解线性方程组的这一方法称为改进平方根法，也叫 $\boldsymbol{L}\boldsymbol{D}\boldsymbol{L}^{\mathrm{T}}$ 法。

[!TIP]

能够进行$\boldsymbol{L}\boldsymbol{L}^{\mathrm{T}}$分解的矩阵需要是方阵，并且矩阵是正定的，如果是半正定矩阵，分解仍然存在，但是对角线上会出现$0$，且分解不唯一，不能用用求解；$\boldsymbol{L}\boldsymbol{D}\boldsymbol{L}^{\mathrm{T}}$分解时，矩阵要求满足$\boldsymbol{A}=\boldsymbol{A}^\mathrm{T}$，且各阶顺序主子式不为零，但是矩阵不必正定，其根本原因在与缩放因子储存在对角矩阵$\boldsymbol{D}$中，计算过程不涉及开方运算。

(3) 代码

下面是对矩阵做$\boldsymbol{L}\boldsymbol{L}^{\mathrm{T}}$分解的MATLAB代码：

function L = decompLL(A)
    L = eye(size(A));
    n = size(A, 1);
    d = zeros(n, 1);
    for k = 1:n
        L(k, k) = sqrt(A(k, k) - L(k, 1:k-1) * (L(k, 1:k-1))');
        L(k+1:n, k) = (A(k+1:n, k) - L(k+1:n, 1:k-1) * (L(k, 1:k-1))') / L(k, k);
    end
end

下面是利用$\boldsymbol{L}\boldsymbol{L}^{\mathrm{T}}$法求解的MATLAB代码：

function x = solveWithLL(L, b)
    n = length(b);
    y = zeros(1, n);
    x = zeros(1, n);
    for k = 1:n
        y(k) = b(k) - L(k, 1:k-1) * (y(1:k-1))';
    end
    for k = n:-1:1
        x(k) = (y(k) - x(k+1:n) * L(k+1:n, k)) / L(k, k);
    end
end

下面是对矩阵做$\boldsymbol{L}\boldsymbol{D}\boldsymbol{L}^{\mathrm{T}}$分解的MATLAB代码：

function [L, D] = decompLDL(A)
    D = zeros(size(A));
    L = eye(size(A));
    U = zeros(size(A));
    n = size(A, 1);
    d = zeros(n, 1);
    for k = 1:n
        d(k) = A(k, k) - U(k, 1:k-1) * (L(k, 1:k-1))';
        D(k, k) = d(k);
        U(k+1:n, k) = A(k+1:n, k) - U(k+1:n, 1:k-1) * (L(k, 1:k-1))';
        L(k+1:n, k) = U(k+1:n, k) / d(k);
    end
end

下面是利用$\boldsymbol{L}\boldsymbol{D}\boldsymbol{L}^{\mathrm{T}}$法求解的MATLAB代码：

function x = solveWithLDL(L, D, b)
    n = length(b);
    y = zeros(1, n);
    x = zeros(1, n);
    y(1) = b(1);
    for k = 2:n
        y(k) = b(k) - L(k, 1:k-1) * (y(1:k-1))';
    end
    x(n) = y(n) / D(n, n);
    for k = n-1:-1:1
        x(k) = y(k) / D(k, k) - x(k+1:n) * L(k+1:n, k);
    end
end

四、误差分析

1、向量范数

(1) 定义

设$\boldsymbol{x}$为 $\mathbb{R}^n$ 中的向量。函数 $\|\cdot\|:\mathbb{R}^n\to\mathbb{R}$ 称为一个范数，当且仅当对任意 $\boldsymbol{x},\boldsymbol{y}\in\mathbb{R}^n$ 和任意标量 $\alpha\in\mathbb{R}$，满足：

1. 正定性（非负性且仅在零向量时为零）：

   $$ \|\boldsymbol{x}\|\ge 0,\qquad \|\boldsymbol{x}\|=0 \iff \boldsymbol{x}=\boldsymbol{0} $$

2. 齐次性（绝对齐次）：

   $$ \|\alpha\boldsymbol{x}\|=|\alpha|\|\boldsymbol{x}\| $$

3. 三角不等式：

   $$ \|\boldsymbol{x}+\boldsymbol{y}\|\le \|\boldsymbol{x}\|+\|\boldsymbol{y}\| $$

这三条构成范数的公理。

(2) 常见向量范数（p-范数及特殊范数）

对 $1\le p<\infty$，定义向量 $\boldsymbol{x}=(x_1,\dots,x_n)^\top$ 的 $p$-范数为

$$ \|\boldsymbol{x}\|_p=\Big(\sum_{i=1}^n |x_i|^p\Big)^{1/p} $$

三个常用的情形：

• 1-范数（曼哈顿范数）：

  $$ \|\boldsymbol{x}\|_1=\sum_{i=1}^n |x_i| $$

• 2-范数（欧几里得范数）：

  $$ \|\boldsymbol{x}\|_2=\Big(\sum_{i=1}^n x_i^2\Big)^{1/2} $$

• $\infty$-范数（切比雪夫范数）：

  $$ \|\boldsymbol{x}\|_\infty=\max_{1\le i\le n}|x_i| $$

这些都是满足范数公理的函数（下文给出证明与必要不等式）。

(3) 有限维空间范数之间的等价性

在有限维向量空间 $\mathbb{R}^n$ 中，存在常数 $c,C>0$ 使得对所有 $\boldsymbol{x}$ 都有：

$$ c\|\boldsymbol{x}\|_a \le \|\boldsymbol{x}\|_b \le C\|\boldsymbol{x}\|_a $$

我们称任意两范数都是等价的。可以证明，1-范数、2-范数和$\infty$-范数是等价的。

2、矩阵范数

(1) 定义

设 $\boldsymbol{A}\in\mathbb{R}^{n\times n}$（此处只研究方阵）。矩阵范数 $\|\cdot\|$ 是定义在矩阵集合上的函数，若对任意方阵 $\boldsymbol{A},\boldsymbol{B}$ 和标量 $\alpha$ 满足：

1. $\|\boldsymbol{A}\|\ge 0$, 且 $\|\boldsymbol{A}\|=0 \iff \boldsymbol{A}=\boldsymbol{0}$
2. $\|\alpha\boldsymbol{A}\|=|\alpha|\,\|\boldsymbol{A}\|$
3. $\|\boldsymbol{A}+\boldsymbol{B}\|\le \|\boldsymbol{A}\|+\|\boldsymbol{B}\|$
4. $\|\boldsymbol{A}\boldsymbol{B}\|\le \|\boldsymbol{A}\|\|\boldsymbol{B}\|$（相容性条件）

上述1-3条为向量范数的直接推广，第四条会使矩阵范数在数值计算中使用更加方便。此外，常用的一类矩阵范数是由向量范数诱导的范数，定义如下：

给定向量范数 $\|\cdot\|$（域为 $\mathbb{R}^n$ ），对矩阵 $\boldsymbol{A}\in\mathbb{R}^{n\times n}$ 定义诱导范数（也称算子范数）：

$$ \|\boldsymbol{A}\|=\max_{\boldsymbol{x}\neq 0}\frac{\|\boldsymbol{A}\boldsymbol{x}\|}{\|\boldsymbol{x}\|} $$

其满足矩阵范数的性质要求。

(2) 常见矩阵向量

• 1-范数（矩阵）：

  $$ \|\boldsymbol{A}\|_1=\max_{1\le j\le n}\sum_{i=1}^n |a_{ij}| $$

  即各列绝对值和的最大值（最大列和）。

• $\infty$-范数（矩阵）：

  $$ \|\boldsymbol{A}\|_\infty=\max_{1\le i\le n}\sum_{j=1}^n |a_{ij}|, $$

  即各行绝对值和的最大值（最大行和）。

• 2-范数（谱范数）：

  $$ \|\boldsymbol{A}\|_2=\max_{\|\boldsymbol{x}\|_2=1}\|\boldsymbol{A}\boldsymbol{x}\|_2=\sqrt{\lambda_{\max}(\boldsymbol{A}^\mathrm{T}\boldsymbol{A})} $$

  这里 $\lambda_{\max}(\boldsymbol{A}^\mathrm{T}\boldsymbol{A})$ 表示对称正定矩阵 $\boldsymbol{A}^\mathrm{T}\boldsymbol{A}$ 的最大特征值。

• Frobenius 范数（F-范数，希尔伯特–施密特范数）：

  $$ \|\boldsymbol{A}\|_F=\sqrt{\sum_{i=1}^n\sum_{j=1}^n a_{ij}^2} $$

  矩阵的F-范数可以认为是向量2-范数的直接推广，将$n$阶方阵看做$n^2$维向量所得到。可以证明，矩阵的F-范数与向量的2-范数相容，但F-范数不是由向量范数诱导出的矩阵范数，因为$\|\boldsymbol{I}\|=\sqrt{n}$，而诱导出的矩阵范数上述值都为$1$。

(3) 等价性与误差

如果将矩阵范数看作 $\mathbb{R}^{n^2}$ 空间上的向量范数，则由向量范数的等价性可得矩阵范数的等价性。

矩阵的误差可用矩阵范数表示。设 $\boldsymbol{A}^*$ 是 $\boldsymbol{A}$ 的近似矩阵，$\lVert \boldsymbol{A}-\boldsymbol{A}^*\rVert$、$\dfrac{\lVert \boldsymbol{A}-\boldsymbol{A}^*\rVert}{\lVert \boldsymbol{A}\rVert}$ 分别称为 $\boldsymbol{A}^*$ 的关于范数 $\lVert\cdot\rVert$ 的绝对误差与相对误差。

3、方程组的状态与条件数

一个实际问题化为数学问题，初始数据往往会有误差，即有扰动，从而使计算结果产生误差，因此需要研究扰动对解的影响。

当一个方程组，由于系数矩阵或右端项的微小扰动，而引起解发生巨大变化时，称该方程组是“病态”的。为了定量刻画方程组“病态”的程度，下面对方程组：

$$ \boldsymbol{A}\boldsymbol{x}=\boldsymbol{b} $$

就系数矩阵或右端项分别有扰动的两种情形进行讨论。首先考察右端项 $\boldsymbol{b}$ 的扰动对解的影响。设 $\boldsymbol{b}$ 有扰动 $\delta\boldsymbol{b}$，相应的解 $\boldsymbol{x}$ 的扰动记为 $\delta\boldsymbol{x}$，即：

$$ \boldsymbol{A}(\boldsymbol{x}+\delta\boldsymbol{x})=\boldsymbol{b}+\delta\boldsymbol{b} $$

由 $\boldsymbol{A}\boldsymbol{x}=\boldsymbol{b}$，得 $\boldsymbol{A}\,\delta\boldsymbol{x}=\delta\boldsymbol{b}$，从而有：

$$ \delta\boldsymbol{x}=\boldsymbol{A}^{-1}\delta\boldsymbol{b} $$

两边取范数，得：

$$ \|\delta\boldsymbol{x}\|=\|\boldsymbol{A}^{-1}\delta\boldsymbol{b}\| \le \|\boldsymbol{A}^{-1}\|\,\|\delta\boldsymbol{b}\| $$

又因为对于任意向量 $\boldsymbol{x}$ 有算子范数的性质：

$$ \|\boldsymbol{A}\boldsymbol{x}\|\le\|\boldsymbol{A}\|\,\|\boldsymbol{x}\| $$

从而：

$$ \|\boldsymbol{x}\|\ge\frac{\|\boldsymbol{A}\boldsymbol{x}\|}{\|\boldsymbol{A}\|} =\frac{\|\boldsymbol{b}\|}{\|\boldsymbol{A}\|} $$

这里用到 $\boldsymbol{A}\boldsymbol{x}=\boldsymbol{b}$。所以将上面两式相除，得到：

$$ \frac{\|\delta\boldsymbol{x}\|}{\|\boldsymbol{x}\|} \le\frac{\|\boldsymbol{A}^{-1}\|\,\|\delta\boldsymbol{b}\|} {\dfrac{\|\boldsymbol{b}\|}{\|\boldsymbol{A}\|}} =\|\boldsymbol{A}\|\,\|\boldsymbol{A}^{-1}\|\,\frac{\|\delta\boldsymbol{b}\|}{\|\boldsymbol{b}\|} $$

上式表明，当右端项有扰动时，解的相对误差不超过右端项的相对误差的 $\|\boldsymbol{A}\|\,\|\boldsymbol{A}^{-1}\|$ 倍。

如果 $\delta\boldsymbol{A}$ 充分小，使得 $\|\boldsymbol{A}^{-1}\|\,\|\delta\boldsymbol{A}\|<1$，则由上式得：

首先由扰动方程出发：

$$ (\boldsymbol{A}+\delta\boldsymbol{A})(\boldsymbol{x}+\delta\boldsymbol{x})=\boldsymbol{b}+\delta\boldsymbol{b} $$

展开得：

$$ \boldsymbol{A}\,\delta\boldsymbol{x}+\delta\boldsymbol{A}\,\boldsymbol{x}+\delta\boldsymbol{A}\,\delta\boldsymbol{x}=\delta\boldsymbol{b} $$

将含 $\delta\boldsymbol{x}$ 的项合并：

$$ (\boldsymbol{A}+\delta\boldsymbol{A})\delta\boldsymbol{x}=\delta\boldsymbol{b}-\delta\boldsymbol{A}\,\boldsymbol{x} $$

两边左乘 $\boldsymbol{A}^{-1}$，得：

$$ (\boldsymbol{I}+\boldsymbol{A}^{-1}\delta\boldsymbol{A})\delta\boldsymbol{x} =\boldsymbol{A}^{-1}\delta\boldsymbol{b}-\boldsymbol{A}^{-1}\delta\boldsymbol{A}\,\boldsymbol{x} $$

因此：

$$ \delta\boldsymbol{x}=(\boldsymbol{I}+\boldsymbol{A}^{-1}\delta\boldsymbol{A})^{-1} \big(\boldsymbol{A}^{-1}\delta\boldsymbol{b}-\boldsymbol{A}^{-1}\delta\boldsymbol{A}\,\boldsymbol{x}\big) $$

当 $\|\boldsymbol{A}^{-1}\|\,\|\delta\boldsymbol{A}\|<1$ 时，利用 Neumann 级数估计有：

$$ \|(\boldsymbol{I}+\boldsymbol{A}^{-1}\delta\boldsymbol{A})^{-1}\|\le\frac{1}{1-\|\boldsymbol{A}^{-1}\|\,\|\delta\boldsymbol{A}\|} $$

于是取范数并利用三角不等式得到：

$$ \begin{align*} \|\delta\boldsymbol{x}\| &\le \frac{1}{1-\|\boldsymbol{A}^{-1}\|\,\|\delta\boldsymbol{A}\|}\, \Big(\|\boldsymbol{A}^{-1}\|\,\|\delta\boldsymbol{b}\|+\|\boldsymbol{A}^{-1}\|\,\|\delta\boldsymbol{A}\|\,\|\boldsymbol{x}\|\Big)\\ &= \frac{\|\boldsymbol{A}^{-1}\|}{1-\|\boldsymbol{A}^{-1}\|\,\|\delta\boldsymbol{A}\|}\, \Big(\|\delta\boldsymbol{b}\|+\|\delta\boldsymbol{A}\|\,\|\boldsymbol{x}\|\Big). \end{align*} $$

两边除以 $\|\boldsymbol{x}\|$，并用不等式：

$$ \|\boldsymbol{x}\|\ge\frac{\|\boldsymbol{b}\|}{\|\boldsymbol{A}\|} $$

（由 $\boldsymbol{A}\boldsymbol{x}=\boldsymbol{b}$ 和算子范数性质）可得：

$$ \begin{align*} \frac{\|\delta\boldsymbol{x}\|}{\|\boldsymbol{x}\|} &\le \frac{\|\boldsymbol{A}^{-1}\|}{1-\|\boldsymbol{A}^{-1}\|\,\|\delta\boldsymbol{A}\|}\, \Big(\frac{\|\delta\boldsymbol{b}\|}{\|\boldsymbol{x}\|}+\|\delta\boldsymbol{A}\|\Big)\\ &\le \frac{\|\boldsymbol{A}^{-1}\|}{1-\|\boldsymbol{A}^{-1}\|\,\|\delta\boldsymbol{A}\|}\, \Big(\frac{\|\boldsymbol{A}\|\,\|\delta\boldsymbol{b}\|}{\|\boldsymbol{b}\|}+\|\delta\boldsymbol{A}\|\Big) \end{align*} $$

将 $\|\delta\boldsymbol{A}\|$ 表示为 $\|\boldsymbol{A}\|\dfrac{\|\delta\boldsymbol{A}\|}{\|\boldsymbol{A}\|}$ 并整理得到：

$$ \frac{\|\delta\boldsymbol{x}\|}{\|\boldsymbol{x}\|} \le \frac{\|\boldsymbol{A}\|\,\|\boldsymbol{A}^{-1}\|}{1-\|\boldsymbol{A}\|\,\|\boldsymbol{A}^{-1}\|\,\dfrac{\|\delta\boldsymbol{A}\|}{\|\boldsymbol{A}\|}} \left(\frac{\|\delta\boldsymbol{A}\|}{\|\boldsymbol{A}\|}+\frac{\|\delta\boldsymbol{b}\|}{\|\boldsymbol{b}\|}\right) $$

在仅有系数矩阵扰动且 $\|\boldsymbol{A}^{-1}\|\,\|\delta\boldsymbol{A}\|<1$ 的情形，令 $\delta\boldsymbol{b}=\boldsymbol{0}$，则上式化为：

$$ \frac{\|\delta\boldsymbol{x}\|}{\|\boldsymbol{x}\|} \le \frac{\|\boldsymbol{A}\|\,\|\boldsymbol{A}^{-1}\|\,\dfrac{\|\delta\boldsymbol{A}\|}{\|\boldsymbol{A}\|}} {1-\|\boldsymbol{A}\|\,\|\boldsymbol{A}^{-1}\|\,\dfrac{\|\delta\boldsymbol{A}\|}{\|\boldsymbol{A}\|}} =\frac{\|\boldsymbol{A}^{-1}\|\,\|\delta\boldsymbol{A}\|} {1-\|\boldsymbol{A}^{-1}\|\,\|\delta\boldsymbol{A}\|} $$

上式表明，当系数矩阵有扰动时，解的扰动仍与 $\|\boldsymbol{A}\|\,\|\boldsymbol{A}^{-1}\|$ 有关，$\|\boldsymbol{A}\|\,\|\boldsymbol{A}^{-1}\|$ 越大，解的扰动也越大。

总的来说，若系数矩阵 $\boldsymbol{A}$ 有扰动 $\delta\boldsymbol{A}$，右端项有扰动 $\delta\boldsymbol{b}$，相应的解 $\boldsymbol{x}$ 有扰动 $\delta\boldsymbol{x}$，且 $\|\boldsymbol{A}^{-1}\|\,\|\delta\boldsymbol{A}\|<1$，则有：

$$ \frac{\|\delta\boldsymbol{x}\|}{\|\boldsymbol{x}\|} \le \frac{\|\boldsymbol{A}\|\,\|\boldsymbol{A}^{-1}\|} {1-\|\boldsymbol{A}\|\,\|\boldsymbol{A}^{-1}\|\,\dfrac{\|\delta\boldsymbol{A}\|}{\|\boldsymbol{A}\|}} \left(\frac{\|\delta\boldsymbol{A}\|}{\|\boldsymbol{A}\|}+\frac{\|\delta\boldsymbol{b}\|}{\|\boldsymbol{b}\|}\right) $$

综合上述各式可以得出，当系数矩阵或右端项有扰动时，数 $\|\boldsymbol{A}\|\,\|\boldsymbol{A}^{-1}\|$ 控制了解的扰动程度。也就是说，数 $\|\boldsymbol{A}\|\,\|\boldsymbol{A}^{-1}\|$ 可以反映方程组的状态。

对于非奇异矩阵$\boldsymbol{A}$，称数 $\|\boldsymbol{A}\|\,\|\boldsymbol{A}^{-1}\|$ 为矩阵$\boldsymbol{A}$的条件数，记为：

$$ \mathrm{cond}(\boldsymbol{A})= \|\boldsymbol{A}\|\,\|\boldsymbol{A}^{-1}\| $$

[!NOTE]

在题目未明确说明的情况下，通常在计算时使用2-范数。