Gauss列主元素法求解方程组Ax=b

原理

在第k次消元前，首先选取此次消元的主元素. 若$ \left| {a_{lk}^{(k)}} \right| = \mathop {\max }\limits_{k \le i \le n} \left| {a_{ik}^{(k)}} \right|; $则选取$ {a_{lk}^{(k)}} $为列主元素。交换增广矩阵$ [{A^{(k)}},{b^{(k)}}] $的第l行和第k行（若l=k则不必交换），然后进行消元。

MATLAB实现

函数gausscol.m

function [A,b,x] = gausscol(A,b)
%% Gauss列主元法
cnt = length(A); % 矩阵A的行数
for i = 1:cnt - 1 % 从第1行到第cnt-1行检索
    temp = A(i,i); % 提取当前对角线元素
    row = i; % 记录当前行数
    for j = i + 1:cnt  % 从i的下一行开始，寻找最大列
        if abs(A(j,i)) > abs(temp)
            temp = A(j,i);
            row = j;
        end
    end
    if row ~= i
        for j = i:cnt % 交换A的第i行与第row行
            m = A(i,j);
            A(i,j) = A(row,j);
            A(row,j) = m;
        end
        m = b(i); % 交换b的第i行与第row行
        b(i) = b(row);
        b(row) = m;
    end
    
    for j = i + 1:cnt
        temp = A(j,i) / A(i,i); %l
        for k = i:cnt
            A(j,k) = A(j,k) - temp * A(i,k); % 消元
        end
        b(j) = b(j) - b(i) * temp;
    end
end

%% 回代
for i = cnt:-1:1 % 从倒数第二行开始回代
    temp = cnt; % 记录开始回代的行号
    while temp > i
        b(i) = b(i) - A(i,temp) * x(temp);
        temp = temp - 1;
    end
    x(i) = b(i)/A(i,i);
end

测试用例

% A = [0.729 0.81 0.9;1 1 1;1.331 1.21 1.1];
% b = [0.8338;0.8338;1];
% A = [3.4129 4.9317 8.7643 1.3142;0.59301 1.7875 2.5330 1.5435 ; 1.1348 3.8326 1.1651 3.4017; 1.2371 4.9998 10.6271 0.0147];
% b = [18.4231;6.3941;9.5342;16.9237];