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

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

演習1

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

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

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

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

%前進消去によるLU分解
for k = 1:n-1  %対角成分の番号
    for i = k+1:n  %行の番号
        alpha = -A(i,k)/A(k,k);                    %Aの(i,k)成分を0にするためのスカラー
        A(i,k+1:n) = A(i,k+1:n)+alpha*A(k,k+1:n);  %Aの第i行に第k行のalpha倍を加える（第1列～第k列の成分の計算はしない）
        A(i,k) = -alpha;                           %下三角行列の成分を格納する
    end
end

A

A =

   2.00000   1.00000  -2.00000   3.00000
   1.00000  -1.00000   6.00000  -2.00000
   1.50000   1.50000  -1.00000   0.50000
   0.50000   1.50000   6.00000  -0.50000

よって、与えられた行列 $A$ は次の下三角行列 $L$ と上三角行列 $U$ を用いて $A=LU$ と分解できることが分かる。

$$ L= \begin{pmatrix} 1 & 0 & 0 & 0\\ 1 & 1 & 0 & 0\\ 1.5 & 1.5 & 1 & 0\\ 0.5 & 1.5 & 6 & 1 \end{pmatrix} ,\quad U= \begin{pmatrix} 2 & 1 & -2 & 3\\ 0 & -1 & 6 & -2\\ 0 & 0 & -1 & 0.5\\ 0 & 0 & 0 & -0.5 \end{pmatrix} $$

演習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$ と $b$ の次数が $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);

%前進消去によるLU分解
for k = 1:n-1
    for i = k+1:n
        alpha = -A(i,k)/A(k,k);
        A(i,k+1:n) = A(i,k+1:n)+alpha*A(k,k+1:n);
        A(i,k) = -alpha;
    end
end

%Aを使ってLとUを作成する必要はない

%前進代入
%L（Aの対角より下）とbを用いて、Ly=bの解yを求める
y = zeros(n,1);
y(1) = b(1);
for i = 2:n
    tmp = b(i);
    for j = 1:i-1
        tmp = tmp-A(i,j)*y(j);
    end
    y(i) = tmp;
end

%後退代入
%U（Aの対角以上）とyを用いて、Ux=yの解xを求める
x = zeros(n,1);
x(n) = y(n)/A(n,n);
for i = n-1:-1:1
    tmp = y(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);

%前進消去によるLU分解
for k = 1:n-1
    for i = k+1:n
        alpha = -A(i,k)/A(k,k);
        A(i,k+1:n) = A(i,k+1:n)+alpha*A(k,k+1:n);
        A(i,k) = -alpha;
    end
end

%前進代入
%L（Aの対角より下）とbを用いて、Ly=bの解yを求める
y = zeros(n,1);
y(1) = b(1);
for i = 2:n
    y(i) = b(i)-A(i,1:i-1)*y(1:i-1);  %内積（行列積）を利用
end

%後退代入
%U（Aの対角以上）とyを用いて、Ux=yの解xを求める
x = zeros(n,1);
x(n) = y(n)/A(n,n);
for i = n-1:-1:1
    x(i) = (y(i)-A(i,i+1:n)*x(i+1:n))/A(i,i);  %内積（行列積）を利用
end

x

x =

   25
  -14
  -10
  -15