連立方程式を計算機に解かせる上で情報系や物理系の学生は必ず目にするのが、ヤコビ法[Jacobi]とガウス・ザイデル法[Gauss-Seidel]かなと思います。今回は大学院の課題でmatlabでやったのを備忘録的に残しておきます。
python版は下記にあります。
忙しい人はgitへどうぞ。matlabファイルがあります。
ヤコビ法[Jacobi]とガウス・ザイデル法[Gauss-Seidel]とは
ヤコビ法[Jacobi]とガウス・ザイデル法[Gauss-Seidel]両方ともに連立方程式を計算機に解かせるためのアルゴリズムです。ガウス・ザイデル法は計算結果をすぐ使うような工夫がされているため、ヤコビ方に比べて計算がとんでもなく早いです。
直説法と反復法の違い
直説法は普通に行列で連立方程式を解く際に用いられる方法で、下記の図のように連立方程式を行列で手計算で解く時に使われる方法です。
直説法だとコンピューターはあまり嬉しくない
直説法は人間にとっては早いかもしれないですが、ゼロでない部分が多く発生してしまい、その分計算を高速化することができません。
また、単調な繰り返しを得意とするコンピューターの性格上からしても複雑であまり得意ではないです。
そこで、ヤコビ法をはじめとする反復法では、対象の行列Aを対角行列D、 下三角行列Lと上三角行列Uに分解して考えます。
コード
以下にヤコビ法[Jacobi]とガウス・ザイデル法[Gauss-Seidel]の比較実験した際のコードを記載しますので、参考にしてみてください。
ただ、matlab初心者なので色々間違いあったらコメント欄で教えていただけると嬉しいです。
\left\{
\begin{array}{ll}
10x_1 - x_2 - 2x_3 &= 72 \\
-1x_1 + 10x_2 - 2x_3 &= 83 \\
-1x_1 - x_2 + 5x_3 &= 42 \\
\end{array}
\right.
$x_0$を初期値として。
x_0 = (0, 0, 0)
以下のコードは同じディレクトリに置くと動きます。不安な方はgitをご覧ください。
初期値の設定
A = [10 -1 -2; -1 10 -2; -1 -1 5];
b = [72 83 42];
ヤコビ法[Jacobi]
function [outputX,outputRes] = jacobi(A, b, torelance)
% Jacobi法の計算
% A:n次正則な行列
% b:定数
% torelance:許容誤差
% 初期値（適当な解）x0と残差を設定
x0 = [0 0 0];
residual = 1e20;
D = transpose(diag(A)); % D := 対角行列
R = A - diag(D); % R := L + U
i = 0;
disp('----------Start iteration----------');
res = [];
% 反復計算の結果、誤差が許容誤差（torelance）よりも小さくなったら終了
while residual > torelance
% 解と残差を出力
% A./B は要素が A(i,j)/B(i,j) の行列を求めるらしい
x = (b - transpose(R*transpose(x0))) ./ D;
% xがbでx0がAxのつもりで、どれだけ解が欲しいターゲットの値と乖離しているかを計算
denominator = sqrt((x - x0).^2);
numerator = sqrt(x.^2);
residual = sum(denominator) / sum(numerator);
res = [res residual];
% 100回ごとに誤差を出力
if rem(i, 100) == 0
fprintf('Iteration = %d \n', i);
fprintf('Residual = %d \n', residual);
fprintf('x = %d \n', x);
end
% 更新した結果を変数に格納して次の処理に持ち越す。
i = i + 1;
x0 = x;
end
% 終了時の結果を出力
disp('----------End iteration----------')
fprintf('Iteration = %d \n', i);
fprintf('Residual = %d \n', residual);
fprintf('x = %d \n', x0);
outputX = x0;
outputRes = res;
end
ガウス・ザイデル法[Gauss-Seidel]
ヤコビ法では問題なかったのですが、whileの中にforがある二重ループになっています。前の計算を待たないと処理できないので、どうにもならないことなのかなと思いますが、実装していて、アルゴリズムで早くしたのに、計算量で特定の問題ではヤコビ法の方が早かったりするのかなと思いました。
例えば、変数がやたらめったら多くて、手数がかかるような問題では、ガウス法は二重のループがボトルネックになって計算が遅くなってしまう可能性があるのかなー？と
もしご存知の方いたら教えてください。
function [outputX,outputRes] = gauss_seidel(A, b, torelance)
% gauss_seidel法の計算
% A:n次正則な行列
% b:定数
% torelance:許容誤差
% 初期値（適当な解）x0と残差を設定
x = [0 0 0];
residual = 1e20;
D = transpose(diag(A)); % D := 対角行列
R = A - diag(D); % R := L + U
iteration = 0;
disp('----------Start iteration----------');
res = [];
% 反復計算の結果、誤差が許容誤差（torelance）よりも小さくなったら終了
while residual > torelance
% 収束判定用に1ステップ前の解を保持
x0 = x;
for i=1:length(x)
% R(i, :)でi番目の行ごと取り出す
x(i) = (b(i) - transpose(R(i, :)*transpose(x))) ./ D(i);
end
% xがbでx0がAxのつもりで、どれだけ解が欲しいターゲットの値と乖離しているかを計算
denominator = sqrt((x - x0).^2);
numerator = sqrt(x.^2);
residual = sum(denominator) / sum(numerator);
res = [res residual];
% 100回ごとに誤差を出力
if rem(iteration, 100) == 0
fprintf('Iteration = %d \n', iteration);
fprintf('Residual = %d \n', residual);
fprintf('x = %d \n', x);
end
% 更新した結果を変数に格納して次の処理に持ち越す。
iteration = iteration + 1;
end
% 終了時の結果を出力
disp('----------End iteration----------')
fprintf('Iteration = %d \n', iteration);
fprintf('Residual = %d \n', residual);
fprintf('x = %d \n', x0);
outputX = x0;
outputRes = res;
end
グラフの描画
A = [10 -1 -2; -1 10 -2; -1 -1 5];
b = [72 83 42];
[x_jacobi, res_jacobi] = jacobi(A, b, 1e-6)
[x_gauss, res_gauss] = gauss_seidel(A, b, 1e-6)
% プロット
x1 = 1:length(res_jacobi);
x2 = 1:length(res_gauss);
set(gca, 'YScale', 'log')
hold on
plot(x1, res_jacobi, 'LineWidth', 2, 'DisplayName', 'jacobi');
plot(x2, res_gauss, 'LineWidth', 2, 'DisplayName', 'gauss-seidel');
hold off
lgd = legend;
lgd.FontSize = 12;
参考
pythonのをmatlabに書き換えただけです。