2023-06-27《计算方法》- 陈丽娟 - 线性方程组的直接解法

Matlab计算方法高斯消元法矩阵分解

线性方程组的解法这一课题我们在高等代数中已经了解过，对于一个非奇异方阵，通过求解或者克莱姆法则均可以直接得到方程的精确解，但是上述方法计算量很大，难以在实际中应用，因此引出了本章的内容。

首先，我们给出问题的形式以方便叙述：

其中

问题中

已知，待求变量为

1. 高斯消元法

高斯消元法是一个经典的算法，其基本思想是先把方程组转化为一个上三角方程组，然后逆次序逐一求出未知量。

针对, 首先记增广矩阵。

Step 1. 从第2行到第n行中消去的系数. 为了做到这一点，我们从高代中知道，一行加上另一行的倍不影响最终的解，因此可以将第一行的加到第行，, 即是
.

Step 2. 依次重复Step 1中的计算，可以得到最终的三角矩阵.

Step 3. 回代得到解.

高斯消元法的时间复杂度为:

在消元的过程中，第步需要计算的除法显然有次，因有行和列，则有乘法和减法各次，因此乘法和除法的时间复杂度为

减法和回代过程的时间复杂度显然比消元的过程小，因此高斯消元法的时间复杂度为, 另外由于没有需要额外存储的数据，其空间复杂度即矩阵的大小.

=我们介绍了高斯消元法的基本步骤，那么当满足什么条件的时候高斯消元法起作用呢？这里给出一个定理来说明这个问题。

定理1:
高斯消元法求解的充分必要条件是的顺序主子式全不为0.

这是由于在高斯消元的过程中总是要求.

虽然理论上上述高斯消元法可以得到问题的解，但是当使用计算机编程时，由于不可避免的误差和有效数字，在某些情况下上述方法不能得到目标解，书中给出了如下例子：

例子1:
在三位有效数字下，给定方程组

精确解为

.
用上述消元法消去第二个式子的

得到

.
而若使用第二个式子消去第一个式子中的

, 则有

因此为了减小误差，给出
列主元高斯消元法:
Step 1. 找到第一列中最大元素所在的行, 交换第1行和第行。然后将第一行的加到第行，.

Step 2. 依次重复Step 1中的计算, 找第列中的最大元素, 可以得到最终的三角矩阵.

Step 3. 回代得到解.

附上代码：

