Matlab 实现解线性方程组的几类直接法

问题

求解线性方程组

$$ \left{ \begin{array}{c} a_{11}x_1+a_{12}x_2+…+a_{1n}x_n=b_1 \ a_{21}x_1+a_{22}x_2+…+a_{2n}x_n=b_2 \
… \ a_{n1}x_1+a_{n2}x_2+…+a_{nn}x_n=b_n \end{array} \right. $$

其中方程组系数 $a_{ij}(i,j=1,2,3,…,n)$ 与右端项 $b_i(i=1,2,3,…,n)$ 均为实数,且 $b_1,b_2,…,b_n$ 不全为零,上述方程组简记为 $$ Ax=b $$ 其中 $$ A=\begin{bmatrix} a_{11} & a_{12} & … & a_{1n} \ a_{21} & a_{22} & … & a_{2n} \ …&…&…&… \ a_{n1} & a_{n2} & … & a_{nn} \ \end{bmatrix}, x=\begin{bmatrix} x_{1}\ x_{2}\ … \ x_{n}\ \end{bmatrix}, b=\begin{bmatrix} b_{1}\ b_{2}\ … \ b_{n}\ \end{bmatrix}. $$ 解线性方程组直接法(也称精确解法)的定义为:

若计算过程中没有舍入误差,则经过有限次的算术运算可以求得方程组的准确解的方法。

其主要有两类方法:

  • Gauss 消去法
    • Gauss 消去法
    • Gauss 列主元消去法
    • Gauss 按比例列主元消去法
    • Gauss-Jordan 消去法
  • 直接三角分解法
    • Corut 方法
    • $L^TDL$ 分解
    • 追赶法

本篇只用 Matlab 实现上述算法,不涉及算法的介绍与证明,如果对证明与原理比较感兴趣,详情参考:

并且事实上这些算法在 Matlab 中都有既成的函数或者工具包,虽然 CS 里面有一句话:“如果有现成的轮子那么绝对不自己造。”很有道理,但是对于算法的理解还是需要亲自尝试,这样才能比较好的理解。

在最初编写这些算法时,我的 Matlab 代码风格还没有确定,以至相比于后续算法的源码,这些源码会显得比较的不规范,我会进行一部分的修正,但是毕竟时间比较久,修正程度可能比较有限,如果有一些问题,欢迎联系。

正文

Gauss 消去法

Gauss 消去法

function x = Gaussian_elimination( A, b )
% 高斯消去法解线性方程组
% Input: A: 系数矩阵; b: 方程组右端项;
% Output: x: 线性方程组解;

n = size(A, 1); % 矩阵 A 阶数

if size(A,1)~=size(A,2) || size(A,1)~=n
    error('Matrix does not match') 
end

for k = 1 : n - 1
    if A(k, k) == 0 % 检测 A 是否为正定矩阵
        for i = k + 1 : n
            if A(i, k) ~= 0
                for j = k : n
                    t = A(k, j); A(k, j) = A(i, j); A(i, j) = t;
                end % for
            else
                error('A is singular.');
            end % if
        end % for
    end % if
    for i = k + 1 : n
        l(i, k) = A(i, k) / A(k, k);
        for j = k : n
            A(i, j) = A(i, j) - l(i, k) * A(k, j);
        end % for
        b(i) = b(i) - A(i, k) * b(k);
    end % for
end % for

x(n) = b(n) / A(n, n);
for k = n - 1 : -1 : 1
    sum = 0;
    for i = k + 1 : n
        sum = sum + A(k, j) * x(j);
    end % for
    x(k) = (b(k) - sum) / A(k, k);
end % for

end % function

Gauss 消去法就非常的基础了,不过是平时我们解方程的固有过程,用来练习一下 Matlab 基本语言的使用。

当初在做这些课题的时候对代码的健壮性考虑的较少,主要是来理解算法的。如果有相应的需求,需要进一步的修改。

Gauss 列主元消去法

function x = Gaussian_principal_element_elimination( A, b )
% 高斯列主元消去法解线性方程组
% Input: A: 系数矩阵; b: 方程组右端项;
% Output: x: 线性方程组解;

n = size(A, 1); % 矩阵 A 阶数

if size(A, 1) ~= size(A, 2) || size(A, 1) ~= n
    error('Matrix does not match') 
end % if

for i = 1 : n
    [ap, p] = max(abs(A(i : n, i)));
    p = p + i - 1;
    
    if ap == 0 
        error('divide by zero!')
    end % if
    if p > i
        t = A(i, :);
        A(i, :) = A(p, :);
        A(p, :) = t;
        
        t = b(i);
        b(i) = b(p);
        b(p) = t; 
    end % if
    %消元计算
    m = A(i+1:n, i) ./ A(i, i);
    A(i+1:n, i+1:n) = A(i+1:n, i+1:n) - m * A(i, i+1:n);
    b(i+1:n) = b(i+1:n) - m * b(i);
    A(i+1:n, i) = zeros(n-i, 1);
end % for

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 % for
end % function

Gauss 按比例列主元消去法

function x = Gauss_proportional_column_principal_elimination( A, b )
% 高斯消去法解线性方程组
% Input: A: 系数矩阵; b: 方程组右端项;
% Output: x: 线性方程组解;

n = size(A, 1); % 矩阵 A 阶数

if size(A, 1) ~= size(A, 2) || size(A, 1) ~= n
    error('Matrix does not match') 
end

A_max = max(A, [], 2); % A 行最大值
for i = 1 : n
    s(i) = A_max(i);
    if s(i) == 0
        error('A is singular.')
    end % if
    p(i) = i;
end % for

for k = 1 : n - 1
    %%%%%%%%%% 待完成
    r = 1; % 解 abs(a(p(r), k)) / s(p(r)) = max_k<=1<=n (abs(a(p(i), k)) / s(p(i));
    %%%%%%%%%%
    if A(p(r), k) == 0
        error('A is singular.')
    end % if
    if k ~= r
        temp = p(r);
        p(k) = p(r);
        p(r) = temp;
    end % if
    for i = k + 1 : n
        A(p(i), k) = A(p(i), k) / A(p(k), k);
        for j = k + 1 : n
            A(p(i), j) = A(p(i), j) - A(p(i), k) * A(p(i), j);
        end % for
    end % for
end % for

b2(p(2)) = b(P(1));
for i = 2 : n
    sum = 0;
    for j = 1 : i - 1
        sum = sum + A(p(i), j) * b2(p(j));
    end % for
    b2(p(i)) = b(p(i)) - sum;
end % for

if A(p(n), n) == 0
    error('A is singular.')
end % if

x(n) = b2(p(n)) / A(p(n), n);
for i = n - 1 : -1 : 1
    sum = 0;
    for j = i + 1 : n
        sum = sum + a(p(i), j) * x(j);
    end % for 
    x(i) = (b2(p(i)) - sum) / a(p(i), i);
end % for 

end % function

Gauss-Jordan 消去法

Gauss-Jordan 消去法就是在上述任一方法之中,消去步骤中对主元行外的所有行进行操作,其他没有任何变化,代码依方法选择而变,在此不进一步考虑。

直接三角形分解法

Crout 方法:

function x = Court( A, b )
% Court  方法解线性方程组
% Input: A: 系数矩阵; b: 方程组右端项;
% Output: x: 线性方程组解;

n = size(A);
L = zeros(n, n);
U = zeros(n, n);
for i = 1 : n
    L(i, i) = 1;
end % for
for k = 1 : n
    for j = k : n
        U(k, j) = A(k, j) - sum(L(k, 1 : k - 1) .* U(1 : k - 1, j)');
    end % for
    for i = k + 1 : n
        L(i, k) = (A(i, k) - sum(L(i, 1 : k - 1) .* U(1 : k - 1 , k)')) / U(k, k);
    end % for
end % for

[L,U] = myLU(A);
[n,m] = size(A);
y(1) = b(1);
for i = 2 : n    
    for j = 1 : i - 1
        b(i) = b(i) - L(i, j) * y(j);
    end
    y(i) = b(i);
end % for
 
x(n) = y(n) / U(n, n);
for i = (n - 1) : -1 : 1
    for j = n : -1 : i + 1
        y(i) = y(i) - U(i, j) * x(j);
    end
    x(i) = y(i) / U(i, i);
end % for

end % function

Court 方法是三角分解法中最简单的一种,计算量会大一些,但该方法适用于任意情况的方程组。

$L^TDL$ 分解法

$L^TDL$ 分解法适用面比较窄:必须是实对称矩阵。此时使用这种方法可以比较好的减少计算量,增加计算速度,但实际原理与 LU 分解法类似。

若是追求比较通用的方法,LDU 分解法是 LU 分解法的优化,但整体上大同小异,在此不占用篇幅专门书写它们的算法了。

追赶法

追赶法的适用范围也相当的窄:三对角实线性方程组。其本质就是 $LU$ 分解法在这里的运用,解这样的方程非常的迅速,需求的计算量也非常的小。

对于追赶法有一篇比较好的文章可以参考:Thomas, L.H., Elliptic Problems in Linear Differential Equations over a Network. Watson Science Computer Laboratory Report, 1949.

其中追赶法在国外常称为“Thomas 方法”,在热传导方程等物理方面有着非常广泛的用处。

结尾

对于实线性方程组,通过恰当的变换,都可以得到任意精确度的数值解或者直接获得解析解,在数值处理上属于最简单的一个类型。

在数值分析中,这些过程需要交给计算机计算,所以就需要着重考虑其中会由计算机计算产生的误差。因此诞生出了诸多的方法并不断优化,避免了 7种情况 带来的计算误差,使得数据处理更加精确可信。

数值分析也研究了常微分方程组、函数逼近、超越方程寻根等问题。在这些问题的算法中更加体现出了数值分析的精妙逻辑与魅力,我们使用精妙的方法通过计算机去一步步逼近真实结果,以此用于实际问题中的复杂方程求解。

这是非常有现实价值与研究价值的,之后我也会分别总结这些方面涉及到的一些算法,本次的线性方程组就是一个简单的入门。