next up previous
Next: 4.2 課題 Up: 4 Octave, Scilab 入門および本日の課題 Previous: 4 Octave, Scilab 入門および本日の課題

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

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

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

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

について調べてみよう。


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%


next up previous
Next: 4.2 課題 Up: 4 Octave, Scilab 入門および本日の課題 Previous: 4 Octave, Scilab 入門および本日の課題
Masashi Katsurada
平成15年5月8日