物理のノート
行列式の数値計算
上三角型行列の行列式
$3\times 3$行列において上三角型の行列式（determinant）は, 対角要素の積で表された.すなわち, \begin{equation} \lvert A \rvert= \begin{vmatrix} a_{1,1} & a_{1,2} & a_{1,3} \\ 0 & a_{2,2} & a_{2,3}\\ 0 & 0 & a_{3,3}\\ \end{vmatrix} \end{equation} の行列式は, \begin{equation} \lvert A \rvert=a_{1,1}a_{2,2} a_{3,3}+0 a_{2,3} a_{1,2}+0 0 a_{1,3} -(a_{1,1}a_{2,3} 0 +a_{1,2}0 a_{3,3} +a_{1,3}a_{2,2} 0)\\ =a_{1,1}a_{2,2} a_{3,3} \end{equation} となる. 導出方法は異なるが,これと同じ計算結果は,ガウスの単純消去法と呼ばれる,連立1次方程式の数値解法の前進消去（Gaussian Elimination）の操作で得ることができる.例えば,次のような$3\times 3$行列を考える. \begin{equation} \begin{pmatrix} 3 & 1 & 2 \\ 5 & 1 & 3 \\ 4 & 2 & 1 \\ \end{pmatrix} \end{equation} この行列を上三角型に変形して,対角要素の積を求める.
ここではそのアルゴリズムをpythonで実装した.
# -*- coding: utf-8 -*-
import numpy as np
print('===========')
print('3 * 3 array')
print('===========')
array2 = np.array([[3.,1.,2.],[5.,1.,3.],[4.,2.,1.]])
print('行列array2')
print(array2)
##前進消去
print('=================')
print('Gauss Elimination')
print('=================')
for p in range(0,2,1):
print('pivot:array2[p,p] = {}'.format(array2[p,p]))
for i in range(p+1,3,1):
m = (array2[i,p]) / (array2[p,p])
print('m = {}'.format(m))
for j in range(p,3,1):
array2[i,j] = array2[i,j] - m * (array2[p,j])
print('行列array2')
print(array2)
##行列式の計算
det = 1.0
for a in range(0,3,1):
det = det * array2[a,a]
print('=====================')
print('Determinant of matrix')
print('=====================')
print('det = {}'.format(det))
Canopyでの実行結果は次のようになった. ターミナルで実行しても同様である.
===========
3 * 3 array
===========
行列array2
[[ 3. 1. 2.]
[ 5. 1. 3.]
[ 4. 2. 1.]]
=================
Gauss Elimination
=================
pivot:array2[p,p] = 3.0
m = 1.66666666667
m = 1.33333333333
pivot:array2[p,p] = -0.666666666667
m = -1.0
行列array2
[[ 3. 1. 2. ]
[ 0. -0.66666667 -0.33333333]
[ 0. 0. -2. ]]
=====================
Determinant of matrix
=====================
det = 4.0
また$m$は, \begin{equation} m_{i,p} = \frac{a_{i,p}}{a_{p,p}} \end{equation} $i = 2,3$（$3\times 3$行列の場合）
で表され, 上三角型行列に変形するために,$pivot$の値を用いてその都度計算で得られる変数である. 上記行列の行列式（$det$）は, \[ 3 \times \left( -\frac{2}{3} \right) \times (-2) = 4 \] である.
掃き出し法による行列式の計算
ガウスの単純消去法では行列を上三角型に変形する計算を行なった.すなわち,$p \times p$行列の場合, 2行目から$p$行目の対角要素の左下にある要素が$0$となるように計算した.これに対して,対角要素に$1$のみが残るように計算する連立1次方程式の数値解法にGauss-Jordanの消去法がある.これを応用して行列式を求める方法を考える. まず,Gauss-Jordanの消去法のアルゴリズムを確かめるため,次のような3元連立1次方程式を例にGauss-Jordanの消去法で解と行列式を求める.
方程式の左辺の係数にだけ着目すると,先にGaussの単純消去法で行列式を求めた行列の要素と同じである. \begin{equation} 3 x_1 + x_2 + 2 x_3 = 13 \\ 5 x_1 + x_2 + 3 x_3 = 20 \\ 4 x_1 + 2 x_2 + x_3 = 13 \end{equation} Gauss-Jordanの消去法のアルゴリズムをPythonで実装する.
#!/usr/bin/python2.7
# -*- coding: utf-8 -*-
import numpy as np
print('=============')
print('3 * 4 array')
print('=============')
array2 = np.array([[3.,1.,2.,13.],[5.,1.,3.,20.],[4.,2.,1.,13.]])
arrayp = np.array([[0.,0.,0.],[0.,0.,0.],[0.,0.,0.]])
arrayq = np.array([[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 #行列式の初期値を1に設定する.
for k in range(0,3,1):
# pivot（枢軸要素）の設定
arrayp[k,k] = array2[k,k]
det = det * arrayp[k,k]
# k行目の各列j+1の要素を,k行j列の要素で割る.
for j in range(k,3,1):
array2[k,j+1] = array2[k,j+1] / arrayp[k,k]
print('array2 =') #計算途中の結果を出力
print(array2) #計算途中の結果を出力
# k行k列を1に設定する.
array2[k,k] = 1
for i in range(0,3,1):
# i行目のk列（左辺の0,すなわち掃き出された要素の右隣の要素）を配列arrayqに収める.
arrayq[k,k] = array2[i,k]
# i行目以外の行の各列の要素を,arrayqに納めた各列の要素とk行j列の積で引く.
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('x_1 = {0}, x_2 = {1}, x_3 = {2}'.format(array2[0,3], array2[1,3], array2[2,3]))
##行列式の計算結果
print('=====================')
print('determinant of matrix')
print('=====================')
print('det = {}'.format(det))
実行結果は次の通り.
=============
3 * 4 array
=============
array2 =
[[ 3. 1. 2. 13.]
[ 5. 1. 3. 20.]
[ 4. 2. 1. 13.]]
========================
Gauss-Jordan Elimination
========================
array2 =
[[ 3. 0.33333333 2. 13. ]
[ 5. 1. 3. 20. ]
[ 4. 2. 1. 13. ]]
array2 =
[[ 3. 0.33333333 0.66666667 13. ]
[ 5. 1. 3. 20. ]
[ 4. 2. 1. 13. ]]
array2 =
[[ 3. 0.33333333 0.66666667 4.33333333]
[ 5. 1. 3. 20. ]
[ 4. 2. 1. 13. ]]
array2 =
[[ 1. 0.33333333 0.66666667 4.33333333]
[ 0. 1. 3. 20. ]
[ 4. 2. 1. 13. ]]
array2 =
[[ 1. 0.33333333 0.66666667 4.33333333]
[ 0. -0.66666667 3. 20. ]
[ 4. 2. 1. 13. ]]
array2 =
[[ 1. 0.33333333 0.66666667 4.33333333]
[ 0. -0.66666667 -0.33333333 20. ]
[ 4. 2. 1. 13. ]]
array2 =
[[ 1. 0.33333333 0.66666667 4.33333333]
[ 0. -0.66666667 -0.33333333 -1.66666667]
[ 4. 2. 1. 13. ]]
array2 =
[[ 1. 0.33333333 0.66666667 4.33333333]
[ 0. -0.66666667 -0.33333333 -1.66666667]
[ 0. 2. 1. 13. ]]
array2 =
[[ 1. 0.33333333 0.66666667 4.33333333]
[ 0. -0.66666667 -0.33333333 -1.66666667]
[ 0. 0.66666667 1. 13. ]]
array2 =
[[ 1. 0.33333333 0.66666667 4.33333333]
[ 0. -0.66666667 -0.33333333 -1.66666667]
[ 0. 0.66666667 -1.66666667 13. ]]
array2 =
[[ 1. 0.33333333 0.66666667 4.33333333]
[ 0. -0.66666667 -0.33333333 -1.66666667]
[ 0. 0.66666667 -1.66666667 -4.33333333]]
array2 =
[[ 1. 0.33333333 0.66666667 4.33333333]
[ 0. -0.66666667 0.5 -1.66666667]
[ 0. 0.66666667 -1.66666667 -4.33333333]]
array2 =
[[ 1. 0.33333333 0.66666667 4.33333333]
[ 0. -0.66666667 0.5 2.5 ]
[ 0. 0.66666667 -1.66666667 -4.33333333]]
array2 =
[[ 1. 0. 0.66666667 4.33333333]
[ 0. 1. 0.5 2.5 ]
[ 0. 0.66666667 -1.66666667 -4.33333333]]
array2 =
[[ 1. 0. 0.5 4.33333333]
[ 0. 1. 0.5 2.5 ]
[ 0. 0.66666667 -1.66666667 -4.33333333]]
array2 =
[[ 1. 0. 0.5 3.5 ]
[ 0. 1. 0.5 2.5 ]
[ 0. 0.66666667 -1.66666667 -4.33333333]]
array2 =
[[ 1. 0. 0.5 3.5 ]
[ 0. 1. 0.5 2.5 ]
[ 0. 0. -1.66666667 -4.33333333]]
array2 =
[[ 1. 0. 0.5 3.5 ]
[ 0. 1. 0.5 2.5 ]
[ 0. 0. -2. -4.33333333]]
array2 =
[[ 1. 0. 0.5 3.5]
[ 0. 1. 0.5 2.5]
[ 0. 0. -2. -6. ]]
array2 =
[[ 1. 0. 0.5 3.5]
[ 0. 1. 0.5 2.5]
[ 0. 0. -2. 3. ]]
array2 =
[[ 1. 0. 0. 3.5]
[ 0. 1. 0.5 2.5]
[ 0. 0. 1. 3. ]]
array2 =
[[ 1. 0. 0. 2. ]
[ 0. 1. 0.5 2.5]
[ 0. 0. 1. 3. ]]
array2 =
[[ 1. 0. 0. 2. ]
[ 0. 1. 0. 2.5]
[ 0. 0. 1. 3. ]]
array2 =
[[ 1. 0. 0. 2.]
[ 0. 1. 0. 1.]
[ 0. 0. 1. 3.]]
x_1 = 2.0, x_2 = 1.0, x_3 = 3.0
=====================
determinant of matrix
=====================
det = 4.0
行列式は$4$で, Gaussの単純消去法の前進消去で求めた結果と一致する.
なお, Gaussの単純消去法で連立1次方程式を解く場合,後退代入（Back Substituion）によって変数の値を求める.しかし,Gauss-Jordan消去法で求める場合は,対角型の要素（$1$である）が得られた時点で解が求められるため,後退代入は不要である.ただし計算量は増える（後述） さて, 行列式を求める場合, 連立方程式の右辺にある列は不要なので左辺の係数のみの要素からなる$3 \times 3$行列の行列式を求める掃き出しアルゴリズムに改良する.下記コードには2カ所の改良点のみコメントで示している.
#!/usr/bin/python2.7
# -*- coding: utf-8 -*-
import numpy as np
print('=============')
print('3 * 3 array')
print('=============')
array2 = np.array([[3.,1.,2.],[5.,1.,3.],[4.,2.,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.]])
print('array2 =')
print(array2)
print('=====================================')
print('Gauss-Jordan Elimination（掃き出し法）')
print('=====================================')
#n = 3 ##元の数はこの場合,3
det = 1.0
for k in range(0,3,1):
arrayp[k,k] = array2[k,k]
det = det * arrayp[k,k]
# j列の範囲を0列目から2列目に変更する.
for j in range(k,2,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列の範囲を0列目から2列目に変更する.
for j in range(k,3,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
=============
array2 =
[[ 3. 1. 2.]
[ 5. 1. 3.]
[ 4. 2. 1.]]
=====================================
Gauss-Jordan Elimination（掃き出し法）
=====================================
array2 =
[[ 3. 0.33333333 2. ]
[ 5. 1. 3. ]
[ 4. 2. 1. ]]
array2 =
[[ 3. 0.33333333 0.66666667]
[ 5. 1. 3. ]
[ 4. 2. 1. ]]
array2 =
[[ 1. 0.33333333 0.66666667]
[ 0. 1. 3. ]
[ 4. 2. 1. ]]
array2 =
[[ 1. 0.33333333 0.66666667]
[ 0. -0.66666667 3. ]
[ 4. 2. 1. ]]
array2 =
[[ 1. 0.33333333 0.66666667]
[ 0. -0.66666667 -0.33333333]
[ 4. 2. 1. ]]
array2 =
[[ 1. 0.33333333 0.66666667]
[ 0. -0.66666667 -0.33333333]
[ 0. 2. 1. ]]
array2 =
[[ 1. 0.33333333 0.66666667]
[ 0. -0.66666667 -0.33333333]
[ 0. 0.66666667 1. ]]
array2 =
[[ 1. 0.33333333 0.66666667]
[ 0. -0.66666667 -0.33333333]
[ 0. 0.66666667 -1.66666667]]
array2 =
[[ 1. 0.33333333 0.66666667]
[ 0. -0.66666667 0.5 ]
[ 0. 0.66666667 -1.66666667]]
array2 =
[[ 1. 0. 0.66666667]
[ 0. 1. 0.5 ]
[ 0. 0.66666667 -1.66666667]]
array2 =
[[ 1. 0. 0.5 ]
[ 0. 1. 0.5 ]
[ 0. 0.66666667 -1.66666667]]
array2 =
[[ 1. 0. 0.5 ]
[ 0. 1. 0.5 ]
[ 0. 0. -1.66666667]]
array2 =
[[ 1. 0. 0.5]
[ 0. 1. 0.5]
[ 0. 0. -2. ]]
array2 =
[[ 1. 0. 0. ]
[ 0. 1. 0.5]
[ 0. 0. 1. ]]
array2 =
[[ 1. 0. 0.]
[ 0. 1. 0.]
[ 0. 0. 1.]]
=====================
determinant of matrix
=====================
det = 4.0
行列式の計算は次のようになる. \[ 3 \times \left(-\frac{2}{3}\right) \times (-2) = 4 \]
部分ピボット選択（行交換）を伴う掃き出し法
ここまでのGaussの単純消去法,およびGauss-Jordan消去法では,正方行列の対角線上に並ぶ要素が$0$の場合,除算ができない.例えば,次の連立1次方程式のような場合である. \begin{equation} \phantom{XXXXX} 36 x_2 + 71 x_3 = 100 \\ -36 x_1 \phantom{XXXX} + 68 x_3 = 50 \\ -75 x_1 -70 x_2 \phantom{XXXX} = 0 \end{equation} この連立方程式を先のGauss-Jordan消去法を用いて計算すると次のような結果になる.
=============
3 * 4 array
=============
array2 =
[[ 0. 36. 71. 100.]
[ -36. 0. 68. 50.]
[ -75. -70. 0. 0.]]
========================
Gauss-Jordan Elimination
========================
array2 =
[[ 0. inf 71. 100.]
[ -36. 0. 68. 50.]
[ -75. -70. 0. 0.]]
array2 =
[[ 0. inf inf 100.]
[ -36. 0. 68. 50.]
[ -75. -70. 0. 0.]]
array2 =
[[ 0. inf inf inf]
[-36. 0. 68. 50.]
[-75. -70. 0. 0.]]
array2 =
[[ 1. inf inf inf]
[ 0. 0. 68. 50.]
[-75. -70. 0. 0.]]
array2 =
[[ 1. inf inf inf]
[ 0. inf 68. 50.]
[-75. -70. 0. 0.]]
array2 =
[[ 1. inf inf inf]
[ 0. inf inf 50.]
[-75. -70. 0. 0.]]
array2 =
[[ 1. inf inf inf]
[ 0. inf inf inf]
[-75. -70. 0. 0.]]
array2 =
[[ 1. inf inf inf]
[ 0. inf inf inf]
[ 0. -70. 0. 0.]]
array2 =
[[ 1. inf inf inf]
[ 0. inf inf inf]
[ 0. inf 0. 0.]]
array2 =
[[ 1. inf inf inf]
[ 0. inf inf inf]
[ 0. inf inf 0.]]
array2 =
[[ 1. inf inf inf]
[ 0. inf inf inf]
[ 0. inf inf inf]]
array2 =
[[ 1. inf inf inf]
[ 0. inf nan inf]
[ 0. inf inf inf]]
array2 =
[[ 1. inf inf inf]
[ 0. inf nan nan]
[ 0. inf inf inf]]
array2 =
[[ 1. nan inf inf]
[ 0. 1. nan nan]
[ 0. inf inf inf]]
array2 =
[[ 1. nan nan inf]
[ 0. 1. nan nan]
[ 0. inf inf inf]]
array2 =
[[ 1. nan nan nan]
[ 0. 1. nan nan]
[ 0. inf inf inf]]
array2 =
[[ 1. nan nan nan]
[ 0. 1. nan nan]
[ 0. nan inf inf]]
array2 =
[[ 1. nan nan nan]
[ 0. 1. nan nan]
[ 0. nan nan inf]]
array2 =
[[ 1. nan nan nan]
[ 0. 1. nan nan]
[ 0. nan nan nan]]
array2 =
[[ 1. nan nan nan]
[ 0. 1. nan nan]
[ 0. nan nan nan]]
array2 =
[[ 1. nan nan nan]
[ 0. 1. nan nan]
[ 0. nan 1. nan]]
array2 =
[[ 1. nan nan nan]
[ 0. 1. nan nan]
[ 0. nan 1. nan]]
array2 =
[[ 1. nan nan nan]
[ 0. 1. nan nan]
[ 0. nan 1. nan]]
array2 =
[[ 1. nan nan nan]
[ 0. 1. nan nan]
[ 0. nan 1. nan]]
x_1 = nan, x_2 = nan, x_3 = nan
=====================
determinant of matrix
=====================
det = nan
$1$行$1$列目のピボットが$0$なので,$1$行目$2$列目から$4$列目までの要素を割ると無限大（infinity）になっている.
無限大を用いて計算するとNaN（Not a Number：非数）という欠損値が表示される.$2$行$3$列にNaNと記されたところは \[ \infty - \infty\times \infty \] という計算がなされている. そこで,要素が$0$ではない行と入れ替えるアルゴリズムの修正を行う.この修正を行なったアルゴリズムは部分ピボット選択による,または行交換を伴う（Gaussの）消去法（掃き出し法）と呼ばれる.
以下に実装例を示した.
ピボットのある列の中で,ピボットのある行以下の要素のうち,絶対値が最大のものを新たにピボットに採用する.次に,旧ピボットを含む行のすべての要素を,新ピボットを含む行の要素と入れ替える処理を行っている.
#!/usr/bin/python2.7
# -*- coding: utf-8 -*-
import numpy as np
print('=============')
print('3 * 4 array')
print('=============')
array2 = np.array([[0.,36.,71.,100.],[-36.,0.,68.0,50.],[-75.,-70.,0.,0.]])
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.]])
print('array2 =')
print(array2)
print('=====================================')
print('Gauss-Jordan Elimination（掃き出し法)')
print('=====================================')
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]):
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]
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,3,1):
arrayq[k,k] = array2[i,k]
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('x_1 = {0}, x_2 = {1}, x_3 = {2}'.format(array2[0,3], array2[1,3], array2[2,3]))
##行列式の計算結果
print('=====================')
print('determinant of matrix')
print('=====================')
print('det = {}'.format(det))
先ほど解が得られなかった連立1次方程式の解と左辺の係数行列の行列式は次のようになった.
=============
3 * 4 array
=============
array2 =
[[ 0. 36. 71. 100.]
[ -36. 0. 68. 50.]
[ -75. -70. 0. 0.]]
=====================================
Gauss-Jordan Elimination（掃き出し法)
=====================================
array2 =
[[ -75. 0.93333333 0. 0. ]
[ 0. 36. 71. 100. ]
[ -36. 0. 68. 50. ]]
array2 =
[[ -75. 0.93333333 -0. 0. ]
[ 0. 36. 71. 100. ]
[ -36. 0. 68. 50. ]]
array2 =
[[ -75. 0.93333333 -0. -0. ]
[ 0. 36. 71. 100. ]
[ -36. 0. 68. 50. ]]
array2 =
[[ 1. 0.93333333 -0. -0. ]
[ 0. 36. 71. 100. ]
[ -36. 0. 68. 50. ]]
array2 =
[[ 1. 0.93333333 -0. -0. ]
[ 0. 36. 71. 100. ]
[ -36. 0. 68. 50. ]]
array2 =
[[ 1. 0.93333333 -0. -0. ]
[ 0. 36. 71. 100. ]
[ -36. 0. 68. 50. ]]
array2 =
[[ 1. 0.93333333 -0. -0. ]
[ 0. 36. 71. 100. ]
[ -36. 0. 68. 50. ]]
array2 =
[[ 1. 0.93333333 -0. -0. ]
[ 0. 36. 71. 100. ]
[ 0. 0. 68. 50. ]]
array2 =
[[ 1. 0.93333333 -0. -0. ]
[ 0. 36. 71. 100. ]
[ 0. 33.6 68. 50. ]]
array2 =
[[ 1. 0.93333333 -0. -0. ]
[ 0. 36. 71. 100. ]
[ 0. 33.6 68. 50. ]]
array2 =
[[ 1. 0.93333333 -0. -0. ]
[ 0. 36. 71. 100. ]
[ 0. 33.6 68. 50. ]]
array2 =
[[ 1. 0.93333333 -0. -0. ]
[ 0. 36. 1.97222222 100. ]
[ 0. 33.6 68. 50. ]]
array2 =
[[ 1. 0.93333333 -0. -0. ]
[ 0. 36. 1.97222222 2.77777778]
[ 0. 33.6 68. 50. ]]
array2 =
[[ 1. 0. -0. -0. ]
[ 0. 1. 1.97222222 2.77777778]
[ 0. 33.6 68. 50. ]]
array2 =
[[ 1. 0. -1.84074074 -0. ]
[ 0. 1. 1.97222222 2.77777778]
[ 0. 33.6 68. 50. ]]
array2 =
[[ 1. 0. -1.84074074 -2.59259259]
[ 0. 1. 1.97222222 2.77777778]
[ 0. 33.6 68. 50. ]]
array2 =
[[ 1. 0. -1.84074074 -2.59259259]
[ 0. 1. 1.97222222 2.77777778]
[ 0. 0. 68. 50. ]]
array2 =
[[ 1. 0. -1.84074074 -2.59259259]
[ 0. 1. 1.97222222 2.77777778]
[ 0. 0. 1.73333333 50. ]]
array2 =
[[ 1. 0. -1.84074074 -2.59259259]
[ 0. 1. 1.97222222 2.77777778]
[ 0. 0. 1.73333333 -43.33333333]]
array2 =
[[ 1. 0. -1.84074074 -2.59259259]
[ 0. 1. 1.97222222 2.77777778]
[ 0. 0. 1.73333333 -25. ]]
array2 =
[[ 1. 0. 0. -2.59259259]
[ 0. 1. 1.97222222 2.77777778]
[ 0. 0. 1. -25. ]]
array2 =
[[ 1. 0. 0. -48.61111111]
[ 0. 1. 1.97222222 2.77777778]
[ 0. 0. 1. -25. ]]
array2 =
[[ 1. 0. 0. -48.61111111]
[ 0. 1. 0. 2.77777778]
[ 0. 0. 1. -25. ]]
array2 =
[[ 1. 0. 0. -48.61111111]
[ 0. 1. 0. 52.08333333]
[ 0. 0. 1. -25. ]]
x_1 = -48.6111111111, x_2 = 52.0833333333, x_3 = -25.0
=====================
determinant of matrix
=====================
det = -4680.0
(3)式で示した$3 \times 3$の正方行列には$0$の要素はない.しかし,行交換を伴うGauss-Jordan消去法によってて行列式を計算してみよう.
コードは次のように作成した.
#!/usr/bin/python2.7
# -*- coding: utf-8 -*-
import numpy as np
print('=============')
print('3 * 3 array')
print('=============')
array2 = np.array([[3.,1.,2.],[5.,1.,3.],[4.,2.,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.]])
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]):
for j in range(k,3,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]
for j in range(k,2,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]
for j in range(k,3,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
=============
array2 =
[[ 3. 1. 2.]
[ 5. 1. 3.]
[ 4. 2. 1.]]
=====================================
Gauss-Jordan Elimination（掃き出し法)
=====================================
array2 =
[[ 5. 0.2 3. ]
[ 3. 1. 2. ]
[ 4. 2. 1. ]]
array2 =
[[ 5. 0.2 0.6]
[ 3. 1. 2. ]
[ 4. 2. 1. ]]
array2 =
[[ 1. 0.2 0.6]
[ 0. 1. 2. ]
[ 4. 2. 1. ]]
array2 =
[[ 1. 0.2 0.6]
[ 0. 0.4 2. ]
[ 4. 2. 1. ]]
array2 =
[[ 1. 0.2 0.6]
[ 0. 0.4 0.2]
[ 4. 2. 1. ]]
array2 =
[[ 1. 0.2 0.6]
[ 0. 0.4 0.2]
[ 0. 2. 1. ]]
array2 =
[[ 1. 0.2 0.6]
[ 0. 0.4 0.2]
[ 0. 1.2 1. ]]
array2 =
[[ 1. 0.2 0.6]
[ 0. 0.4 0.2]
[ 0. 1.2 -1.4]]
array2 =
[[ 1. 0.2 0.6 ]
[ 0. 1.2 -1.16666667]
[ 0. 0.4 0.2 ]]
array2 =
[[ 1. 0. 0.6 ]
[ 0. 1. -1.16666667]
[ 0. 0.4 0.2 ]]
array2 =
[[ 1. 0. 0.83333333]
[ 0. 1. -1.16666667]
[ 0. 0.4 0.2 ]]
array2 =
[[ 1. 0. 0.83333333]
[ 0. 1. -1.16666667]
[ 0. 0. 0.2 ]]
array2 =
[[ 1. 0. 0.83333333]
[ 0. 1. -1.16666667]
[ 0. 0. 0.66666667]]
array2 =
[[ 1. 0. 0. ]
[ 0. 1. -1.16666667]
[ 0. 0. 1. ]]
array2 =
[[ 1. 0. 0.]
[ 0. 1. 0.]
[ 0. 0. 1.]]
=====================
determinant of matrix
=====================
det = 4.0
行列式は同じく$4$であるが,行交換によって行列式のもとになるピボットが変化し,計算式は次のように変わった. \begin{equation} 5 \times 1.2 \times \frac{2}{3} = 4 \end{equation} 行列式は,行交換（行基本変形）によって変化しない. ピボットの選択は,対角要素が$0$になって計算処理の中に無限大（$\infty$）が入り込むことを防ぐほかに,
- 入力データの誤差の伝播を抑える.
- 計算途中の丸め誤差の影響を小さくする.
- ピボットとしてできるだけ絶対値の大きいものを選ぶことで途中の除算におけるオーバーフローの発生の予防が期待される.