求线性方程组解的直接方法 - 老官童鞋gogo的博客

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

其中，为常数，为未知量。上式也可写成矩阵形式：

如果线性方程组的系数矩阵的行列式不为，即，线性方程组存在唯一的解。当较小时，Cramer法则给出了线性方程组的公式解：

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

一、消去法

1、高斯（Gauss）消去法

(1) 高斯消去法原理

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

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

其中：

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

第一步：设，记，将方程组中第个方程减去第个方程乘以，完成第一次消元，得同解方程组

其中：

方程组简记为

第二步：设 , 记。将方程组(2-4)中第个方程减去第个方程乘以，完成第二次消元。

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

简记为：

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

(2) 高斯消去法代码

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

1
function x = gauss_solve(A, b)
2
    n = length(b);
3
    for k = 1:n-1
4
        for i = k+1:n
5
            if abs(A(i,k)) > abs(A(k,k))
6
                A([k i], :) = A([i k], :);
7
                b([k i]) = b([i k]);
8
            end
9
        end
10
        if abs(A(k,k)) < eps
11
            x = [];
12
            return;
13
        end
14
        m = A(k+1:n, k) / A(k,k);
15
        A(k+1:n, k) = m;
16
        b(k+1:n) = b(k+1:n) - m * b(k);
17
        A(k+1:n, k+1:n) = A(k+1:n, k+1:n) - m * A(k, k+1:n);
18
    end
19
    if abs(A(n,n)) < eps
20
        x = [];
21
        return;
22
    end
23
    x = zeros(n,1);
24
    x(n) = b(n) / A(n,n);
25
    for i = n-1:-1:1
26
        x(i) = (b(i) - A(i, i+1:n) * x(i+1:n)) / A(i,i);
27
    end
28
end

(3) 高斯消去法计算量

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

乘法次数：

除法次数：

在回代过程中，计算需要次乘除法，整个回代过程需要乘除运算的总量为：

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

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

2、主元素消去法

(1) 主元素消去法原理

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

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

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

第一步：首先在矩阵的第列中选取绝对值最大的元，比如，则。将增广矩阵中第行与第行互换。为方便起见，记行互换后的增广矩阵为，然后进行第一次消元，得矩阵：

第二步：在矩阵的第列中选主元，比如，使：

将矩阵的第行与第行互换，再进行第二次消元，得矩阵。

第三步：在矩阵的第列中选主元，如使：

将的第行与第行互换，进行第次消元。如此经过步，增广矩阵被化成上三角形矩阵，最后由回代过程求解。

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

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

1
function x = gauss_solve(A, b)
2
    n = length(b);
3
    for k = 1:n-1
4
        [~, p] = max(abs(A(k:n, k)));    % 寻找最大主元所在行数
5
        p = p + k - 1;
6
        if p ~= k              % 交换
7
            A([k p], :) = A([p k], :);
8
            b([k p]) = b([p k]);
9
        end
10
        tol = eps * norm(A, inf);      % 防止因为数值稳定性问题对是否有解的判定造成影响
11
        if abs(A(k,k)) < tol
12
            x = [];
13
            return;
14
        end
15
        m = A(k+1:n, k) / A(k,k);
16
        A(k+1:n, k) = m;
17
        b(k+1:n) = b(k+1:n) - m * b(k);
18
        A(k+1:n, k+1:n) = A(k+1:n, k+1:n) - m * A(k, k+1:n);
19
    end
20
    if abs(A(n,n)) < tol
21
        x = [];
22
        return;
23
    end
24
    x = zeros(n,1);
25
    x(n) = b(n) / A(n,n);
26
    for i = n-1:-1:1
27
        x(i) = (b(i) - A(i, i+1:n) * x(i+1:n)) / A(i,i);
28
    end
29
end

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

二、直接三角分解法

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

(1) LU分解原理

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

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

那么第次消元可以使用下面增广矩阵左乘得到：

经过次消元后，得到：

因为为非奇异矩阵，存在逆矩阵，可以得到：

于是有：

令：

可以得到：

则有：

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

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

高斯消去法要求主元元素，将这一条件与矩阵自身的性质联系起来，可以得到定理如下：

设为阶方阵，若的顺序主子式均不为零，则矩阵存在唯一的LU分解。

(2) LU分解代码

下面是计算得出矩阵的MATLAB代码：

1
function [L, U] = lu_decomp(A)
2
    n = size(A,2);
3
    U = A;
4
    L = eye(n);
5
    for k = 1:n-1
6
        L(k+1:n, k) = U(k+1:n, k) / U(k, k);
