今日では様々な科学技術計算に「行列」が使われています。行列は線形代数学の分野です。行列の計算は手計算だと少々やっかいですが、退屈な単純計算はPythonに任せましょう。ここでは線形代数による連立1次方程式の解法をPythonでコーディングします。
こんにちは。wat(@watlablog)です。
真面目に科学技術計算をしようと思ったら線形代数の理解が不十分とわかりました。ここでは線形代数を利用した連立1次方程式の解法をPythonで学びます!
僕は線形代数の勉強を「中井悦司, 技術者のための線形代数学, 翔泳社※Amazon」でしています。書籍ではプログラミングについては触れていませんが、ここではPythonによるコーディングにフォーカスして学んだ事の実践を行います。
線形代数学の内容について、より理解を深めたい方は本書の購入をすることで、当ブログで紹介するPythonコードの意味がよくわかるようになると思います。
当ブログでは、基本的に難しい説明は面倒なのでできるだけしないようにし、結果がPythonで簡単に得られることを優先するという方針で進めます!
連立1次方程式の解法
方程式の行列表現
連立方程式とは、以下の式(1)に示す未知数\(x, y\)について記述した方程式が連なっているものを指します。
\[ \begin{cases}
ax+by=p \\
cx+dy=q
\end{cases} (1)
\]
この式(1)を解くためには、中学校では筆算の形式にしたり、代入の方法をとったりして解いた記憶のある方も沢山いらっしゃると思います。
しかし、線形代数として行列を覚えた人は式(1)を式(2)と書き換えて計算を行うことができます。
\[
\begin{bmatrix}
a &b \\
c &d
\end{bmatrix}
\begin{bmatrix}
x\\
y\end{bmatrix}
=
\begin{bmatrix}
p\\
q
\end{bmatrix} (2)
\]
式(2)は、式(1)の左辺\(x\)と\(y\)にかかっている係数を行列形式で記述した係数行列を作っています。この係数行列を行列\([x, y]\)にかけていますが、この\(x, y\)が入った行列を変数ベクトルと呼びます。そして定数\(p, q\)が入った右辺の行列を定数ベクトルと呼びます。
もっと式を一般化すると、変数がいくつあっても式(3)のように行列のサイズをどんどん拡大していくことで表現することが可能です。
ある変数が無い式があっても、そこには0という係数がかかっていると考えれば行列が成立します。
\[
\begin{bmatrix}
a_{11} &a_{12} &\cdots &a_{1n} \\
a_{21} &a_{22} &\cdots &a_{2n} \\
\vdots &\vdots &\ddots &\vdots \\
a_{n1} &a_{n2} &\cdots &a_{nn}
\end{bmatrix}
\begin{bmatrix}
x_{1}\\
x_{2}\\
\vdots \\
x_{n}\end{bmatrix}
=
\begin{bmatrix}
b_{1}\\
b_{2}\\
\vdots \\
b_{n}
\end{bmatrix} (3)
\]
係数行列をA, 変数ベクトルをx, 定数ベクトルをbと書けば、式(3)は式(4)と書くことができます。これらは全て行列です。
$$\mathbf{A}\mathbf{x}= \mathbf{b} (4)$$
式(4)は変数ベクトルとして解を出したい変数リストが左辺にまとまっているので、式(5)と定数行列bに左からAの逆行列をかけることで等式が成立して解を求めることができます。
$$\mathbf{x}= \mathbf{A}^{-1}\mathbf{b} (5)$$
通常未知数の数だけ方程式を用意する問題を扱う場合がほとんどだと思いますので、ここでは係数行列が正方形になる正方行列しか扱いません(ある式に変数が無かったら、0を使って係数行列を作ると解決します)。
これで行列を使った連立方程式の解法のイメージが固まったと思います。早速Pythonでコーディングして確かめてみましょう!
Pythonによる連立1次方程式の解法コード
まじめに解く
式(5)の関数を作成して、まじめに解いてみましょう。係数行列Aと定数行列bには適当に値を入れています。
逆行列の計算だけは「Python/NumPyで線形代数!linalgで逆行列を求める方法」で使ったnp.linalg.inv関数で手抜きをしています。行列の演算をするためにnp.dot関数を使っているのも、前回の記事で紹介した通りです。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
import numpy as np def sol(A, b): inv_A = np.linalg.inv(A) # 逆行列 sol = np.dot(inv_A, b) # 行列の乗算 return sol # 係数行列 A = np.array([[1, 2], [3, 5]]) # 定数ベクトル b = np.array([4, 6]) sol_manual = sol(A, b) print(sol_manual) |
結果、「-8, 6」という解を得ます。Pythonで一度コーディングしておけば、連立1次方程式なんて一瞬で解けますね!
solve関数でもっと手を抜く
上記方法でもかなり簡単に解くことができましたが、まだまだ先人の知恵を活用して巨人の肩に乗ることができます。
PythonのNumPyには逆行列の時に使ったlinalgの中に、solve関数というものがあり、この関数を使うとさらにコードの行数を減らして解を求めることが出来ます。
1 2 3 4 5 6 7 8 9 10 |
import numpy as np # 係数行列 A = np.array([[1, 2], [3, 5]]) # 定数ベクトル b = np.array([4, 6]) sol_auto = np.linalg.solve(A, b) print(sol_auto) |
上記コードを実行しても、「-8, 6」という同様の解を得ることができます。
まとめ
本ページでは様々な所で活用されている線形代数を使った連立方程式の解法をPythonで実装してみました。
連立方程式は行列形式で記述することで、非常にエレガントな表記ができ、理解の助けにもなるし計算も簡便になります。
今回は高校生や大学生初年度クラスの内容でしたが、Pythonでコーディングすることで面倒な行列の演算からくるストレスから解放されます!
ストレスになる計算はPythonにやらせましょう!
Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています!
コメント