[Matlab]Upper Triangularization & Back Substitution代码

用来做数值计算作业的代码，来自Numerical Methods Using Matlab (4th Edition) [John H. Mathews, Kurtis K. Fink]，改了一下注释，并添加输出消元过程的代码

uptrbk.m

function X = uptrbk(A, B)
    % initialization
    [N N] = size(A);
    X = zeros(N, 1);
    C = zeros(1, N+1);
    
    % construct augmented matrix
    Aug = [A B];
    
    disp('Aug before elimination');
    disp(Aug);
    
    for p=1:N-1
        % partial pivoting
        [Y,j] = max(abs(Aug(p:N, p)));
        C = Aug(p,:);
        Aug(p,:) = Aug(j+p-1,:);
        Aug(j+p-1,:) = C;
        
        disp('Aug after pivoting');
        disp(Aug);
        
        if Aug(p, p) == 0
            'A was singular';
            break
        end
        
        % Elimination
        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
        
        disp(['Aug after ', num2str(p), ' elimination']);
        disp(Aug);
    end
    
    X = backsub(Aug(1:N, 1:N), Aug(1:N, N+1));
end

backsub.m

function X = backsub(A, B)
    n = length(B);
    X = zeros(n, 1);
    X(n) = B(n)/A(n, n);
    
    for k = n-1:-1:1
        X(k) = (B(k) - A(k, k+1:n) * X(k+1:n)) / A(k,k);
    end
end

Note: 如果要使用分数输出结果，需要用format rat设置输出格式，默认是short（上次作业为了这个精度问题还换了C++……原来完全没必要，设format long就行了- -b）

 参考连接： http://www.mathworks.cn/cn/help/matlab/ref/format.html