7
        U(k+1:n, k+1:n) = U(k+1:n, k+1:n) - L(k+1:n, k) * U(k, k+1:n);
8
        U(k+1:n, k) = 0;
9
    end
10
end

下面是通过矩阵求解方程组的MATLAB代码：

1
function x = LU_solve(L, U, b)
2
    n = size(b, 1);
3
    y = zeros(n, 1);
4
    x = zeros(n, 1);
5
    for i = 1:n
6
        y(i) = b(i) - L(i, 1:i-1) * y(1:i-1);
7
    end
8
    for i = n:-1:1
9
        x(i) = (y(i) - U(i, i+1:n) * x(i+1:n)) / U(i, i);
10
    end
11
end

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

从直接三角分解公式知，当时，计算将中断；或当很小时，按公式计算会使舍入误差扩散。但若矩阵非奇异，可增加列选主元过程，即先进行一系列的行交换，保证LU分解能够实现。从而有如下定理：

设矩阵非奇异, 则存在置换矩阵，使得有唯一的 Doolittle 分解，其中是单位下三角形矩阵，是上三角形矩阵，且。即存在置换矩阵、单位下三角矩阵和上三角矩阵，使得。

我们一般不使用这种方法，但是下面给出一份代码，在已知的情况下求解线性方程组：

1
function x = solveWithPLU(P, L, U, b)
2
    n = size(b, 1);
3
    y = zeros(n, 1);
4
    x = zeros(n, 1);
5
    Pb = P * b;
6
    for i = 1:n
7
        y(i) = Pb(i) - L(i, 1:i-1) * y(1:i-1);
8
    end
9
    for i = n:-1:1
10
        x(i) = (y(i) - U(i, i+1:n) * x(i+1:n)) / U(i, i);
11
    end
12
end

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

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

(1) 原理

方程组：

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

且

则可唯一分解为：

其中由原矩阵给出。

NOTE

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

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

当矩阵三角分解后, 解线性方程组化为求解方程组和。

解，得：

再解，得三对角方程组的解：

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

(2) 代码

由于矩阵的大部分位置都是，因此系统不会给你整个矩阵，而是只给您传入一个长为的向量，以及长为的向量。分别代表上文中的数组，即：

实现求解的MATLAB代码如下：

1
function x = solve_tridiag(a, p, q, b)
2
    n = length(a);
3
    l = zeros(n, 1);
4
    u = zeros(n, 1);
5
    u(1) = a(1);
6
    for i = 1:n
7
        if i > 1
8
            l(i) = q(i-1) / u(i-1);
9
            u(i) = a(i) - l(i) * p(i-1);
10
        end
11
    end
12
    y = zeros(n, 1);
13
    y(1) = b(1);
14
    for i = 2:n
15
        y(i) = b(i) - l(i) * y(i-1);
16
    end
17
    x = zeros(n, 1);
18
    x(n) = y(n) / u(n);
19
    for i = n-1:-1:1
20
        x(i) = (y(i) - p(i) * x(i+1)) / u(i);
21
    end
22
end

2、平方根法

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

NOTE

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

且的对角元素皆为正数。矩阵的这种分解叫做Cholesky分解。

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

比较与的相应元素，可得：

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

当矩阵完成 Cholesky 分解后，求解方程组就转化为依次求解方程组：

它们的解分别为：

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

(2) 改进平方根法（法）

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

其中为单位下三角形矩阵，为对角矩阵。记：

然后利用待定系数法，导出分解的计算公式，对：

计算顺序如下：

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

则上式可改写成：

按上式进行分解，乘除运算量与 Cholesky 分解相当，且避免了开方运算。

矩阵作分解后，解方程组可分两步进行：先解方程组，再由求。具体计算公式为

求解线性方程组的这一方法称为改进平方根法，也叫法。

TIP

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

(3) 代码

下面是对矩阵做分解的MATLAB代码：

1
function L = decompLL(A)
2
    L = eye(size(A));
3
    n = size(A, 1);
4
    d = zeros(n, 1);
5
    for k = 1:n
6
        L(k, k) = sqrt(A(k, k) - L(k, 1:k-1) * (L(k, 1:k-1))');
7
        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);
8
    end
9
end

下面是利用法求解的MATLAB代码：

1
function x = solveWithLL(L, b)
2
    n = length(b);
3
    y = zeros(1, n);
4
    x = zeros(1, n);
5
    for k = 1:n
6
        y(k) = b(k) - L(k, 1:k-1) * (y(1:k-1))';
7
    end
8
    for k = n:-1:1
9
        x(k) = (y(k) - x(k+1:n) * L(k+1:n, k)) / L(k, k);
10
    end
