基本方法

设有线性方程组
$$Ax=b$$
当 $detA\ne 0$ 时，方程组解存在且唯一，对增广矩阵 $(A,b)=(A^{(1)},b^{(1)})$ 作行初等变换，化 $A^{(1)}$ 为上三角矩阵 $A^{(n)}$ ，同时 $b^{(1)}化为b^{(n)}$ ，这时与增广矩阵 $(A^{(n)},b^{(n)})$ 相应的线性方程组为上三角形方程组
$$A^{(x)}x=b^{(n)}$$
记上标为矩阵的行数，下标为(行,列)坐标，设 $a_{ii}^{(i)}\ne 0(i=1,2,\cdots ,n)$ ，则上式的解为
$$\{\begin{array}{l}{x}_n=\frac{b_n^{(n)}}{a_{nn}^{(n)}},\\ \\ x_i=\frac{-{\sum }_{j=i+1}^n{a}_{ij}^{(i)}{x}_j+b_i^{(i)}}{a_{ii}^{(i)}}(i=n-1,n-2,\cdots ,2,1)\end{array}$$
它便是原方程组的解，实现上述求解过程的方法称为Gauss消元法

Gauss列主元素消去法

为避免小主元做除数，在第k步 $k=(1,\cdots ,n-1)$ 消元时，在 $A^{(k)}$ 的第k列主对角元以下(含主对角元)的元素中挑选绝对值最大者，记对应行的行数为idx，则先交换第k行和第idx行后再进行消元计算，每一步都按列选主元的Gauss消元法称为Gauss列主元素消去法。

Matlab实现

求解系数矩阵为上三角或下三角的函数

function X = linear_solve(A,b)
%LINEAR_SOLVE 求解线性方程组 AX=b
%   A只能是上三角或者下三角矩阵
if istriu(A)
    %如果A是上三角
    s = size(A);
    r = s(1, 1); % 获取行数
    c = s(1, 2); % 获取列数
    X=zeros(1, c);% 初始化X
    for rn=r:-1:1
        a = A(rn, :);
        mask_a = ones(1, c);
        mask_a(rn) = 0;
        x = (b(rn) - sum(mask_a .* a .* X)) / a(rn);
        X(rn) = x;
    end
elseif istril(A)
    %如果A是下三角
    s = size(A);
    r = s(1, 1); % 获取行数
    c = s(1, 2); % 获取列数
    X=zeros(1, c);% 初始化X
    for rn=1:r
        a = A(rn, :);
        mask_a = ones(1, c);
        mask_a(rn) = 0;
        x = (b(rn) - sum(mask_a .* a .* X)) / a(rn);
        X(rn) = x;
    end
else
    error('A需要为上三角或者下三角');
end

end

Gauss消去法解线性方程组

function X = Gauss_solve(A, b, pri)
%GAUSS_SOLVE 利用高斯消去法解线性方程组
%   inputs:   
%   A, b为系数
%   pri:是否选主元，1为选主元，0为不选住院
%   outputs:
%   X:解的矩阵
%
b = reshape(b, [length(b), 1]); % 确保b是列矩阵
Ab = [A b]; %将A和b列合并
s = size(Ab);

if pri
   % 选主元
   for r=2:s(1)
       [m, idx_m] = max(abs(Ab(r-1:end, r-1))); % 获取要消去的当前列的主元索引
       idx_m = idx_m + r - 2;
       % 将主元行与当前行做交换
       temp = Ab(r-1, 1:end);
       Ab(r-1, 1:end) = Ab(idx_m, 1:end);
       Ab(idx_m, 1:end) = temp;
       % 进行消去
       ratio = Ab(r:end, r-1) / Ab(r-1, r-1);
       Ab(r:end, 1:end) = Ab(r:end, 1:end) - ratio * Ab(r-1, 1:end);
   end
else
   % 不选主元 
   for r=2:s(1)
       ratio = Ab(r:end, r-1) / Ab(r-1, r-1);
       Ab(r:end, 1:end) = Ab(r:end, 1:end) - ratio * Ab(r-1, 1:end);
   end
end

A = Ab(1:end, 1:end-1);
b = Ab(1:end, s(2));

X = linear_solve(A, b);

end