解线性方程组——直接解法：(Gauss)高斯消去法、列主元、全主元 | 北太天元

一、问题描述

对于线性方程组
$$Ax=b,A=\left(\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ ⋮& ⋮& & ⋮\\ a_{n1}& a_{n2}& \cdots & a_{nn}\end{array}\right),b=\left(\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_n\end{array}\right)$$

A为非奇异矩阵，求向量$x$

二、Gauss消去法

思路：

系数矩阵逐步消元得到上三角矩阵，然后使用回代算法求解。

1. (耿直版)消元法

$$\{\begin{array}{l}{a}_{11}^{(1)}{x}_1+a_{12}^{(1)}{x}_2+...+a_{1n}^{(1)}{x}_n=b_1^{(1)}\\ a_{21}^{(1)}{x}_1+a_{22}^{(1)}{x}_2+...+a_{2n}^{(1)}{x}_n=b_2^{(1)}\\ \cdots \\ a_{n1}^{(1)}{x}_1+a_{n2}^{(1)}{x}_2+...+a_{nn}^{(1)}{x}_n=b_n^{(1)}\end{array}$$

简记为$A^{(1)}x=b^{(1)}$，其中$A^{(1)}=A,b^{(1)}=b$

若$a_{11}^{(1)}\ne 0$，令$l_{i1}=\frac{a_{i1}^{(1)}}{a_{11}^{(1)}}$，

令$-l_{i1}$乘第$1$行加到第$i$行$(i=2,3,\cdots ,n)$并保留第$1$个方程。

此时 $A^{(1)}\to A^{(2)},b^{(1)}\to b^{(2)}$
$$\begin{gathered} a_{ij}^{(2)}=a_{ij}^{(1)}-l_{i1}{a}_{1j}^{(1)}(i,j=2,3,\cdots ,n) \\ {b}_i^{(2)}=b_i^{(1)}-l_{i1}{b}_1^{(1)}(i=2,3,\cdots ,n) \end{gathered}$$

若$a_{kk}^{(k)}\ne 0$，令$l_{ik}=\frac{a_{ik}^{(k)}}{a_{kk}^{(k)}}$，
令$-l_{ik}$乘第$k$行加到第$i$行$(i=k+1,k+2,\cdots ,n)$
$$a_{ij}^{(k+1)}=a_{ij}^{(k)}-l_{ik}{a}_{kj}^{(k)}(i,j=k+1,k+2,\cdots ,n)$$ $$b_i^{(k+1)}=b_i^{(k)}-l_{ik}{b}_k^{(k)}(i=k+1,k+2,\cdots ,n)$$

缺点

$a_{kk}^{(k)}\ne 0$ 若为0，不能使用

$\left|a_{kk}^{(k)}\right|\approx 0$ 或相对其他元素较小时，舍入误差增加，可能导致误差剧增

进行每一步消元前，先找到列中最大元素所在的那一行，并与原来的行交换位置

2. 列主元

进行每一步消元前，先找到列中最大元素所在的那一行，并与原来的行交换位置

3. 全主元

进行每一步消元前，先找到剩余矩阵中最大的元素，做对应的行交换与列交换，
由于进行了列交换，最后要把$X$的位置对应回去

三、北太天元源代码

耿直版

function [X,Ae,be] = gsem_base(A,b)
% Gauss消去法耿直版
% A : 系数矩阵
% b : 右端常数 列向量
% X : 求得的解向量
% Ae: 得到的上三角矩阵
% be: 对应的右端项
% 使用要求，A1(k,k)在计算过程中不出现0
% 缺点：当A1(k,k)较小或近似于0时，计算误差可能增大很多
%
%   Version:            1.0
%   last modified:      09/09/2023
A1 = [A,b];%  A和b的增广矩阵
n = length(A);
n1 = length(A1);
% 系数矩阵 化为上三角的过程如下
t = 0;
for i = 1:n-1
    if A1(i,i) != 0
    	  % 对每一行进行操作
	       for k = i+1:1:n
            lik = A1(k,i)/A1(i,i);    % 算子
            A1(k,i:n1) = A1(k,i:n1)-lik*A1(i,i:n1);
%            T=A1 % 显示运算步骤
	       end
    else
        t = 1;
        fprintf('出现零元素，不能使用该函数');
        break;
    end
end
% 得到上三角后 进行回代
if t == 0
Ae = A1(:,1:n);
be = A1(:,n1);
X  = reg_utm(Ae,be);
end
end

保存为 gsem_base.m文件

列主元

function [X,Ae,be] = gsem_column(A,b)
% 列主元消去法
% A : 表示系数矩阵
% b : 表示右端常数 [列向量]
% X : 求得的解向量
% Ae: 得到的上三角矩阵
% be: 对应的右端项
% 避免了Gauss消去法中出现0的情况
%
%   Version:            1.0
%   last modified:      09/09/2023
A1 = [A,b];%  A和b的增广矩阵
n = length(A);
n1 = length(A1);
% 系数矩阵 化为上三角的过程如下
for i = 1:n-1
	 %[a,s]= max(A1(i:n,i))
    [~,s]= max(A1(i:n,i)); % s表示对应的位置,由于a用不到，用~替代
    s = s+i-1; % 由于s是相A对位置，故有此步
    A1([i,s],:) = A1([s,i],:); % 对应行交换位置
    % 对每一行进行操作
     for k = i+1:1:n
        lik = A1(k,i)/A1(i,i);    % 算子
        A1(k,i:n1) = A1(k,i:n1)-lik*A1(i,i:n1);
%        T=A1 % 显示运算步骤
     end
end
% 得到上三角后 进行回代
Ae = A1(:,1:n);
be = A1(:,n1);
X = reg_utm(Ae,be);
end

保存为 gsem_column.m文件

全主元

写全主元时,
注意,进行列的交换后, 对应解的位置发生了互换,最后要换回来.

function [X,Ae,be] = gsem_complete(A,b)
% 全主元消去法
% A : 系数矩阵
% b : 右端常数 [列向量]
% X : 求得的解向量
% Ae: 得到的上三角矩阵
% be: 对应的右端项
% 比列主元更加稳定了
%
%   Version:            1.0
%   last modified:      09/09/2023
A1 = [A,b];    %A和b的增广矩阵
n = length(A);
n1 = length(A1);
X = zeros(n,1);
t = 1:n; % t 用于修正解的位置与原方程不匹配
% 系数矩阵化为上三角的过程如下
for i = 1:n-1
	   A = A1(i:n,i:n);
    [m,a,b] = max_loc(A); %这里不用A1，避免误判b中元素
    % a,b是相对的位置
    a = a+i-1;
    b = b+i-1;
    A1([i,a],:) = A1([a,i],:); % 对应行交换位置
    A1(:,[i,b]) = A1(:,[b,i]); % 对应列交换位置
    [t(i),t(b)]= deal(t(b),t(i));
    % 对每一行进行操作
	    for k = i+1:1:n
        lik = A1(k,i)/A1(i,i);    % 算子
        A1(k,i:n1) = A1(k,i:n1)-lik*A1(i,i:n1);
%        T=A1 % 显示运算步骤
     end
end
% 得到上三角后进行回代
Ae = A1(:,1:n);
be = A1(:,n1);
X([t])= reg_utm(Ae,be); %使解匹配原方程位置
end

保存为gsem_complete.m文件
其中用到了我写的 max_loc.m

function [m,a,b] =max_loc(A)
% 获取矩阵中最大元素的位置
% m：矩阵中最大元素的值
% [a,b] 最大值在矩阵中的位置
%
%   Version:            1.0
%   last modified:      05/16/2023
    [r,l] = size(A);
    if r == 1    % 行向量
        a = 1;
        [m,b] = max(A);
    elseif l == 1    % 列向量
        b = 1;
        [m,a] = max(A);
    else
        [M,I] = max(A);
        [m,b] = max(M); % b = 最大元素所在的列
        a = I(b); % a = 最大元素所在的行
    end
end

保存为max_loc.m文件

四、简单测试一下

% 线性方程组的例子
%   last modified:      09/09/2023
clc;clear all;
%%
A = [10 1 2 3 4; 1 9 -1 2 -3; 2 -1 7 3 -5; 3 2 3 12 -1; 4 2 3 4 15];
b = [12; -27; 14; -17; 10];

X2 = gsem_column(A,b)
X1 = gsem_base(A,b)
X3 = gsem_complete(A,b)

文中三个函数都用到了解上三角的回代算法，
参考之前写的：解上三角、回代算法|北太天元