物理のノート
逆行列
逆行列の数値計算($2\times2$正方行列)
逆行列$A^{-1}$は,行列式と
余因子(Cofactor)を用いて,
次のように定義される.
次の$2\times2$正方行列の場合,
\begin{equation}
A=
\begin{pmatrix}
a & b\\
c & d
\end{pmatrix}
,
\end{equation}
逆行列$A^{-1}$は次のようになる.
\begin{equation}
A^{-1}
=\frac{1}{\lvert A \rvert }
\begin{pmatrix}
C_{1,1} & C_{2,1}\\
C_{1,2} & C_{2,2}
\end{pmatrix}
=\frac{1}{\lvert A \rvert }
\begin{pmatrix}
(-1)^{1+1}\Delta_{1,1} & (-1)^{2+1}\Delta_{2,1}\\
(-1)^{1+2}\Delta_{1,2} & (-1)^{2+2}\Delta_{2,2}
\end{pmatrix}
=
\frac{1}{ad-bc }
\begin{pmatrix}
d & (-1)b\\
(-1)c & a
\end{pmatrix}\\
=
\frac{1}{ad-bc }
\begin{pmatrix}
d & -b\\
-c & a
\end{pmatrix}
\end{equation}
この定義が示すように,逆行列が存在するには, 行列式$\lvert A \rvert \ne 0$が条件となる.
定義に従って,次の$2\times 2$正方行列の逆行列を求める.
\begin{equation}
A=
\begin{pmatrix}
1 & 3\\
2 & 4
\end{pmatrix}
\end{equation}
とすると,
\begin{equation}
A^{-1}=\frac{1}{1\cdot 4-3\cdot 2}
\begin{pmatrix}
4 & (-1)3\\
(-1)2 & 1
\end{pmatrix}\\
=-\frac{1}{2}
\begin{pmatrix}
4 & -3\\
-2 & 1
\end{pmatrix}
=
\begin{pmatrix}
-2 & \frac{3}{2}\\
1 & -\frac{1}{2}
\end{pmatrix}
\end{equation}
ところで,この$2\times 2$正方行列の逆行列は,
Gauss-Jordan消去法(掃き出し法)で用いた数値計算を次のように改良して求めることができる.
=============
2 * 4 array
=============
array2 =
[[ 1. 3. 1. 0.]
[ 2. 4. 0. 1.]]
=====================================
Gauss-Jordan Elimination(掃き出し法)
=====================================
array2 =
[[ 2. 2. 0. 1.]
[ 1. 3. 1. 0.]]
array2 =
[[ 2. 2. 0. 1.]
[ 1. 3. 1. 0.]]
array2 =
[[ 2. 2. 0. 0.5]
[ 1. 3. 1. 0. ]]
array2 =
[[ 1. 2. 0. 0.5]
[ 0. 3. 1. 0. ]]
array2 =
[[ 1. 2. 0. 0.5]
[ 0. 1. 1. 0. ]]
array2 =
[[ 1. 2. 0. 0.5]
[ 0. 1. 1. 0. ]]
array2 =
[[ 1. 2. 0. 0.5]
[ 0. 1. 1. -0.5]]
array2 =
[[ 1. 2. 0. 0.5]
[ 0. 1. 1. -0.5]]
array2 =
[[ 1. 2. 0. 0.5]
[ 0. 1. 1. -0.5]]
array2 =
[[ 1. 0. 0. 0.5]
[ 0. 1. 1. -0.5]]
array2 =
[[ 1. 0. -2. 0.5]
[ 0. 1. 1. -0.5]]
array2 =
[[ 1. 0. -2. 1.5]
[ 0. 1. 1. -0.5]]
=====================
determinant of matrix
=====================
det = 2.0
Pythonで実装したアルゴリズムは次の通りである.
#!/usr/bin/python2.7
# -*- coding: utf-8 -*-
import numpy as np
print('=============')
print('2 * 2 array')
print('=============')
array2 = np.array([[1.,3.,1.,0.],[2.,4.,0.,1.]])
arrayp = np.array([[0.,0.],[0.,0.]])
arrayq = np.array([[0.,0.],[0.,0.]])
tmp = np.array([[0.,0.,0.,0.],[0.,0.,0.,0.]])
print('array2 =')
print(array2)
print('=====================================')
print('Gauss-Jordan Elimination(掃き出し法)')
print('=====================================')
#n = 2 元の数
det = 1.0
for k in range(0,2,1):
for i in range(k,2,1):
if i !=k:
if abs(array2[k,k]) < abs(array2[i,k]):
# j列の範囲をk列目から3列目に変更する.
for j in range(k,4,1):
tmp[i,j] = array2[k,j]
array2[k,j] = array2[i,j]
array2[i,j] = tmp[i,j]
arrayp[k,k] = array2[k,k]
det = det * arrayp[k,k]
# j列の範囲をk列目から3列目に変更する.
for j in range(k,3,1):
array2[k,j+1] = array2[k,j+1] / arrayp[k,k]
print('array2 =')
print(array2)
array2[k,k] = 1
for i in range(0,2,1):
arrayq[k,k] = array2[i,k]
# j列の範囲をk列目から3列目に変更する.
for j in range(k,4,1):
if i != k:
array2[i,j] = array2[i,j] - array2[k,j] * arrayq[k,k]
print('array2 =')
print(array2)
else:
array2[i,j] = array2[i,j]
##行列式の計算結果
print('=====================')
print('determinant of matrix')
print('=====================')
print('det = {}'.format(det))
このアルゴリズムでは,行列$A$の列を右に拡張し,$2\times2$単位行列を加えている.
この状態で,掃き出し法(Gauss-Jordan消去法)を実行すると,実行結果には$2\times2$単位行列と逆行列$A^{-1}$が並んだ$2\times4$行列が前述のように得られる.
念のため,検算をしてみよう.
ある行列とその逆行列の積は単位行列になる.
\[
AA^{-1}=I
\]
行列同士の積の定義に沿ってコードを実装する.
#!/usr/bin/python2.7
# -*- coding: utf-8 -*-
import numpy as np
print('==========================')
print('2 * 2 array multiplication')
print('==========================')
a = np.array([[1.,3.],[2.,4.]])
b = np.array([[-2.,1.5],[1.,-0.5]])
c = np.array([[0.,0.],[0.,0.]])
#行列Aをaとする.
print('a =')
print(a)
#行列Aの逆行列をbとする.
print('b =')
print(b)
for i in range(0,2,1):
for j in range(0,2,1):
for k in range(0,2,1):
c[i,j] += a[i,k]*b[k,j]
#行列aと行列bの乗算結果を行列cとする.
print('c =')
print(c)
実行結果は次の通り.
==========================
2 * 2 array multiplication
==========================
a =
[[ 1. 3.]
[ 2. 4.]]
b =
[[-2. 1.5]
[ 1. -0.5]]
c =
[[ 1. 0.]
[ 0. 1.]]
単位行列が得られた.
逆行列の数値計算($3\times 3$正方行列)
$3\times3$正方行列についても同様に,逆行列を求めてみる.
行列式の数値計算で利用した次の行列を用いる.
\begin{equation}
\begin{pmatrix}
3 & 1 & 2 \\
5 & 1 & 3 \\
4 & 2 & 1 \\
\end{pmatrix}
\end{equation}
結果は次のようになる.
[[ 1. 0. 0. -1.25 0.75 0.25]
[ 0. 1. 0. 1.75 -1.25 0.25]
[ 0. 0. 1. 1.5 -0.5 -0.5 ]]
途中の計算プロセスと行列式を記す.
=============
3 * 4 array
=============
array2 =
[[ 3. 1. 2. 1. 0. 0.]
[ 5. 1. 3. 0. 1. 0.]
[ 4. 2. 1. 0. 0. 1.]]
=====================================
Gauss-Jordan Elimination(掃き出し法)
=====================================
array2 =
[[ 5. 0.2 3. 0. 1. 0. ]
[ 3. 1. 2. 1. 0. 0. ]
[ 4. 2. 1. 0. 0. 1. ]]
array2 =
[[ 5. 0.2 0.6 0. 1. 0. ]
[ 3. 1. 2. 1. 0. 0. ]
[ 4. 2. 1. 0. 0. 1. ]]
array2 =
[[ 5. 0.2 0.6 0. 1. 0. ]
[ 3. 1. 2. 1. 0. 0. ]
[ 4. 2. 1. 0. 0. 1. ]]
array2 =
[[ 5. 0.2 0.6 0. 0.2 0. ]
[ 3. 1. 2. 1. 0. 0. ]
[ 4. 2. 1. 0. 0. 1. ]]
array2 =
[[ 5. 0.2 0.6 0. 0.2 0. ]
[ 3. 1. 2. 1. 0. 0. ]
[ 4. 2. 1. 0. 0. 1. ]]
array2 =
[[ 1. 0.2 0.6 0. 0.2 0. ]
[ 0. 1. 2. 1. 0. 0. ]
[ 4. 2. 1. 0. 0. 1. ]]
array2 =
[[ 1. 0.2 0.6 0. 0.2 0. ]
[ 0. 0.4 2. 1. 0. 0. ]
[ 4. 2. 1. 0. 0. 1. ]]
array2 =
[[ 1. 0.2 0.6 0. 0.2 0. ]
[ 0. 0.4 0.2 1. 0. 0. ]
[ 4. 2. 1. 0. 0. 1. ]]
array2 =
[[ 1. 0.2 0.6 0. 0.2 0. ]
[ 0. 0.4 0.2 1. 0. 0. ]
[ 4. 2. 1. 0. 0. 1. ]]
array2 =
[[ 1. 0.2 0.6 0. 0.2 0. ]
[ 0. 0.4 0.2 1. -0.6 0. ]
[ 4. 2. 1. 0. 0. 1. ]]
array2 =
[[ 1. 0.2 0.6 0. 0.2 0. ]
[ 0. 0.4 0.2 1. -0.6 0. ]
[ 4. 2. 1. 0. 0. 1. ]]
array2 =
[[ 1. 0.2 0.6 0. 0.2 0. ]
[ 0. 0.4 0.2 1. -0.6 0. ]
[ 0. 2. 1. 0. 0. 1. ]]
array2 =
[[ 1. 0.2 0.6 0. 0.2 0. ]
[ 0. 0.4 0.2 1. -0.6 0. ]
[ 0. 1.2 1. 0. 0. 1. ]]
array2 =
[[ 1. 0.2 0.6 0. 0.2 0. ]
[ 0. 0.4 0.2 1. -0.6 0. ]
[ 0. 1.2 -1.4 0. 0. 1. ]]
array2 =
[[ 1. 0.2 0.6 0. 0.2 0. ]
[ 0. 0.4 0.2 1. -0.6 0. ]
[ 0. 1.2 -1.4 0. 0. 1. ]]
array2 =
[[ 1. 0.2 0.6 0. 0.2 0. ]
[ 0. 0.4 0.2 1. -0.6 0. ]
[ 0. 1.2 -1.4 0. -0.8 1. ]]
array2 =
[[ 1. 0.2 0.6 0. 0.2 0. ]
[ 0. 0.4 0.2 1. -0.6 0. ]
[ 0. 1.2 -1.4 0. -0.8 1. ]]
array2 =
[[ 1. 0.2 0.6 0. 0.2 0. ]
[ 0. 1.2 -1.16666667 0. -0.8 1. ]
[ 0. 0.4 0.2 1. -0.6 0. ]]
array2 =
[[ 1. 0.2 0.6 0. 0.2 0. ]
[ 0. 1.2 -1.16666667 0. -0.8 1. ]
[ 0. 0.4 0.2 1. -0.6 0. ]]
array2 =
[[ 1. 0.2 0.6 0. 0.2 0. ]
[ 0. 1.2 -1.16666667 0. -0.66666667 1. ]
[ 0. 0.4 0.2 1. -0.6 0. ]]
array2 =
[[ 1. 0.2 0.6 0. 0.2 0. ]
[ 0. 1.2 -1.16666667 0. -0.66666667 0.83333333]
[ 0. 0.4 0.2 1. -0.6 0. ]]
array2 =
[[ 1. 0. 0.6 0. 0.2 0. ]
[ 0. 1. -1.16666667 0. -0.66666667 0.83333333]
[ 0. 0.4 0.2 1. -0.6 0. ]]
array2 =
[[ 1. 0. 0.83333333 0. 0.2 0. ]
[ 0. 1. -1.16666667 0. -0.66666667 0.83333333]
[ 0. 0.4 0.2 1. -0.6 0. ]]
array2 =
[[ 1. 0. 0.83333333 0. 0.2 0. ]
[ 0. 1. -1.16666667 0. -0.66666667 0.83333333]
[ 0. 0.4 0.2 1. -0.6 0. ]]
array2 =
[[ 1. 0. 0.83333333 0. 0.33333333 0. ]
[ 0. 1. -1.16666667 0. -0.66666667 0.83333333]
[ 0. 0.4 0.2 1. -0.6 0. ]]
array2 =
[[ 1. 0. 0.83333333 0. 0.33333333 -0.16666667]
[ 0. 1. -1.16666667 0. -0.66666667 0.83333333]
[ 0. 0.4 0.2 1. -0.6 0. ]]
array2 =
[[ 1. 0. 0.83333333 0. 0.33333333 -0.16666667]
[ 0. 1. -1.16666667 0. -0.66666667 0.83333333]
[ 0. 0. 0.2 1. -0.6 0. ]]
array2 =
[[ 1. 0. 0.83333333 0. 0.33333333 -0.16666667]
[ 0. 1. -1.16666667 0. -0.66666667 0.83333333]
[ 0. 0. 0.66666667 1. -0.6 0. ]]
array2 =
[[ 1. 0. 0.83333333 0. 0.33333333 -0.16666667]
[ 0. 1. -1.16666667 0. -0.66666667 0.83333333]
[ 0. 0. 0.66666667 1. -0.6 0. ]]
array2 =
[[ 1. 0. 0.83333333 0. 0.33333333 -0.16666667]
[ 0. 1. -1.16666667 0. -0.66666667 0.83333333]
[ 0. 0. 0.66666667 1. -0.33333333 0. ]]
array2 =
[[ 1. 0. 0.83333333 0. 0.33333333 -0.16666667]
[ 0. 1. -1.16666667 0. -0.66666667 0.83333333]
[ 0. 0. 0.66666667 1. -0.33333333 -0.33333333]]
array2 =
[[ 1. 0. 0.83333333 0. 0.33333333 -0.16666667]
[ 0. 1. -1.16666667 0. -0.66666667 0.83333333]
[ 0. 0. 0.66666667 1.5 -0.33333333 -0.33333333]]
array2 =
[[ 1. 0. 0.83333333 0. 0.33333333 -0.16666667]
[ 0. 1. -1.16666667 0. -0.66666667 0.83333333]
[ 0. 0. 0.66666667 1.5 -0.5 -0.33333333]]
array2 =
[[ 1. 0. 0.83333333 0. 0.33333333 -0.16666667]
[ 0. 1. -1.16666667 0. -0.66666667 0.83333333]
[ 0. 0. 0.66666667 1.5 -0.5 -0.5 ]]
array2 =
[[ 1. 0. 0. 0. 0.33333333 -0.16666667]
[ 0. 1. -1.16666667 0. -0.66666667 0.83333333]
[ 0. 0. 1. 1.5 -0.5 -0.5 ]]
array2 =
[[ 1. 0. 0. -1.25 0.33333333 -0.16666667]
[ 0. 1. -1.16666667 0. -0.66666667 0.83333333]
[ 0. 0. 1. 1.5 -0.5 -0.5 ]]
array2 =
[[ 1. 0. 0. -1.25 0.75 -0.16666667]
[ 0. 1. -1.16666667 0. -0.66666667 0.83333333]
[ 0. 0. 1. 1.5 -0.5 -0.5 ]]
array2 =
[[ 1. 0. 0. -1.25 0.75 0.25 ]
[ 0. 1. -1.16666667 0. -0.66666667 0.83333333]
[ 0. 0. 1. 1.5 -0.5 -0.5 ]]
array2 =
[[ 1. 0. 0. -1.25 0.75 0.25 ]
[ 0. 1. 0. 0. -0.66666667 0.83333333]
[ 0. 0. 1. 1.5 -0.5 -0.5 ]]
array2 =
[[ 1. 0. 0. -1.25 0.75 0.25 ]
[ 0. 1. 0. 1.75 -0.66666667 0.83333333]
[ 0. 0. 1. 1.5 -0.5 -0.5 ]]
array2 =
[[ 1. 0. 0. -1.25 0.75 0.25 ]
[ 0. 1. 0. 1.75 -1.25 0.83333333]
[ 0. 0. 1. 1.5 -0.5 -0.5 ]]
array2 =
[[ 1. 0. 0. -1.25 0.75 0.25]
[ 0. 1. 0. 1.75 -1.25 0.25]
[ 0. 0. 1. 1.5 -0.5 -0.5 ]]
=====================
determinant of matrix
=====================
det = 4.0
Pythonコードを以下に記す.
#!/usr/bin/python2.7
# -*- coding: utf-8 -*-
import numpy as np
print('=============')
print('3 * 6 array')
print('=============')
array2 = np.array([[3.,1.,2.,1.,0.,0.],[5.,1.,3.,0.,1.,0.],[4.,2.,1.,0.,0.,1.]])
arrayp = np.array([[0.,0.,0.],[0.,0.,0.],[0.,0.,0.]])
arrayq = np.array([[0.,0.,0.],[0.,0.,0.],[0.,0.,0.]])
tmp = np.array([[0.,0.,0.,0.,0.,0.],[0.,0.,0.,0.,0.,0.],[0.,0.,0.,0.,0.,0.]])
print('array2 =')
print(array2)
print('=====================================')
print('Gauss-Jordan Elimination(掃き出し法)')
print('=====================================')
#n = 3 ##元の数はこの場合,3
det = 1.0
for k in range(0,3,1):
for i in range(k,3,1):
if i !=k:
if abs(array2[k,k]) < abs(array2[i,k]):
# j列の範囲をk列目から5列目に変更する.
for j in range(k,6,1):
tmp[i,j] = array2[k,j]
array2[k,j] = array2[i,j]
array2[i,j] = tmp[i,j]
arrayp[k,k] = array2[k,k]
det = det * arrayp[k,k]
# j列の範囲をk列目から5列目に変更する.
for j in range(k,5,1):
array2[k,j+1] = array2[k,j+1] / arrayp[k,k]
print('array2 =')
print(array2)
array2[k,k] = 1
for i in range(0,3,1):
arrayq[k,k] = array2[i,k]
# j列の範囲をk列目から5列目に変更する.
for j in range(k,6,1):
if i != k:
array2[i,j] = array2[i,j] - array2[k,j] * arrayq[k,k]
print('array2 =')
print(array2)
else:
array2[i,j] = array2[i,j]
##行列式の計算結果
print('=====================')
print('determinant of matrix')
print('=====================')
print('det = {}'.format(det))
行列とその逆行列の積が単位行列になるか確かめてみる.
==========================
3 * 3 array multiplication
==========================
a =
[[ 3. 1. 2.]
[ 5. 1. 3.]
[ 4. 2. 1.]]
b =
[[-1.25 0.75 0.25]
[ 1.75 -1.25 0.25]
[ 1.5 -0.5 -0.5 ]]
c =
[[ 1. 0. 0.]
[ 0. 1. 0.]
[ 0. 0. 1.]]
単位行列が得られることが確かめられた.
上記,$3\times3$行列同士を乗算するコードを以下に記す.
#!/usr/bin/python2.7
# -*- coding: utf-8 -*-
import numpy as np
print('==========================')
print('3 * 3 array multiplication')
print('==========================')
a = np.array([[3.,1.,2.],[5.,1.,3.],[4.,2.,1.]])
b = np.array([[-1.25,0.75,0.25],[1.75,-1.25,0.25],[1.5,-0.5,-0.5]])
c = np.array([[0.,0.,0.],[0.,0.,0.],[0.,0.,0.]])
print('a =')
print(a)
print('b =')
print(b)
for i in range(0,3,1):
for j in range(0,3,1):
for k in range(0,3,1):
c[i,j] += a[i,k]*b[k,j]
print('c =')
print(c)
index.htmlに戻る