2022年度プログラミング演習A・B

第6回レポート課題の解説

演習1

ガウスの消去法における前進消去のアルゴリズムを実現するコードを書き、

$$ A= \begin{pmatrix} 2 & 1 & -2 & 3\\ 2 & 0 & 4 & 1\\ 3 & 0 & 5 & 2\\ 1 & -1 & 2 & 1 \end{pmatrix} ,\quad b= \begin{pmatrix} 11\\ -5\\ -5\\ 4 \end{pmatrix} $$

に対する変形の結果を表示してください。ただし、$A$ が $4$ 次以外の正方行列の場合にも正しく動作するコードにすること。

A = [2,1,-2,3; 2,0,4,1; 3,0,5,2; 1,-1,2,1];
b = [11;-5;-5;4];
n = size(A,1);  %4ではなく変数nを使えば、4次以外の正方行列の場合にも対応できる

%前進消去
for k = 1:n-1  %対角成分の番号
    for i = k+1:n  %行の番号
        alpha = -A(i,k)/A(k,k);        %Aの(i,k)成分を0にするためのスカラー
        A(i,:) = A(i,:)+alpha*A(k,:);  %Aの第i行に第k行のalpha倍を加える
        b(i) = b(i)+alpha*b(k);        %bの第i成分に第k成分のalpha倍を加える
    end
end

A
b

A =

   2.00000   1.00000  -2.00000   3.00000
   0.00000  -1.00000   6.00000  -2.00000
   0.00000   0.00000  -1.00000   0.50000
   0.00000   0.00000   0.00000  -0.50000

b =

   11.0000
  -16.0000
    2.5000
    7.5000

演習2

演習1のコードにガウスの消去法における後退代入のアルゴリズムを実現するコードを追加し、

$$ A= \begin{pmatrix} 2 & 1 & -2 & 3\\ 2 & 0 & 4 & 1\\ 3 & 0 & 5 & 2\\ 1 & -1 & 2 & 1 \end{pmatrix} ,\quad b= \begin{pmatrix} 11\\ -5\\ -5\\ 4 \end{pmatrix} $$

に対して $Ax=b$ の厳密解を求めてください。ただし、$A$ が $4$ 次以外の正方行列の場合にも正しく動作するコードにすること。

A = [2,1,-2,3; 2,0,4,1; 3,0,5,2; 1,-1,2,1];
b = [11;-5;-5;4];
n = size(A,1);

%前進消去
for k = 1:n-1
    for i = k+1:n
        alpha = -A(i,k)/A(k,k);
        A(i,:) = A(i,:)+alpha*A(k,:);
        b(i) = b(i)+alpha*b(k);
    end
end

%後退代入
x = zeros(n,1);  %空の配列でもよいが、他の多くのプログラミング言語では次行でエラーが発生する
x(n) = b(n)/A(n,n);  %実はi=nの場合として次のfor文に含めてもよい（内側のfor文で変数の範囲が空の配列だと何も起こらない）
for i = n-1:-1:1  %n-1から1まで1ずつ減らす
    tmp = b(i);
    for j = i+1:n
        tmp = tmp-A(i,j)*x(j);
    end
    x(i) = tmp/A(i,i);
end

x

x =

   25
  -14
  -10
  -15

（別解）

後退代入の際の分子に現れる積和の形は、ベクトルの内積（横ベクトルと縦ベクトルの行列としての積）を利用すると簡単に計算できます。

A = [2,1,-2,3; 2,0,4,1; 3,0,5,2; 1,-1,2,1];
b = [11;-5;-5;4];
n = size(A,1);

%前進消去
for k = 1:n-1
    for i = k+1:n
        alpha = -A(i,k)/A(k,k);
        A(i,:) = A(i,:)+alpha*A(k,:);
        b(i) = b(i)+alpha*b(k);
    end
end

%後退代入
x = zeros(n,1);
x(n) = b(n)/A(n,n);
for i = n-1:-1:1
    x(i) = (b(i)-A(i,i+1:n)*x(i+1:n))/A(i,i);  %内積の計算にはdot関数を使ってもよい
end

x

x =

   25
  -14
  -10
  -15

演習3（オプション）

ピボット選択を行うガウスの消去法のアルゴリズムを実現するコードを書き、ピボット選択が必要な $A$ と $b$ の例を与えて $Ax=b$ の厳密解を求めてください。

ピボット選択が必要な例として

$$ A= \begin{pmatrix} 0 & 1 & -2 & 3\\ 1 & 2 & 4 & 3\\ 2 & 0 & 5 & 2\\ 4 & 8 & 0 & 4 \end{pmatrix} ,\quad b= \begin{pmatrix} 1\\ 1\\ 1\\ 1 \end{pmatrix} $$

とすると、ピボット選択を行わない上記のコードでは次に示すように解を求めることができません。

A = [0,1,-2,3; 1,2,4,3; 2,0,5,2; 4,8,0,4];
b = [1;1;1;1];
n = size(A,1);

%前進消去
for k = 1:n-1
    for i = k+1:n
        alpha = -A(i,k)/A(k,k);
        A(i,:) = A(i,:)+alpha*A(k,:);
        b(i) = b(i)+alpha*b(k);
    end
end

A
b

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

x

warning: division by zero
warning: division by zero
warning: division by zero
A =

     0     1    -2     3
   NaN  -Inf   Inf  -Inf
   NaN   NaN   NaN   NaN
   NaN   NaN   NaN   NaN

b =

     1
  -Inf
   NaN
   NaN

warning: division by zero
x =

   NaN
   NaN
   NaN
   NaN

