Numpyだけで書いたガウスの消去法で連立1次方程式を解いてみた

  • このエントリーをはてなブックマークに追加

ガウスの消去法は掃き出し法とも呼ばれ、機械的に連立1次方程式を解く方法で有名です。Pythonには様々な便利ライブラリがありますが、ここではNumpyだけを使ってガウスの消去法を書き、アルゴリズムの理解をすることを目標とします。

こんにちは。wat(@watlablog)です。ここでは連立1次方程式をプログラムで解くために、ガウスの消去法をゼロからPythonでコーディングします

ガウスの消去法の概要

ガウスの消去法とは?

ガウスの消去法(Gaussian Elimination)とは、掃き出し法(Row Reduction)とも呼ばれ、連立1次方程式を解くアルゴリズムの1つです。

ガウスの消去法は以下のような連立1次方程式に対して使います。

\[\left\{\begin{matrix} a_{11}x_{1}+a_{12}x_{2}+a_{13}x_{3}+a_{14}x_{4}=b_{1}\\ a_{21}x_{1}+a_{22}x_{2}+a_{23}x_{3}+a_{24}x_{4}=b_{2}\\ a_{31}x_{1}+a_{32}x_{2}+a_{33}x_{3}+a_{34}x_{4}=b_{3}\\ a_{41}x_{1}+a_{42}x_{2}+a_{43}x_{3}+a_{44}x_{4}=b_{4} \end{matrix}\right.\]

この連立1次方程式は係数行列\(\mathbf{A}\)と定数ベクトル\(\mathbf{b}\)を使うと、以下のように表現することが可能です。

\[
\mathbf{A}=\begin{bmatrix}
a_{11}& a_{12}& a_{13}& a_{14}\\
a_{21}& a_{22}& a_{23}& a_{24}\\
a_{31}& a_{32}& a_{33}& a_{34}\\
a_{41}& a_{42}& a_{43}& a_{44}\\
\end{bmatrix},
\mathbf{b}=\begin{bmatrix}
b_{1}\\
b_{2}\\
b_{3}\\
b_{4}
\end{bmatrix}
\]

また、拡大係数行列として以下の表記も可能です。

\[\mathbf{[Ab]}=\begin{bmatrix} a_{11}& a_{12}& a_{13}& a_{14}&b_{1}\\ a_{21}& a_{22}& a_{23}& a_{24}&b_{2}\\ a_{31}& a_{32}& a_{33}& a_{34}&b_{3}\\ a_{41}& a_{42}& a_{43}& a_{44}&b_{4}\\ \end{bmatrix}\]

ガウスの消去法はこの拡大係数行列に対し、前進消去後退代入という2つのステップを使って連立方程式の解を機械的に求めます。

前進消去

前進消去(Forward elimination)とは、先ほどの拡大係数行列に対し行基本操作を行い、係数行列の対角成分を1に、下三角を0に変形して階段行列(0が階段のように見えることが由来)を作ることです。

前進消去の説明図

行基本操作とは、
①ある行を定数倍する。
②ある行を定数倍したものを他の行に加える。
③行を入れ替える。

といった操作で、これらの操作は方程式の意味を変えない変形であるため同値変形と呼ばれます。

学校で最初に連立方程式を教わった時に学ぶ方法ですね!この同値変形は方程式がいくつになろうと成立します。

後退代入

後退代入(Backward substitution)とは、前進消去で作り上げた階段行列に対して行う操作で、連立1次方程式の求解部分のことです。

階段行列の最下段は係数行列の左側が全て0、右端が1になっています。その右隣は定数項なので、既に解が得られている段となります。

最下段の解がわかれば、その1つ上の段の未知数も計算可能になるため、以下の図のようにドミノ倒し的に解を求めることができます。下から代入していって解を求めることから、後退代入という名称が付いています。

後退代入の説明図

よく使われる所

教科書や受験レベルの数学問題は手計算によって解くことが前提の場合が多いため、多くの人にとってはあまりガウスの消去法を使うことは無いと思われます。

ですが、実は大量の連立1次方程式を解く問題は世の中に沢山あり、その代表的な例がコンピュータシミュレーションの分野となります。

有限要素法(FEM:Finite Element Method)は特に有名で、構造物を有限の要素に離散化して連立方程式を立てます。

その方程式の数は100万、1000万はざらで、とても人間が手計算で解くべきものではありません。

ガウスの消去法は機械的に(プログラム的に)連立1次方程式を計算することができるため、コンピュータシミュレーションの分野では基本の手法となっています。

手計算で解いてみよう!

Pythonプログラミングに行く前に、簡単な例題を手計算で解いてみましょう。

問題

例題として以下の4元連立1次方程式を用意しました。

\[
\left\{\begin{matrix}
x_{1}-x_{2}-2x_{3}+2x_{4}=5\\
2x_{1}-x_{2}-3x_{3}+3x_{4}=10\\
-x_{1}+3x_{2}+3x_{3}-2x_{4}=2\\
x_{1}+2x_{2}-x_{4}=-10
\end{matrix}\right.
\]

拡大係数行列は以下の\([\mathbf{Ab}]\)になります。

\[
[\mathbf{Ab}]=\begin{bmatrix}
1&-1&-2&2&5\\
2&-1&-3&3&10\\
-1&3&3&-2&2\\
1&2&0&-1&-10\\
\end{bmatrix}
\]

解法(前進消去)

前進消去の第1ステップとして、\(a_{11}\)のピボットと呼ばれる成分で、その行の全ての要素を割ります(今回は\(a_{11}\)が既に1なので特に何もしませんが…)。

ガウスの消去法プロセスのstep1

次に2行目以降について、先頭列(\(a_{11}\)ピボットと同じ列)を0にするよう、\(a_{11}\)行を定数倍して引くという同値変形を行います

ガウスの消去法プロセスのstep2

次に\(a_{22}\)をピボットとして、先ほどと同様にピボットの成分で割ります(といってもここもピボットが1なので変化はありませんが…)。

ガウスの消去法プロセスのstep3

続いて3行目以降について、\(a_{22}\)ピボットと同じ列を0にするよう同値変形を行います。

ガウスの消去法プロセスのstep4

ここまで来れば「もうわかった」と言われそうですが、1つずつ説明していきます。

3行目もピボットで割り、\(a_{33}\)を1にします。

ガウスの消去法プロセスのstep5

そしてピボット行の定数倍を引く同値変形を行い…

ガウスの消去法プロセスのstep6

最終行の\(a_{44}\)も1にするようにピボットで割ると…係数行列の下三角が0になる階段行列を得ることが出来ました。

ガウスの消去法プロセスのstep7

解法(後退代入)

後退代入は前進消去よりさらに簡単です。最下段は、

\[x_{4}=11\]

を意味しているので、3行目は、

\[x_{3}=-7+2x_{4}=15\]

となります。
同様に2行と1行目は、

\[x_{2}=-x_{3}+x_{4}=-4\\
x_{1}=5+x_{2}+2x_{3}-2x_{4}=9\]

と、全ての解\(x_{1}=9, x_{2}=-4, x_{3}=15, x_{4}=11\)を得ることが出来ました。

ガウスの消去法(掃き出し法)は手計算でやってみるのが最速の理解法と思います。是非手を動かしてみて下さい。

行列の基本操作があやしいという方は、以下の「技術者のための線形代数学」をオススメします。数学書にありがちな難しい言い回しはなく、数学を専門としていない人のために平易な言葉で書いてあり、入門書としては読みやすいです。

中井悦司著【大学の数学を本気で学ぶ】

Python/Numpyによるガウスの消去法プログラム

基本コード

以下にPython3.7でコーディングしたガウスの消去法プログラムを紹介します。
詳細説明はコード内にコメントを記載しましたので、そちらを参照下さい。

実は以下の書籍「解析塾秘伝 有限要素法のつくり方」のP76にあるVBAコードを参考にしました。

石川博幸, 青木伸輔, 日比学 著【FEMプログラミングの手順とノウハウ】

書籍内ではExcel VBAで書かれていますが、Pythonだと行の割り算でforループを使わなくても良いといった違いがありましたので、Python用にいくつか変更しています。

本記事手計算の説明では拡大係数行列を作ってから計算をしていましたが、プログラム的には定数ベクトルはそのままにしておいた方が後退代入がしやすいというメリットもあるようです。

このプログラムを実行すると、すべての解がコンソールに表示されます。

Pivotエラー付きガウスの消去法プログラム

先ほどのプログラムで使用していた行列は問題ありませんでしたが、pivotの値が0になってしまうことはよくあると思います。そのため、まずはエラーが発生したら正常にエラー終了させるように先ほどのコードを改修してみました。

エラー処理はtryexceptで書きました。if文でpivotの値を調べ、閾値(ここでは1e-6)未満であればraiseで強制的に0割エラー(ZeroDivisionError)フラグを立たせています。このエラーフラグが立った時点でfor文は実行されなくなり、エラー文を出力させて関数を終了させます。

以下が実行結果です。上記コードでは最初の\(a_{11}\)を0と変更しエラーとなるようにしてみました。 (「Error termination」は商用コードである某LS-DYNAをイメージしてみました。)

計算が正常終了すると安心の「Normal termination」文が表示されるので、是非お試し下さい。

部分ピボット選択付きガウスの消去法プログラム

先ほどの0割エラーは「部分ピボット選択」という手法を使えばほとんどの場合回避可能です。

手計算では無意識にやることだと思いますが、連立方程式なので拡大係数行列の行を入れ替える操作も同値変形となり、方程式の意味を崩しません。

以下のコードはpivotに値を格納する前に、pivot列のpivot行以降で絶対値が最大になる行を常にpivotに使うよう上記コードを改修したものになります。

今回はAとbに分けて書いているので、入れ替えはそれぞれで行う必要がありますが、これで先ほどの連立方程式もNormal terminationになりました。
(Aとbの中の行を入れ替えて検算はしてみましたが、もしお使いになる場合はご自身で検証下さい)

コード検証方法

本当に計算結果があっているのかを検証するには、最後に係数行列Aと解xの内積を計算するのが良いと思います。結果が定数ベクトルbになっていれば正解。

メイン部分のみですが、以下に検証のためのコードを示します。
本記事のプログラムはAとbの中身が計算中に変わってしまうため、検証のためにはnp.copy()で値を保持しておく必要があります。

また、このコードではデータ型を「dtype=float」としています。これをつけ忘れると値が変わってしまうのでご注意(関数の中で処理すれば良いかも?)。

上記コードを追記して実行すると、最後に内積(A@b)とbが出力されます。この値がイコールになっていればOK。

連立方程式を解くその他のPythonコード

Pythonで線形代数!連立1次方程式を解く」では、逆行列を使ったり、Pythonのsolveを使ったりした連立方程式の解き方も紹介していました。今回はFEMでよく使われている方法であるガウスの消去法を理解するためにNumpyで書いてみましたが、ご興味があればこちらの記事もご覧ください。

まとめ

本記事ではガウスの消去法(畳み込み法)の概要として、前進消去と後退代入を説明しました。そして実際に簡単な問題を手計算でステップを追いながらアルゴリズムを体感しました。

最後に、PythonのNumpyのみでガウスの消去法コードを書き、プログラム的に連立1次方程式を解く方法を紹介しました。プログラムは基本コードに加え、ピボットエラー検出コード、部分ピボット選択付きコードと計3例作ってみました。

補足

連立方程式と解について

1つ補足します。連立1次方程式の解の個数は係数行列\(\mathbf{A}\)と定数ベクトル\(\mathbf{b}\)の状態によって変わります。

以下は直線の方程式の場合の図ですが、\(\mathbf{A}\)が正則行列であれば2直線の交点というただ1つの解が存在することになります。

連立方程式の解が1つに決まる時

また、以下のように2直線が平行である場合は\(\mathbf{A}\)が正則行列にならず、定数ベクトル\(\mathbf{b}\)によって解の個数が変わります。

連立方程式の解が無い場合と無数にある場合

直線が離れている場合は交点が無いので解無しとなります。一方、直線が一致している場合は全ての点が解となり、無限に解が存在することを意味します。

もし今回のプログラムでエラーばかり出てしまった場合は、係数行列が正則かどうかを調べることをオススメします。

数値計算の基礎中の基礎、ガウスの消去法をプログラミングしてみました!今後は関連記事を増やす予定です!
Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています!

  • このエントリーをはてなブックマークに追加

SNSでもご購読できます。

コメントを残す

*