数值分析之直接法解线性方程组：高斯消去法、高斯乔丹消去法

高斯消去法（Gauss消去法）和高斯-乔丹消去法（Gauss-Jordan消去法）都很简单。代码如下：

高斯消去法：

function [x] = Gaussian_elim(a, b)
% a is a coefficient matrix, b is a column vector storing values of a*x,
% and here the return value x should be a column vector storing values of x.
% We should suppose that matrix a is not a singular matrix.

X = [a b];
n = size(X,1);     % the number of rows
x = zeros(n,1);    % x stores the values of x

for i = 1:n
    [c,p] = max(abs(X(i:end,i)), [], 1);    % find the index p of the max element in column i
                                            % notice that (abs(X(i:end,i)) is a new matrix)
    X([i,p+i-1],:) = X([p+i-1,i],:);        % swap column i and column p+i
    for j = i+1:n                           % to eliminate the ith element in rows i to n
        m = X(j,i) / X(i,i);                % get the multiplier
        X(j,:) = X(j,:) - m*X(i,:);         % eliminate the first non-zero element in row j
    end
end

% 将a,b分开
a = X(:,1:size(a,2));
b = X(:,size(X,2));

%回代求解x1,x2,...xn
x(n) = b(n,1)/a(n,n);
for i = n-1:-1:1
    x(i) = (b(i) - a(i,i+1:size(a,2)) * x(i+1:end,:)) / a(i,i);
end

end

高斯-乔丹消去法：

function [x] = GaussianJordan_elim(a, b)
% a is a coefficient matrix, b is a column vector storing values of a*x,
% and here the return value x should be a column vector storing values of x.
% We should suppose that matrix a is not a singular matrix.

X = [a b];
n = size(X,1);     % the number of rows
x = zeros(n,1);    % x stores the values of x

for i = 1:n
    [c,p] = max(abs(X(i:end,i)), [], 1);    % find the index p of the max element in column i
                                            % notice that (abs(X(i:end,i)) is a new matrix)
    X([i,p+i-1],:) = X([p+i-1,i],:);        % swap column i and column p+i
    for j = 1:n                             % to eliminate the ith element in rows i to n
        if j ~= i
            m = X(j,i) / X(i,i);            % get the multiplier
            X(j,:) = X(j,:) - m*X(i,:);     % eliminate the first non-zero element in row j
        else
            X(j,:) = X(j,:) / X(j,i);       % set the first non-zero element to 1
        end
    end
end

%得到 x1,x2,...xn
x = X(:,size(X,2));

end