概要
Gaussの消去法は連立一次方程式の解法の1つです。連立方程式を行列式として解く手法の中でも最も基本的なものです。
大学で習う線形代数では手計算でやっていましたが、なかなか面倒くさかったのを覚えています。
実務では表計算ソフトやシミュレータで連立方程式を解くことがほとんどだと思いますので、ある意味楽になりましたね。
ここでは化学反応式の係数決定問題を例に、PythonでGaussの消去法を計算してみます。
Gaussの消去法
原理
ここでは次の3元連立一次方程式を解くことを考えます。
$$a_{11}x_{1}+a_{12}x_{2}+a_{13}x_{3}=b_{1}・・・(1)$$
$$a_{21}x_{1}+a_{22}x_{2}+a_{23}x_{3}=b_{2}・・・(2)$$
$$a_{31}x_{1}+a_{32}x_{2}+a_{33}x_{3}=b_{3}・・・(3)$$
Gaussの消去法(掃き出し法)は、(1)～(3)式に方程式であれば成り立つ変形を行なうことで未知数を消去する手法です。
まず(1)式を-a21/a11倍すると
$$-a_{21}x_{1}-\frac{a_{21}a_{12}}{a_{11}}x_{2}-\frac{a_{21}a_{13}}{a_{11}}x_{3}=-\frac{a_{21}b_{1}}{a_{11}}・・・(4)$$
(4)式となります。
次に(4)式と(2)式を辺々足し合わせることでx1の項を消去することができます。
$$(a_{22}-\frac{a_{21}a_{12}}{a_{11}})x_{2}+(a_{23}-\frac{a_{21}a_{13}}{a_{11}})x_{3}=(b_{2}-\frac{a_{21}b_{1}}{a_{11}})・・・(5)$$
同じような操作を(3)式に対しても行ないます。
(3)式の第一項を消せるように、(1)式を-a31/a11倍します。
$$-a_{31}x_{1}-\frac{a_{31}a_{12}}{a_{11}}x_{2}-\frac{a_{31}a_{13}}{a_{11}}x_{3}=-\frac{a_{31}b_{1}}{a_{11}}・・・(6)$$
次に(6)式と(3)式を辺々足し合わせることでx1の項を消去します。
$$(a_{32}-\frac{a_{31}a_{12}}{a_{11}})x_{2}+(a_{33}-\frac{a_{31}a_{13}}{a_{11}})x_{3}=(b_{3}-\frac{a_{31}b_{1}}{a_{11}})・・・(7)$$
(5),(7)式について係数がややこしくなったので、簡単のために
$${a_{22}}'x_{2}+{a_{23}}'x_{3}={b_{2}}'・・・(8)$$
$${a_{32}}'x_{2}+{a_{33}}'x_{3}={b_{3}}'・・・(9)$$
と係数を書き換えます。
続いて(8)式を-a32'/a22'倍します。
$$-{a_{32}}'x_{2}-\frac{{a_{32}}'{a_{23}}'}{{a_{22}}'}x_{3}=-\frac{{a_{32}}'{b_{2}}'}{{a_{22}}'}・・・(10)$$
次に(10)式と(9)式を辺々足し合わせます。
$$({a_{33}}'-\frac{{a_{32}}'{a_{23}}'}{{a_{22}}'})x_{3}=({b_{3}}'-\frac{{a_{32}}'{b_{2}}'}{{a_{22}}'})・・・(11)$$
(11)式のようにx2の項が消去できました。
(11)式の係数を簡単のために、
$${a_{33}}''x_{3}={b_{3}}''・・・(12)$$
(12)式に書き換えます。
以上、(1),(8),(12)式をまとめて書くと、
$$a_{11}x_{1}+a_{12}x_{2}+a_{13}x_{3}=b_{1}・・・(1)$$
$${a_{22}}'x_{2}+{a_{23}}'x_{3}={b_{2}}'・・・(8)$$
$${a_{33}}''x_{3}={b_{3}}''・・・(12)$$
となります。
行列式で書くと
\(\begin{pmatrix}a_{11}&a_{12}&a_{13}\\0&{a_{22}}'&{a_{23}}'\\0&0&{a_{33}}''\end{pmatrix}
\begin{pmatrix}x_{1}\\x_{2}\\x_{3}\end{pmatrix}
=\begin{pmatrix}b_{1}\\{b_{2}}'\\{b_{3}}''\end{pmatrix}\)・・・(13)
(13)式となります。
ちなみに、(13)式のように左下部分の係数を0にする操作を前進消去といいます。このような操作を行なうメリットはx1～x3を簡単に求めることができることです。
x1,x2の項を消去したことで(12)式からx3を簡単に求めることができます。
求めたx3を(8)式に代入することでx2も算出することができます。最後にx2,x3を(1)式に代入することでx1も求めることができます。
このようにx3→x2→x1と下側から解を求めて代入する操作を後退代入といいます。(別に用語は覚える必要はありません。)
以上が大まかなガウスの消去法の原理です。
部分ピポット選択
上で説明した原理には1つ欠点があります。対角成分であるa11,a22,a33で割り算をするので、いずれかの値がゼロであればゼロ割になって計算できないことです。
そこで行列式の行を入れ替えてゼロ割を回避する手法を部分ピポット選択といいます。
\(\begin{pmatrix}0&a_{12}&a_{13}\\a_{21}&a_{22}&a_{23}\\a_{31}&a_{32}&a_{33}\end{pmatrix}
\begin{pmatrix}x_{1}\\x_{2}\\x_{3}\end{pmatrix}
=\begin{pmatrix}b_{1}\\b_{2}\\b_{3}\end{pmatrix}\)・・・(14)
仮に(14)式のような状態だとすれば、
\(\begin{pmatrix}a_{21}&a_{22}&a_{23}\\0&a_{12}&a_{13}\\a_{31}&a_{32}&a_{33}\end{pmatrix}
\begin{pmatrix}x_{1}\\x_{2}\\x_{3}\end{pmatrix}
=\begin{pmatrix}b_{2}\\b_{1}\\b_{3}\end{pmatrix}\)・・・(15)
(15)式のように入れ替えればゼロ割を回避できます。
行列式だと難しいことをしているように見えますが、実際は
$$a_{12}x_{2}+a_{13}x_{3}=b_{1}・・・(16)$$
$$a_{21}x_{1}+a_{22}x_{2}+a_{23}x_{3}=b_{2}・・・(17)$$
$$a_{31}x_{1}+a_{32}x_{2}+a_{33}x_{3}=b_{3}・・・(18)$$
これら(16),(17),(18)式を
$$a_{21}x_{1}+a_{22}x_{2}+a_{23}x_{3}=b_{2}・・・(17)$$
$$a_{12}x_{2}+a_{13}x_{3}=b_{1}・・・(16)$$
$$a_{31}x_{1}+a_{32}x_{2}+a_{33}x_{3}=b_{3}・・・(18)$$
(17),(16),(18)式に並び替えただけです。何の問題もない操作ですよね。
また、入れ替える行は0でなければどこでも良いわけではなく、なるべく絶対値が大きい値の行と入れ替えた方が計算誤差が少なくなるようです。
Pythonでの計算法
Pythonの記述はこちらのサイト様の記事を参考にさせていただきました。
計算条件
逆行列の記事で扱った連立方程式をGaussの消去法で解いてみましょう。
$$C_{10}H_{20}+aO_{2}+bNH_{3} → cCH_{1.84}O_{0.41}N_{0.18}+dCO_{2}+eH_{2}O・・・(19)$$
$$\frac{生成したCO_{2}のモル数}{消費されたO_{2}のモル数}=0.48・・・(20)$$
(19),(20)式から、(19)の反応式の係数a,b,c,d,eを決定します。
細かい導出式は逆行列の記事で解説しているので省きますが、最終的に行列式は
\(\begin{pmatrix}0&0&1&1&0\\0&-3&1.84&0&2\\-2&0&0.41&2&1\\0&-1&0.18&0&0\\0.48&0&0&-1&0\end{pmatrix}
\begin{pmatrix}a\\b\\c\\d\\e\end{pmatrix}
=\begin{pmatrix}10\\20\\0\\0\\0\end{pmatrix}\)・・・(21)
(21)式となります。この式をGaussの消去法で解いてみましょう。
ライブラリ・モジュールのインポート
ここからPythonでの記述となります。Numpyをインストールします。
#モジュールインポート
import numpy as np前提条件
行列AとベクトルBをそれぞれ入力します。
#前提条件の入力
A=np.array([[0, 0, 1, 1, 0], [0, -3, 1.84, 0, 2], [-2, 0, 0.41, 2, 1],[0, -1, 0.18, 0, 0],[0.48, 0, 0, -1, 0]])
B=np.array([10, 20, 0, 0, 0])部分ピポット選択の入力
部分ピポット選択の計算式をdef文で入力します。
#部分ピポット選択
def select_pivot(xs, i):
k = np.abs(xs[i:,i]).argmax() + i
if k != i:
temp = xs[i].copy()
xs[i] = xs[k]
xs[k] = temp2行目で対称列で最も絶対値が高い値を検索しています。
3行目以降で、最も絶対値が高い値が対角成分以外の位置にあれば入れ替える操作を行なっています。
Gaussの消去法の入力
Gaussの消去法をdef文で入力します。
#Gaussの消去法
def gauss1(xs, ys):
# 拡大係数行列の生成
n = len(xs)
zs = np.c_[xs.astype(np.float_), ys.astype(np.float_)]
# 前進消去
for i in range(n - 1):
select_pivot(zs, i)
for j in range(i + 1, n):
temp = zs[j, i] / zs[i, i]
zs[j, i:] -= temp * zs[i, i:]
print(zs)
# 後退代入
for i in range(n - 1, -1, -1):
zs[i, n] -= zs[i, i+1:n] @ zs[i+1:, n]
zs[i, n] /= zs[i, i]
return zs[:, n]3～5行目では拡大係数行列を作成しています。
拡大係数行列とはAx=bの形から、Aとbをくっつけた行列のことをいいます。今回では、
\(\begin{pmatrix}0&0&1&1&0\\0&-3&1.84&0&2\\-2&0&0.41&2&1\\0&-1&0.18&0&0\\0.48&0&0&-1&0\end{pmatrix}\)
上の行列Aに、
\(\begin{pmatrix}0&0&1&1&0&10\\0&-3&1.84&0&2&20\\-2&0&0.41&2&1&0\\0&-1&0.18&0&0&0\\0.48&0&0&-1&0&0\end{pmatrix}\)
ベクトルbをくっつけて表した行列です。
その後のGaussの消去法の操作自体は変わりません。
6～12行目で前進消去の操作を入力しています。
13～17行目で後退代入の操作を入力しています。
Gaussの消去法の実行
定義したGaussの消去法の実行を記述しています。
#結果の出力
print(gauss1(A,B))計算結果の確認
入力に間違いがなければ下のように結果が表示されます。
[[-2. 0. 0.41 2. 1. 0. ]
[ 0. -3. 1.84 0. 2. 20. ]
[ 0. 0. 1. 1. 0. 10. ]
[ 0. 0. 0. -0.6184 0.24 -0.984 ]
[ 0. 0. 0. 0. -0.49849073 -3.02285468]]
[8.21799308 1.0899654 6.05536332 3.94463668 6.06401384]1～6行目に表示されている行列が前進消去が最後まで終わった後の値となります。
きれいに下三角部分が0になっていますね。
7行目の値が解になります。
a=8.22
b=1.09
c=6.06
d=3.94
e=6.06
となりました。
全体のコードを下に載せておきます。
#モジュールインポート
import numpy as np
#前提条件の入力
A=np.array([[0, 0, 1, 1, 0], [0, -3, 1.84, 0, 2], [-2, 0, 0.41, 2, 1],[0, -1, 0.18, 0, 0],[0.48, 0, 0, -1, 0]])
B=np.array([10, 20, 0, 0, 0])
#部分ピポット選択
def select_pivot(xs, i):
k = np.abs(xs[i:,i]).argmax() + i
if k != i:
temp = xs[i].copy()
xs[i] = xs[k]
xs[k] = temp
# Gaussの消去法
def gauss1(xs, ys):
# 拡大係数行列の生成
n = len(xs)
zs = np.c_[xs.astype(np.float_), ys.astype(np.float_)]
# 前進消去
for i in range(n - 1):
select_pivot(zs, i)
for j in range(i + 1, n):
temp = zs[j, i] / zs[i, i]
zs[j, i:] -= temp * zs[i, i:]
print(zs)
# 後退代入
for i in range(n - 1, -1, -1):
zs[i, n] -= zs[i, i+1:n] @ zs[i+1:, n]
zs[i, n] /= zs[i, i]
return zs[:, n]
#結果の出力
print(gauss1(A,B))まとめ
Gaussの消去法の解説をしました。逆行列で求める手法よりは計算負荷が低く、手計算でもできる点が魅力です。