Moneytizerによる広告を掲載しています。概要Gaussの消去法は連立一次方程式の解法の1つです。連立方程式を行列式として解く手法の中でも最も基本的なものです。大学で習う線形代数では手計算でやっていましたが、なかなか面倒くさかったのを覚えています。実務では表計算ソフトやシミュレータで連立方程式を解くことがほとんどだと思いますので、ある意味楽になりましたね。ここでは化学反応式の係数決定問題を例に、Pythonで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でなければどこでも良いわけではなく、なるべく絶対値が大きい値の行と入れ替えた方が計算誤差が少なくなるようです。記事で扱った連立方程式をGaussの消去法で解いてみましょう。本記事では連立一次方程式の計算として、逆行列を使用した解法を解説しています。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をインストールします。行列AとベクトルBをそれぞれ入力します。部分ピポット選択の計算式をdef文で入力します。temp2行目で対称列で最も絶対値が高い値を検索しています。3行目以降で、最も絶対値が高い値が対角成分以外の位置にあれば入れ替える操作を行なっています。Gaussの消去法の入力Gaussの消去法をdef文で入力します。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の消去法の実行を記述しています。入力に間違いがなければ下のように結果が表示されます。6.06401384]1~6行目に表示されている行列が前進消去が最後まで終わった後の値となります。きれいに下三角部分が0になっていますね。7行目の値が解になります。a=8.22b=1.09c=6.06d=3.94e=6.06となりました。全体のコードを下に載せておきます。n]#結果の出力print(gauss1(A,B))まとめGaussの消去法の解説をしました。逆行列で求める手法よりは計算負荷が低く、手計算でもできる点が魅力です。コメントをキャンセルメールアドレスが公開されることはありません。※サイト次回のコメントで使用するためブラウザーに自分の名前、メールアドレス、サイトを保存する。date*日本語が含まれない投稿は無視されますのでご注意ください。【Secant法(割線法)】を解説:非線型方程式の解法Secant法(割線法)は非線形方程式を解く手法の1つです。method)とも呼ばれています。Python【Python】線形最小二乗法の計算例:アレニウスプロット本記事ではPythonの線形最小二乗法の計算例としてアレニウスプロットををやってみました。Pythonで最小二乗法をやる場合は計算手法まできちんと記述しないといけませんから、Excelよりも難易度が上がります。本記事では非線形最小二乗法の計算手法の1つである最急降下法について解説します。他の手法と比較して簡単に計算できますが、繰り返し計算回数が多くなり収束しにくいのが欠点です。【Runge-Kutta法】をわかりやすく解説:Euler法より精度の良い数値解法Runge-Kutta法(ルンゲクッタ法)は常微分方程式を解く手法の1つです。本記事ではRunge-Kutta法の計算手法について解説します。加えて、回分式反応器の一次反応に関する例題を解いてみて、Euler法とRunge-Kutta法の計算精度を比較してみました。Plus等の化学工学計算シミュレータを使用するにあたって悩ましい問題の1つが収束計算です。この記事では化学工学計算シミュレータでよく使用される収束計算手法について解説しています。加えて、どのように設定すれば収束しやすいか、というコツも紹介しています。Pythonでの化学工学計算例を紹介しています。今回はEuler法を使用して、タンク内の原料濃度変化を計算します。シミュレータの収束計算手法をわかりやすく解説プロフィールこのブログを運営しているルートです。化学メーカーで働くケミカルエンジニアです。