11
end

下面是对矩阵做分解的MATLAB代码：

1
function [L, D] = decompLDL(A)
2
    D = zeros(size(A));
3
    L = eye(size(A));
4
    U = zeros(size(A));
5
    n = size(A, 1);
6
    d = zeros(n, 1);
7
    for k = 1:n
8
        d(k) = A(k, k) - U(k, 1:k-1) * (L(k, 1:k-1))';
9
        D(k, k) = d(k);
10
        U(k+1:n, k) = A(k+1:n, k) - U(k+1:n, 1:k-1) * (L(k, 1:k-1))';
11
        L(k+1:n, k) = U(k+1:n, k) / d(k);
12
    end
13
end

下面是利用法求解的MATLAB代码：

1
function x = solveWithLDL(L, D, b)
2
    n = length(b);
3
    y = zeros(1, n);
4
    x = zeros(1, n);
5
    y(1) = b(1);
6
    for k = 2:n
7
        y(k) = b(k) - L(k, 1:k-1) * (y(1:k-1))';
8
    end
9
    x(n) = y(n) / D(n, n);
10
    for k = n-1:-1:1
11
        x(k) = y(k) / D(k, k) - x(k+1:n) * L(k+1:n, k);
12
    end
13
end

四、误差分析

1、向量范数

(1) 定义

设为中的向量。函数称为一个范数，当且仅当对任意和任意标量，满足：

正定性（非负性且仅在零向量时为零）：
齐次性（绝对齐次）：
三角不等式：

这三条构成范数的公理。

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

对，定义向量的 -范数为

三个常用的情形：

1-范数（曼哈顿范数）：
2-范数（欧几里得范数）：
-范数（切比雪夫范数）：

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

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

在有限维向量空间中，存在常数使得对所有都有：

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

2、矩阵范数

(1) 定义

设（此处只研究方阵）。矩阵范数是定义在矩阵集合上的函数，若对任意方阵和标量满足：

, 且
（相容性条件）

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

给定向量范数（域为），对矩阵定义诱导范数（也称算子范数）：

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

(2) 常见矩阵向量

1-范数（矩阵）：

即各列绝对值和的最大值（最大列和）。
-范数（矩阵）：

即各行绝对值和的最大值（最大行和）。
2-范数（谱范数）：

这里表示对称正定矩阵的最大特征值。
Frobenius 范数（F-范数，希尔伯特–施密特范数）：

矩阵的F-范数可以认为是向量2-范数的直接推广，将阶方阵看做维向量所得到。可以证明，矩阵的F-范数与向量的2-范数相容，但F-范数不是由向量范数诱导出的矩阵范数，因为，而诱导出的矩阵范数上述值都为。

(3) 等价性与误差

如果将矩阵范数看作空间上的向量范数，则由向量范数的等价性可得矩阵范数的等价性。

矩阵的误差可用矩阵范数表示。设是的近似矩阵，、分别称为的关于范数的绝对误差与相对误差。

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

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

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

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

由，得，从而有：

两边取范数，得：

又因为对于任意向量有算子范数的性质：

从而：

这里用到。所以将上面两式相除，得到：

上式表明，当右端项有扰动时，解的相对误差不超过右端项的相对误差的倍。

如果充分小，使得，则由上式得：

首先由扰动方程出发：

展开得：

将含的项合并：

两边左乘，得：

因此：

当时，利用 Neumann 级数估计有：

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

两边除以，并用不等式：

（由和算子范数性质）可得：

将表示为并整理得到：

在仅有系数矩阵扰动且的情形，令，则上式化为：

上式表明，当系数矩阵有扰动时，解的扰动仍与有关，越大，解的扰动也越大。

总的来说，若系数矩阵有扰动，右端项有扰动，相应的解有扰动，且，则有：

综合上述各式可以得出，当系数矩阵或右端项有扰动时，数控制了解的扰动程度。也就是说，数可以反映方程组的状态。

对于非奇异矩阵，称数为矩阵的条件数，记为：

NOTE

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

老官童鞋gogo

一、消去法

1、高斯（Gauss）消去法

(1) 高斯消去法原理

(2) 高斯消去法代码

(3) 高斯消去法计算量

2、主元素消去法

(1) 主元素消去法原理

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

二、直接三角分解法

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

(1) LU分解原理

(2) LU分解代码

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

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

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

(1) 原理

(2) 代码

2、平方根法

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

(2) 改进平方根法（法）

(3) 代码

四、误差分析

1、向量范数

(1) 定义

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

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

2、矩阵范数

(1) 定义

(2) 常见矩阵向量

(3) 等价性与误差

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

目录