%% 顺序、列主元高斯消元法 
% 输入系数矩阵A, 向量b, 是否使用列主元方法iscol 
% 返回求解向量x, 最终的上三角矩阵A1, 以及对应的向量b1 
function [x, A1, b1] = GaussLE(A, b, iscol) 
    x = b; 
    if iscol == 1 
        % 主元高斯消元法 
        AA = [A b]; 
        n = length(b); 
        for i = 1:(n-1) 
            if AA(i,i) == 0 
                % 主元素为0，不能进行高斯消元返回并报错 
                x = 0; 
                A1 = 0; 
                b1 = 0; 
                print("The element is 0."); 
                break 
            else 
                % 找到最大的主元 
                [M I] = max(AA(i:n,i)); 
                % 交换位置 
                I = I + i - 1; 
                temp = AA(I,:); 
                AA(I,:) = AA(i,:); 
                AA(i,:) = temp; 
                for j = (i+1):n 
                    mij = - AA(j,i) / AA(i,i); 
                    AA(j,:) = AA(j,:) + mij .* AA(i,:); 
                end 
            end 
        end 
        % 回代过程 
        x(n) = AA(n, n+1) / AA(n, n); 
        for i = flip(1:(n-1)) 
            x(i) = (AA(i, n+1) - sum(AA(i, (i+1):n) .* x((i+1):n)')) / AA(i, i); 
        end 
        A1 = AA(:,1:n); 
        b1 = AA(:,n+1); 
    else 
        % 顺序高斯消元法 
        AA = [A b]; 
        n = length(b); 
        for i = 1:(n-1) 
            if AA(i,i) == 0 
                % 主元素为0，不能进行高斯消元返回并报错 
                x = 0; 
                A1 = 0; 
                b1 = 0; 
                print("The element is 0."); 
                break 
            else 
                for j = (i+1):n 
                    mij = - AA(j,i) / AA(i,i); 
                    AA(j,:) = AA(j,:) + mij .* AA(i,:); 
                end 
            end 
        end 
        % 回代过程 
        x(n) = AA(n, n+1) / AA(n, n); 
        for i = flip(1:(n-1)) 
            x(i) = (AA(i, n+1) - sum(AA(i, (i+1):n) .* x((i+1):n)')) / AA(i, i); 
        end 
        A1 = AA(:,1:n); 
        b1 = AA(:,n+1); 
    end 
end 

虽然可以将判断是否列主元消元置于循环当中以减少代码的重复度，但是这样需要增加n次判断，因此还是分开写着两类消元法。
然后我们给出一个例子：

format long 
 
A = rand(100,100)*1000; 
b = rand(100,1)*1000; 
% 顺序消元 
[x1, A1, b1] = GaussLE(A, b, 0); 
% 列主元消元 
[x2, A2, b2] = GaussLE(A, b, 1); 
% 内置求解器 
x3 = linsolve(A, b); 
plot(abs(x1-x3)) 
hold on 
plot(abs(x2-x3)) 
legend(["顺序消元" "列主元消元"]) 

和matlab内置求解器的误差为
enter description here
这个不严谨的例子可以部分说明列主元是要比顺序消元更好。

2. 高斯消元法的矩阵形式和LU分解

由于在矩阵的某一行加上另一行的某倍实际上等价于矩阵左乘, 为基坐标, 为倍数, 为单位阵。因此前述的一系列变换相当于矩阵左乘了一系列矩阵：, 同时由于均为下三角矩阵，因此也是下三角矩阵，则原问题可表述为

其中

是消元后的系数（显然是上三角矩阵）。
利用三角矩阵的优良形式可以快速求得原问题的解：令

, 先求解

, 解得

, 再求解

那么如何通过分解得到呢？首先假设的对角线元素均为1, 则立即可得.

我们有

依次看，即对

的第2行有

对

的第3行有

对

的第r行有

对

的第r列有

上面我们说明了对于顺序主子式都大于0的矩阵，存在分解，使得是单位下三角矩阵，是上三角矩阵，下面我们说明这个分解是唯一的。设存在两个分解, 可得, 由于是上三角矩阵，是单位下三角矩阵，因此, 即, .

下面给出三角分解算法求解线性方程组的程序：

%% 求解线性方程组的LU分解 
% 输入系数矩阵A，向量b， 如果b为空，则只返回分解 
function [x, L, U] = LULE(A, b) 
    %% 先对A进行LU分解 
    x = 0; 
    n = size(A,1); 
    L = eye(n); 
    U = zeros(n, n); 
    U(1,:) = A(1,:); 
    L(:,1) = A(:,1) / U(1,1); 
    for i = 2:n 
        for j = i:n 
            U(i,j) = A(i,j) - L(i,1:(i-1)) * U(1:(i-1),j); 
            L(j,i) = (A(j,i) - L(j,1:(i-1)) * U(1:(i-1),i)) / U(i,i); 
        end 
    end 
    if length(b) == n 
        %% 然后计算 y 和 x 
        x = zeros(n,1); 
        y = x; 
        y(1) = b(1); 
        for i = 2:n 
            y(i) = b(i) - L(i,1:(i-1)) * y(1:(i-1)); 
        end 
        x(n) = y(n) / U(n,n); 
        for i = flip(1:(n-1)) 
            x(i) = (y(i) - U(i,(i+1):n) * x((i+1):n))/U(i,i); 
        end 
    end 
end 

测试程序

format long 
 
A = rand(100,100)*1000; 
b = rand(100,1)*1000; 
% 顺序消元 
[x1, A1, b1] = GaussLE(A, b, 0); 
% 列主元消元 
[x2, A2, b2] = GaussLE(A, b, 1); 
% 内置求解器 
x3 = linsolve(A, b); 
 
[x4, L, U] = LULE(A, b); 
plot(abs(x1-x3)) 
hold on 
plot(abs(x2-x3)) 
hold on 
plot(abs(x4-x3)) 
legend(["顺序消元" "列主元消元" "LU"]) 