ガウスの消去法（MATLAB バージョン）

どうも，筆者です．部分ピボット選択付きのガウスの消去法を C 言語と Fortran で書いてきたが，MATLAB のスクリプトを探している人もいた．

目的

なので，今回は，MATLAB 用にガウスの消去法を書き直すことにした．これも久しぶりに使うので，もっと効率化できる部分があるかもしれない．まぁ，MATLAB 用にわざわざ 1 から書く人はいない（組み込みの関数を使うべきである）と思うので，勉強用として考えれば，速度と精度はそこまで必要にならないだろう．

C 言語と Fortran は，以下でコードを公開している．

workspacememory.hatenablog.com workspacememory.hatenablog.com

MATLAB 用のコード

MATLAB 用に書いていて気付いたが，Fortran も同様なコードに落とせるのではないか（確か，ベクトル演算できたよなぁ）．

%===================================%
%                                   %
% gaussSolver                       %
%                                   %
% A: coefficient matrix, A(row,col) %
% b: right side vector,  b(row)     %
%                                   %
% return: solution vector x         %
%===================================%
function [x] = gaussSolver(A,b)
    [n, ~] = size(A);

    for i = 1:(n - 1)
        % i 列目のうち，最大のものを探す
        [~, id] = max(abs(A(i:n, i)));

        if id ~= 1
            pivot = id + i - 1;
            % 行の入れ替え
            tmpVec = A(i, :);
            A(i, :) = A(pivot, :);
            A(pivot, :) = tmpVec;

            % 右辺ベクトルの入れ替え
            tmpScalar = b(i);
            b(i) = b(pivot);
            b(pivot) = tmpScalar;
        end
        diagVal = A(i, i);

        % 対角成分が eps よりも小さい場合は，エラー処理
        if (abs(diagVal) < eps)
            error('Error: diagonal component is less than eps=%.15e\n', eps);
        end

        % ========
        % 前進消去
        % ========
        for row = (i + 1):n
            idList = i:n;
            weight = A(row, i) / diagVal;
            A(row, idList) = A(row, idList) - weight * A(i, idList);
            b(row) = b(row) - weight * b(i);
        end
    end

    % ========
    % 後退代入
    % ========
    b(n) = b(n) / A(n, n);
    for row = (n - 1):(-1):1
        idList = (row + 1):n;
        b(row) = (b(row) - A(row, idList) * b(idList)) / A(row, row);
    end
    x = b;
end

動作確認

筆者は，MATLAB を持っていないので，Octave を利用して，検証を行った．その際，テストスクリプトとして，以下のようなものを作成した．

n = 10;
A = rand(n);
xexact = (1:n)';
b = A * xexact;

xhat = gaussSolver(A, b)

これを，Octave で実行したところ，真の解の近似解が得られたので，正しく動作していると考えられる．