ガウス消去法(II) LU分解

ガウス消去法によって、行列のLU分解ができます。以下、$4\times 4$の行列$A$を例にして説明します。

Aの第i列の対角成分以下の要素を消去するために、以下の作業行列$E_i$が用意されます。

$$ E_1= \left( \begin{array}{cccc} 1 & 0 & 0 & 0\\ f_{21} & 1 & 0 & 0 \\ f_{31} & 0 & 1 & 0\\ f_{41} & 0 & 0 & 1 \end{array} \right) $$$$ E_2= \left( \begin{array}{cccc} 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 \\ 0 & f_{32} & 1 & 0\\ 0 & f_{42} & 0 & 1 \end{array} \right) $$$$ E_3= \left( \begin{array}{cccc} 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0\\ 0 & 0& f_{43} & 1 \end{array} \right) $$

ガウス消去法によって、以下の分解ができます。 $$ E_3E_2E_1 A = U $$ よって、 $$ A = E_1^{-1} E_2^{-1} E_3^{-1} U = LU \quad (L := E_1^{-1} E_2^{-1} E_3^{-1} ) $$

Lの各成分を計算すると、以下の形が分かります。

$$ L= \left( \begin{array}{cccc} 1 & 0 & 0 & 0\\ -f_{21} & 1 & 0 & 0 \\ -f_{31} & -f_{32} & 1 & 0\\ -f_{41} & -f_{42} & -f_{43} & 1 \end{array} \right) $$

上記の計算過程によって、ガウス消去法の各列の作業とともに、Lの各列の成分を計算でき、Aの対角成分以下に格納することができます。 即ち、ガウス消去法の計算後、以下の行列が得られます。

$$ \left( \begin{array}{cccc} u_{11} & u_{12} & u_{13} & u_{14}\\ -f_{21} & u_{22} & u_{23} & u_{24} \\ -f_{31} & -f_{32} & u_{33} & u_{34}\\ -f_{41} & -f_{42} & -f_{43} & u_{44} \end{array} \right) $$

演習１

行列のLU分解を行い、上記の説明のようにLとUを求めてください。

$$ A= \left( \begin{array}{cccc} 2& 1 & -2 & 3\\ 2 & 0 & 4 & 1 \\ 3 & 0 & 5& 2\\ 1 & -1 & 2 & 1 \end{array} \right) $$

Step 1.

Aの行列の第１行を使って、行列の一番目の列の対角以下の成分を０に変形します。

まず、第１行のalpha倍を第2,3,4行に加えます。 作業行列$E_1$の成分 $f_{j1}$ (j=2,3,4)$をAの第１列格納する。

In [7]:
A = [2, 1, -2, 3; 2, 0, 4, 1; 3, 0, 5, 2; 1, -1, 2, 1];
n = 4;
%対角の第１成分からの作業
i = 1;

%----------- 第１行によるAの第１列の対角以下の成分を消去する ----------- 
for j=2:n
%作業の行の番号
alpha = -A(j,1)/A(1,1);
for k=2:n %ｋ：各列の番号
  A(j,k) = A(j,k) + alpha*A(1,k);
end
%作業行列のf_ｊ1をAの第１列格納する。
A(j,1) = - alpha ;
end

%変形後のA
A

A =

   2.00000   1.00000  -2.00000   3.00000
   1.00000  -1.00000   6.00000  -2.00000
   1.50000  -1.50000   8.00000  -2.50000
   0.50000  -1.50000   3.00000  -0.50000

Step 2.

Aの行列の第2行を使って、2番目の列の対角以下の成分を０に変形します。

第2行のalpha倍を第3,4行に加えます。 作業行列$E_2$の成分 $f_{j2}$ (j=3,4)$をAの第2列の対角以下に格納する。

In [ ]:
i=2;

%----------- 第２行によるAの第２列の対角以下の成分を消去する ----------- 

for j=2:n
%作業の行の番号
alpha = ??;
for k=??:n %ｋ：各列の番号
  A(j,k) = A(j,k) + ??;
end
%作業行列のf_ｊ2をAの第2列に格納する。
A(??,??) = ?? ;
end

%変形後のA
A

Step 3.

Aの行列の第3行を使って、3番目の列の対角以下の成分を０に変形します。

第2行のalpha倍を第4行に加えます。 作業行列$E_3$の成分 $f_{43}$をAの第3列の対角以下に格納する。

In [9]:
i=3;

%----------- 第3行によるAの第3列の対角以下の成分を消去する ----------- 

%ここにコードを書いてください。

コードのまとめ

上記の計算をここにまとめて書いてください。行の番号iに関するループを使ってください。

In [10]:
%ここにコードを書いてください。

LU分解による連立一次方程式を解く

AのLU分解が分かると、以下の２つの三角行列に関する方程式を解いてから、Ax=bの解は容易に計算できます。

$$ Ly=b, Ux=y $$

$Ly=b$の計算方法

Lは以下の形を持っています。 $$ L= \left( \begin{array}{cccc} 1 & 0 & 0 & 0\\ l_{21} & 1 & 0 & 0 \\ l_{31} & l_{32} & 1 & 0\\ l_{41} & l_{42} & l_{43} & 1 \end{array} \right) $$

よって、$Ly=b$の解$y=(y_1, y_2, y_3, y_4)$は以下のように$y_i$を順番に算出することができます。

  • $y_1 = b_1$
  • $y_1 l_{21} + y_2 = b_2 \to y_2 = b_2 - y_1 l_{21}$
  • $y_1 l_{31} + y_2 l_{32}+y_3 = b_3 \to y_3 = b_3 - y_1 l_{31} - y_2 l_{32}$
  • ...

for loopでまとめます。

In [ ]:
for i=1:n
   y(i) = b(i)
   for j = 1:(i-1)
     y(i) = y(i) - y(j)*L(i,j)
   end
end

講義番号：１８１S1520　　水準　０４

演習２

上記のAのLU分解に利用して、以下の$b_1$, $b_2$に関する連立一次方程式$Ax=b_1, Ax=b_2$を解いてください。 $$ b_1= \left( \begin{array}{cccc} 11\\ -5\\ -5\\ -4 \end{array} \right),\quad b_2= \left( \begin{array}{cccc} 10\\ 18\\ 26\\ 9 \end{array} \right) $$

演習３

以下の行列$A_1$,$A_2$のLU分解を計算してください。

  • $A_1$は以下のように生成されます。

    A1=diag(2*ones(10,1))+diag(ones(9,1),1)+diag(ones(9,1),-1)

  • $A_2$は$6\times 6$の行列であり、第i行、第j列の成分は $$ A_2(i,j) = \frac{1}{i+j-1} $$ です。
In [16]:
A1=diag(2*ones(10,1))+diag(ones(9,1),1)+diag(ones(9,1),-1)

A1 =

   2   1   0   0   0   0   0   0   0   0
   1   2   1   0   0   0   0   0   0   0
   0   1   2   1   0   0   0   0   0   0
   0   0   1   2   1   0   0   0   0   0
   0   0   0   1   2   1   0   0   0   0
   0   0   0   0   1   2   1   0   0   0
   0   0   0   0   0   1   2   1   0   0
   0   0   0   0   0   0   1   2   1   0
   0   0   0   0   0   0   0   1   2   1
   0   0   0   0   0   0   0   0   1   2

演習4

$A_2 x=b$を解いてください。ただし、$b$はすべて成分が１である$6\times 1$ベクトルです。

In [ ]: