物理のノート
逆行列
逆行列の数値計算（$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)