Python NumPy SciPy サンプルコード: 線形連立方程式, LU 分解

これから何回かに分けて、 線形問題 の解法についてまとめていきます。

Python で線形方程式系を解く

Python で正方系の線形方程式系(線形連立方程式)の数値解を求めてみます。

行列表現の線形方程式系を解く

簡単な以下の例で Python による線形方程式系の解を求めます。

\begin{eqnarray} \left\{ \begin{array}{l} 6x_1 + 4x_2 + 1x_3 = 7 \\ x_1 + 8x_2 - 2x_3 = 6 \\ 3x_1 + 2x_2 = 8 \end{array} \right. \end{eqnarray}

これを行列表現で表すと、

\begin{equation} \mathbf{A} = \begin{bmatrix} {6} & {4} & {1}\\ {1} & {8} & {-2}\\ {3} & {2} & {0} \end{bmatrix}, \boldsymbol{x} = \begin{pmatrix} {x_1}\\ {x_2}\\ {x_3} \end{pmatrix}, \boldsymbol{b} = \begin{pmatrix} {7}\\ {6}\\ {8} \end{pmatrix} \end{equation}
\begin{equation} \mathbf{A}\boldsymbol{x} = \boldsymbol{b} \end{equation}

となり、この例では \(\mathbf{A}\) が正則で逆行列を持つので、解は以下のように求めることができます。

\begin{equation} \boldsymbol{x} = \mathbf{A}^{-1}\boldsymbol{b} \end{equation}

解を計算するプログラムは、NumPy の配列を用いれば次のように書くことができます。

import numpy as np

A = np.array([[6, 4, 1],
              [1, 8, -2],
              [3, 2, 0]])
b = np.array([7, 6, 8])

Ainv = np.linalg.inv(A)         # A の逆行列を計算
x = np.dot(Ainv, b)             # A の逆行列と b の内積

print(x)
[ 4. -2. -9.]

ただ、NumPy には線形方程式系を解くメソッドがあるので、通常はそれを使います(こちらの方が計算量が少なくて速い)。

import numpy as np

A = np.array([[6, 4, 1],
              [1, 8, -2],
              [3, 2, 0]])
b = np.array([7, 6, 8])

x = np.linalg.solve(A, b)

print(x)
[ 4. -2. -9.]

LU 分解

NumPy と SciPy の linalg.solve 関数は LAPACK の GESV ルーチンを使用しており、LU 分解による方法で効率的に解が計算されます。

\(\mathbf{A}\) が同じで、\(\boldsymbol{b}\) が異なるような方程式系を解く場合は、lu_factor 関数と lu_solve 関数を使用すれば効率良く計算することができます。 lu_factor 関数は \(\mathbf{L}+\mathbf{U}\) から \(\mathbf{L}\) の単位対角成分を除いた行列と、置換行列を表す指数を返します。 こうして求めた行列を使い回せば、計算時間を短縮することができます。

import scipy.linalg as linalg
import numpy as np

A = np.array([[6, 4, 1],
              [1, 8, -2],
              [3, 2, 0]])
b = np.array([7, 6, 8])

LU = linalg.lu_factor(A)
x = linalg.lu_solve(LU, b)

print(x)
[ 4. -2. -9.]

コメント

Comments powered by Disqus
書籍更新情報
2017-02-18
Pythonによる科学技術計算 基礎編
1.4版への更新が可能になりました。
サポートページはこちら
電子書籍
Pythonによる科学技術計算 基礎編
電子書籍
線形代数(1): Pythonによる科学技術計算 実践編