ピボット選択を行うガウスの消去法のコードは、次のようになります。

A = [0,1,-2,3; 1,2,4,3; 2,0,5,2; 4,8,0,4];
b = [1;1;1;1];
n = size(A,1);

%前進消去
for k = 1:n-1
    %ピボット選択
    if abs(A(k,k))<1E-15  %浮動小数点数の計算誤差を考慮して「== 0」の判定は避けるべき
        %第k列の対角より下で絶対値が最大の成分の行番号を求める
        l = k+1;  %求める行番号の候補
        max_value = abs(A(l,k));  %最大値の候補
        for i = k+2:n
            if abs(A(i,k)) > max_value  %最大値の候補より大きなものが見つかったら候補を更新する
                l = i;
                max_value = abs(A(l,k));
            end
        end
        %Aの第k行と第l行、bの第k成分と第l成分を入れ替える
        tmp = A(k,:);  %一時変数を用いてそれぞれの入れ替えを行う
        A(k,:) = A(l,:);
        A(l,:) = tmp;
        tmp = b(k);
        b(k) = b(l);
        b(l) = tmp;
    end
    for i = k+1:n
        alpha = -A(i,k)/A(k,k);
        A(i,:) = A(i,:)+alpha*A(k,:);
        b(i) = b(i)+alpha*b(k);
    end
end

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

x

x =

   0.12500
  -0.12500
   0.00000
   0.37500

（別解）

ピボット選択の処理をより簡潔に書くと、次のようになります。

A = [0,1,-2,3; 1,2,4,3; 2,0,5,2; 4,8,0,4];
b = [1;1;1;1];
n = size(A,1);

%前進消去
for k = 1:n-1
    %ピボット選択
    if abs(A(k,k))<1E-15
        %第k列の対角より下で絶対値が最大の成分の行番号を求める
        [max_value,l] = max(abs(A(k+1:n,k)));  %max関数の特別な使い方で、最大値とそれをとる番号を変数に格納する
        l = l+k;  %行番号に直す
        %Aの第k行と第l行、bの第k成分と第l成分を入れ替える
        A([k,l],:) = A([l,k],:);  %Aの第k行と第l行でできる行列に、Aの第l行と第k行でできる行列を代入する
        b([k,l]) = b([l,k]);  %bの第k成分と第l成分でできるベクトルに、bの第l成分と第k成分でできるベクトルを代入する
    end
    for i = k+1:n
        alpha = -A(i,k)/A(k,k);
        A(i,:) = A(i,:)+alpha*A(k,:);
        b(i) = b(i)+alpha*b(k);
    end
end

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

x

x =

   0.12500
  -0.12500
   0.00000
   0.37500

（参考）

ガウスの消去法の前進消去では、注目する対角成分が $0$ かどうかに関わらず常にピボット選択を行う方が、得られる解の計算誤差が大抵の場合は少なくなることが知られています。そのため、そちらの方が一般的です。

ただし、それは行の入れ替え回数が増えることを意味し、行列の次数が大きい場合には特にメモリ上での値の入れ替えを非常に多く行うことになります。

この非効率への対策として、行列 $A$ の行およびベクトル $b$ の成分をメモリ上で実際には入れ替えずに、配列を用意してその成分を入れ替えることで本来の行番号とメモリ上の行番号の対応関係を管理するという方法があります。

$Ax=b$ が解けない場合の対応もついでに含めて上記を実装すると、次のようになります（複雑なので、無理に理解しなくて大丈夫です）。

A = [0,1,-2,3; 1,2,4,3; 2,0,5,2; 4,8,0,4];
b = [1;1;1;1];
n = size(A,1);

index = 1:n;  %行番号の対応関係を管理する配列（置換写像）
%Aの行およびbの成分を入れ替える代わりに、この対応する成分を入れ替える
%本来の行番号iに対してメモリ上の行番号はindex(i)になる
%Aの行およびbの成分にアクセスする際は、この配列を通してアクセスする

error_flag = 0;  %エラー判定のためのフラグ（Ax=bが解けない場合は1にする）

%前進消去
for k = 1:n-1
    %ピボット選択
    %第k列の対角以下で絶対値が最大の成分の行番号を求める
    [max_value,l] = max(abs(A(index(k:n),k)));
    l = l+k-1;
    if max_value < 1E-15  %絶対値が最大の成分がほぼ0ならエラー
        error_flag = 1;
        break;  %for文を終了する
    end
    if l ~= k
        %Aの第k行と第l行、bの第k成分と第l成分を入れ替える代わりに、indexの対応する成分を入れ替える
        index([k,l])=index([l,k]);
    end
    for i = k+1:n
        alpha = -A(index(i),k)/A(index(k),k);
        A(index(i),:) = A(index(i),:)+alpha*A(index(k),:);
        b(index(i)) = b(index(i))+alpha*b(index(k));
    end
end

if abs(A(index(n),n)) < 1E-15  %Aの(n,n)成分がほぼ0ならエラー
    error_flag = 1;
end

if error_flag == 1
    display("エラー：Ax=bは唯一の解をもたない。")
else
    %後退代入
    x = zeros(n,1);
    x(n) = b(index(n))/A(index(n),n);
    for i = n-1:-1:1
        x(i) = (b(index(i))-A(index(i),i+1:n)*x(i+1:n))/A(index(i),i);
    end
    x
end

x =

   0.12500
  -0.12500
   0.00000
   0.37500

なお、今回扱ったのはピボット選択の中でも行の入れ替えのみを伴う部分ピボット選択と呼ばれる方法であり、他に列の入れ替えも伴う完全ピボット選択と呼ばれる方法があります。