pianofisica

Mathematics & Physics, Maxima, a bit Python & Wolfram, and Arts

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])

と入力します。

例えば2次方程式

 ~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

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

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

 ~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

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

連立1次方程式の数値解


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

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

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

連立1次方程式

 \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

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

例えば、2変数の連立1次方程式

 ~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)

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


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

pianofisica.hatenablog.com

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

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


キーワードPython、NumPy、JupyterNotebook