読者です 読者をやめる 読者になる 読者になる

作業中のメモ

よく「計算機」を使って作業をする.知らなかったことを中心にまとめるつもり.

ガウスの消去法(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 で実行したところ,真の解の近似解が得られたので,正しく動作していると考えられる.