MATLAB使用高斯消元法计算方程组的解

function X = uptrbk(A,B)
% A,B是系数矩阵和列向量
% 求方阵A 含多少行(列) X是N*1列向量解
[~, N] = size(A);
X = zeros(N,1);
% C一行,N+1列
C = zeros(1,N+1);
% 增广矩阵Aug
Aug = [A,B];
%循环从第一列到倒数第二列
for p=1:N-1
    [~,j] = max(abs(Aug(p:N,p)));%返回每一列中的绝对值最大的数的索引j,Y是最大值
    C = Aug(p,:);%把p行所有元素存储到C中,1,2,3,...,N-1.
    %进行交换
    Aug(p,:)=Aug(j+(p-1),:); %Aug第p行改为j+p-1行,(相对索引).
    Aug(j+(p-1),:)=C;
    if Aug(p,p) == 0   %如果交换后,Aii为0
        'A was singular. No unique solution';
        break
    end
    %交换后 进行计算
    for k = p+1:N
        m = Aug(k,p)/Aug(p,p);
        Aug(k,p:N+1) = Aug(k,p:N+1)-m*Aug(p,p:N+1);
    end
end
% 通过回代计算方程的解
    for i = N:-1:1  
        X(i) = (Aug(i, N+1) - Aug(i, i+1:N) * X(i+1:N)) / Aug(i, i);  
    end  
end  


% function X = uptrbk(A, B)  
%     % uptrbk Solves the linear equation system A*X = B using back substitution  
%     % A: Upper triangular matrix  
%     % B: Column vector  
%     
%     % 获取矩阵 A 的维度  
%     [N, ~] = size(A);  
%     % 初始化解向量 X  
%     X = zeros(N, 1);  
%     % 初始化增广矩阵 Aug  
%     Aug = [A, B];  
%     
%     % 高斯消元过程  
%     for p = 1:N-1  
%         % 寻找列 p 从当前行到最后一行的最大值  
%         [Y, j] = max(abs(Aug(p:N, p)));  
%         
%         % 交换当前行与最大值所在行  
%         if j+p-1 ~= p  
%             Aug([p, j+p-1], :) = Aug([j+p-1, p], :);  % 进行行交换  
%         end  
%         
%         % 检查非零主元  
%         if Aug(p, p) == 0  
%             error('Matrix A is singular. No unique solution exists.');  
%         end  
%         
%         % 消元过程  
%         for k = p+1:N  
%             m = Aug(k, p) / Aug(p, p);  
%             Aug(k, p:N+1) = Aug(k, p:N+1) - m * Aug(p, p:N+1);  
%         end  
%     end  
%     
%     % 回代过程，解出 X  
%     X(N) = Aug(N, N+1) / Aug(N, N);  
%     for i = N-1:-1:1  
%         X(i) = (Aug(i, N+1) - Aug(i, i+1:N) * X(i+1:N)) / Aug(i, i);  
%     end  
% end

其中,A是系数矩阵,B是右边列向量.