Python (NumPy) で方程式・連立方程式を解く

プログラミング言語 Pythonを使って方程式・連立方程式を解いてみたいと思います。

今回は数値計算ライブラリNumPyを使って数値的に解く方法をみていきます。

代数的に厳密に解く方法は、数式処理ライブラリSymPyを使っている次の記事

pianofisica.hatenablog.com

で紹介しています。

数学定数の値
代数方程式の数値解
連立１次方程式の数値解
行列の数値的取り扱い

数学定数の値

まずNumPyに慣れるためにいくつかの数学定数をみてみましょう：

import numpy
numpy.sqrt(2)

$~\sqrt{2}=1.4142135623730951$

numpy.e

$~e=2.718281828459045$

numpy.pi

$~\pi=3.141592653589793$

numpy.log(2)

$~\log 2=0.6931471805599453$

もちろんいずれも無限小数ですから誤差を含んでいることには注意が必要です：たとえば

numpy.sqrt(2)**2

$~\sqrt{2}^2\overset{?}{\sim}2.0000000000000004$

と計算されて、最後の桁に誤差が現れていることがわかります。

以下ではその都度numpyとつけるのが面倒なので

import numpy as np

と省略しておきます：例えば

np.pi

代数方程式の数値解

代数方程式

${\displaystyle{~a_1 x^n+a_2x^{n-1}+\dots+a_{n-1}x+a_{n}=0}}$

の解をPythonで数値的に求めるには

import numpy as np
np.roots([a_1, a_2, ..., a_n-1, a_n])

と入力します。

例えば２次方程式

$~x^2+2x+3=0$

を考えてみましょう。解の公式から

$~x=-1\pm i\sqrt{2}$

ですが、実際

np.roots([1, 2, 3])

と入力することにより

$~x_1=-1.+1.41421356i \quad ~x_2 =-1.-1.41421356i$

が出力され、確かに数値解が得られていることがわかります。

また別な例としては、一般に解の公式が存在しない５次方程式

$~x^5+6x^4+x^3+x+1=0$

を考えてみましょう。

この方程式は数式処理をするSymPyでは解くことができません：

import sympy
sympy.var('x')
sympy.solve (x**5+6*x**4+x**3+x+1, x)

ところがNumPyを用いると

np.roots([1, 6, 1, 0, 1, 1])

によって、その数値解を

$~x_1=-5.8241071 \\ ~x_2 =0.42084104+0.52540961\, i \\ ~x_3= 0.42084104-0.52540961\,i \\ ~x_4=-0.50878749+0.34645119\,i \\~x_5=-0.50878749-0.34645119\,i$

と求めることができます。

<br />

連立１次方程式の数値解

（独立な）連立１次方程式は原理的に厳密に解を求めることができるので、

ことさら数値解と強調する必要はないかもしれませんが、

ここではNumPyを用いて解く方法としてまとめておきます。

連立１次方程式

$\displaystyle{~\begin{pmatrix} a_{11} & a_{12} & \dots & a_{1n} \\ a_{21} & a_{22} & \dots & a_{2n} \\ \vdots & \vdots & \dots & \vdots \\ a_{n1} & a_{n2} & \dots & a_{nn} \end{pmatrix}\begin{pmatrix} x_1 \\ x_2 \\ \dots \\ x_n \end{pmatrix}=\begin{pmatrix} y_1 \\ y_2 \\ \dots \\ y_n \end{pmatrix}}$

の解は、行列

$\displaystyle{~A=\begin{pmatrix} a_{11} & a_{12} & \dots & a_{1n} \\ a_{21} & a_{22} & \dots & a_{2n} \\ \vdots & \vdots & \dots & \vdots \\ a_{n1} & a_{n2} & \dots & a_{nn} \end{pmatrix}}$

の逆行列 $A^{-1}$ （いま逆行列は存在するものと仮定します）を使って

$\displaystyle{~\begin{pmatrix} x_1 \\ x_2 \\ \dots \\ x_n \end{pmatrix}=A^{-1}\begin{pmatrix} y_1 \\ y_2 \\ \dots \\ y_n \end{pmatrix}}$

として求められます。

これをNumPyで計算するには、まず行列を

A = np.matrix([
[a_11, a_12, ..., a_1n],
[a_21, a_22, ..., a_2n],
[...., ...., ...., ....]
[a_n1, a_n2, ..., a_nn]
])

で入力し、その逆行列を

np.linalg.inv(A)

によって求めます。そして列ベクトル $Y$ を

Y = np.matrix([
[y_1],
[y_2],
[...],
[y_n]
])

によって入力し、 $A^{-1}$ と掛け合わせます：

np.linalg.inv(A)*Y

これで連立１次方程式の解が得られたことになります。

例えば、２変数の連立１次方程式

$~2x+y=3 \\ ~x+3y=4$

に対しては

A = np.matrix([
[2, 1],
[1, 3]
])
Y = np.matrix([
[3],
[4]
])
np.linalg.inv(A)*Y

より解

$~x=1.0, \quad y=1.0$

を得ます。または

np.linalg.solve(A,Y)

と入力することによっても解を得ることができます。

行列の数値的取り扱い

連立１次方程式は行列の理論（線型代数学）と密接に関係しています。NumPyによる行列のとりあつかいについては次の記事でまとめて解説しています。

pianofisica.hatenablog.com

本記事の基本的な内容としては、Maximaで同様のことをしている過去記事

pianofisica.hatenablog.com

の中のMaximaの入力を、PythonのNumPyライブラリの書法で置き換えたものになっています。

あわせて読んでみると、PythonとMaximaの違いを見比べられるかと思います。

キーワード：Python、NumPy、JupyterNotebook