4.1 三重対角行列の逆行列と LU 分解

行列 $A$ が疎行列であっても、その逆行列 $A^{-1}$ は疎行列とは限らない。 しかし $P A=L U$ と LU 分解したときの¹¹、$L$, $U$ は疎行列である。

このことを $3$ 重対角行列

$$A= \left( \begin{array}{rrrrrr} 2 & -1 \\ -1 & 2 & -1\... ...ots \\ & & -1 & 2 & -1\\ & & & -1 & 2 \end{array} \right)$$

について調べてみよう。

mathpc% octave
GNU Octave, version 2.0.16 (i386-unknown-freebsd3.4).
Copyright (C) 1996, 1997, 1998, 1999, 2000 John W. Eaton.
This is free software with ABSOLUTELY NO WARRANTY.
For details, type `warranty'.

octave:1> n=5
n = 5
octave:2> ones(n-1,1)
ans =

  1
  1
  1
  1

octave:3> J=diag(ones(n-1,1),1)
J =

  0  1  0  0  0
  0  0  1  0  0
  0  0  0  1  0
  0  0  0  0  1
  0  0  0  0  0

octave:4> J+J'
ans =

  0  1  0  0  0
  1  0  1  0  0
  0  1  0  1  0
  0  0  1  0  1
  0  0  0  1  0

octave:5> a=2*eye(n,n)-J-J'
a =

   2  -1   0   0   0
  -1   2  -1   0   0
   0  -1   2  -1   0
   0   0  -1   2  -1
   0   0   0  -1   2

octave:6> inv(a)
ans =

  0.83333  0.66667  0.50000  0.33333  0.16667
  0.66667  1.33333  1.00000  0.66667  0.33333
  0.50000  1.00000  1.50000  1.00000  0.50000
  0.33333  0.66667  1.00000  1.33333  0.66667
  0.16667  0.33333  0.50000  0.66667  0.83333

octave:7> [L U P]=lu(a)
L =

   1.00000   0.00000   0.00000   0.00000   0.00000
  -0.50000   1.00000   0.00000   0.00000   0.00000
   0.00000  -0.66667   1.00000   0.00000   0.00000
   0.00000   0.00000  -0.75000   1.00000   0.00000
   0.00000   0.00000   0.00000  -0.80000   1.00000

U =

   2.00000  -1.00000   0.00000   0.00000   0.00000
   0.00000   1.50000  -1.00000   0.00000   0.00000
   0.00000   0.00000   1.33333  -1.00000   0.00000
   0.00000   0.00000   0.00000   1.25000  -1.00000
   0.00000   0.00000   0.00000   0.00000   1.20000

P =

  1  0  0  0  0
  0  1  0  0  0
  0  0  1  0  0
  0  0  0  1  0
  0  0  0  0  1

octave:8> U\(L\P)        あるいは  inv(U)*inv(L)*P
ans =

  0.83333  0.66667  0.50000  0.33333  0.16667
  0.66667  1.33333  1.00000  0.66667  0.33333
  0.50000  1.00000  1.50000  1.00000  0.50000
  0.33333  0.66667  1.00000  1.33333  0.66667
  0.16667  0.33333  0.50000  0.66667  0.83333

octave:9> n=8;J=diag(ones(n-1,1),1);a=2*eye(n,n)-J-J';[L U P]=lu(a);inv(a)
ans =

  0.88889  0.77778  0.66667  0.55556  0.44444  0.33333  0.22222  0.11111
  0.77778  1.55556  1.33333  1.11111  0.88889  0.66667  0.44444  0.22222
  0.66667  1.33333  2.00000  1.66667  1.33333  1.00000  0.66667  0.33333
  0.55556  1.11111  1.66667  2.22222  1.77778  1.33333  0.88889  0.44444
  0.44444  0.88889  1.33333  1.77778  2.22222  1.66667  1.11111  0.55556
  0.33333  0.66667  1.00000  1.33333  1.66667  2.00000  1.33333  0.66667
  0.22222  0.44444  0.66667  0.88889  1.11111  1.33333  1.55556  0.77778
  0.11111  0.22222  0.33333  0.44444  0.55556  0.66667  0.77778  0.88889

octave:10> U\(L\P)     あるいは   inv(U)*inv(L)*P
ans =

  0.88889  0.77778  0.66667  0.55556  0.44444  0.33333  0.22222  0.11111
  0.77778  1.55556  1.33333  1.11111  0.88889  0.66667  0.44444  0.22222
  0.66667  1.33333  2.00000  1.66667  1.33333  1.00000  0.66667  0.33333
  0.55556  1.11111  1.66667  2.22222  1.77778  1.33333  0.88889  0.44444
  0.44444  0.88889  1.33333  1.77778  2.22222  1.66667  1.11111  0.55556
  0.33333  0.66667  1.00000  1.33333  1.66667  2.00000  1.33333  0.66667
  0.22222  0.44444  0.66667  0.88889  1.11111  1.33333  1.55556  0.77778
  0.11111  0.22222  0.33333  0.44444  0.55556  0.66667  0.77778  0.88889

octave:11> quit
mathpc%