Pythonを用いた連立１次方程式の数値解法（LU分解） - モデリングと数値計算（主に化学工学向け）

はじめに

Python (Numpy, Scipy) を用いて連立１次方程式を数値的に解く方法についてメモを記す。

例題

$\begin{aligned} &x+2y-z=2 \\ &2x-y-z=-3 \\ &5x+y-3z=-2 \end{aligned}$

これを行列で表すと、

$\mathbf{A}= \begin{bmatrix} 1 & 2 & -1 \\ 2 & -1 & -1 \\ 5 &1 &-3 \\ \end{bmatrix}, \quad \mathbf{x}= \begin{bmatrix} x \\ y \\ z \\ \end{bmatrix}, \quad \mathbf{b}= \begin{bmatrix} 2 \\ -3 \\ -2 \\ \end{bmatrix}$

$\mathbf{A} \mathbf{x} = \mathbf{b}$

解き方

numpy.linalg.solveを使って次のように計算できる。

import numpy

A = numpy.array([[1, 2, -1],
                 [2, -1, -1],
                 [5, 1, -3]])

b = numpy.array([2,
                 -3,
                 -2])

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

print x

結果は次のとおりとなる。

array([1., 2., 3.])

1組の連立方程式を単発で解く場合は上の方法で問題ないが、 $\textbf{A}$ の部分は変わらずに $\textbf{b}$ が変わる場合を数多く解く場合には次のLU分解を使うやり方のほうが計算時間を節約できる。一度計算して求めたlu_factorを使いまわして $\textbf{b}$ が変わった場合を計算すればよい。

import numpy
import scipy.linalg

A = numpy.array([[1, 2, -1],
                 [2, -1, -1],
                 [5, 1, -3]])

b = numpy.array([2,
                 -3,
                 -2])

lu = scipy.linalg.lu_factor(A)

x = scipy.linalg.lu_solve(lu, b)