pianofisica

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

Python (SymPy) による行列の計算

数学の具体的な計算にPythonを使って、数学もPythonも同時に学んでしまいましょう。

今回はPythonを使って行列の計算をしてみたいと思います。

Pythonのごく基本的な使い方については以下の記事を参照してください:

pianofisica.hatenablog.com



行列の諸操作

行列を入力する

Python (SymPy) で行列を入力するには

import sympy
sympy.init_printing()
sympy.var('m11, m12, n11, n21')
sympy.var('a11, a12, a21, a22, b11, b12, b21, b22') 
M = sympy.Matrix([
[m11,m12]
])
N = sympy.Matrix([
[n11],
[n21]
])
A = sympy.Matrix([
[a11,a12],
[a21,a22]
])
B = sympy.Matrix([
[b11,b12],
[b21,b22]
])

などとすれば良いです。

ここで  M (1,2)-行列、 N (2,1)-行列、 A B (2,2)-行列としました:

 \displaystyle{M=\begin{bmatrix} m_{11} & m_{12}  \end{bmatrix} \qquad N=\begin{bmatrix} n_{11} \\ n_{21}  \end{bmatrix} \qquad A=\begin{bmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{bmatrix} \qquad B=\begin{bmatrix} b_{11} & b_{12} \\ b_{21} & b_{22} \end{bmatrix}}

行列の要素を参照する

行列の要素を参照するには

A[0,1]   #  つまり行列 A の (k ,l) 成分は A[ k-1 , l-1 ] です(ラベルは0から数えます)
B[1,1]

などとします。

行列の演算(積、和、定数倍、ベキ乗)

行列の積の定義から、行列の積  AB MA BN MBN などが許されます:

A*B    #  行列の積はアスタリスク(通常の積と同じ)です
M*A
B*N
M*B*N

などとして計算します。行列のベキ乗は

A**2

です。

同じ型の2つの行列の線形結合は

sympy.var('a, b')
a*A+b*B

となります。

逆行列をとる

逆行列

A.inv()

または

A**(-1)

によって求められます。

連立1次方程式を解く

求めた逆行列を応用する問題として、連立1次方程式

 ~ a_{11}x_1 + a_{12} x_2 + \dots + a_{1n} x_n=y_1 \\ ~ a_{21}x_1 + a_{22} x_2 + \dots + a_{2n} x_n=y_2 \\ ~~ \vdots \\ ~ a_{n1}x_1 + a_{n2} x_2 + \dots + a_{nn} x_n=y_n

を考えましょう。これは行列を使って

 \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_n=\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_n^{-1}(いま逆行列は存在するものと仮定します)を使って

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

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

例えば3変数  x y z の連立1次方程式

 ~ ax+by+cz=j \\ ~ dx+ey+fz=k \\ ~ gx+hy+iz=l

は行列

 \displaystyle{~A_3=\begin{pmatrix} a & b & c \\ d & e & f \\ g & h & i \end{pmatrix}}

逆行列をもとめれば良いわけです:

import sympy
sympy.init_printing()
sympy.var('a, b, c, d, e, f, g, h, i') 
A3 = sympy.Matrix([
[a, b, c],
[d, e, f],
[g, h, i]
])
sympy.cancel(A3**(-1))

により、逆行列

 \displaystyle{~A_3^{-1}=\frac{1}{\det A_3}\begin{pmatrix} ei-fh & -(bi-ch) & bf-ce \\ -(𝑑𝑖−𝑓𝑔) & 𝑎𝑖−𝑐𝑔 & -(𝑎𝑓−𝑐𝑑) \\ 𝑑ℎ−𝑒𝑔 & -(𝑎ℎ−𝑏𝑔) & 𝑎𝑒−𝑏𝑑 \end{pmatrix}}

と求まります。ここで

 \displaystyle{~\det A_3=a(ei-fh)-b(di-fg)+c(dh-eg)}

です。最後はこの逆行列にベクトル

 \displaystyle{~Y=\begin{pmatrix} j \\ k \\ l \end{pmatrix}}

を掛け合わせれば解が得られます:

 \displaystyle{~\begin{pmatrix} x \\ y \\ z \end{pmatrix} = A_3^{-1}Y= \frac{1}{\det A_3}\begin{pmatrix} j(ei-fh)  -k(bi-ch) +l(bf-ce) \\ -j(𝑑𝑖−𝑓𝑔) +k( 𝑎𝑖−𝑐𝑔)  -l(𝑎𝑓−𝑐𝑑) \\ j(𝑑ℎ−𝑒𝑔)  -k(𝑎ℎ−𝑏𝑔) +l( 𝑎𝑒−𝑏𝑑) \end{pmatrix}}

これをsolveを使って求めた解と見比べてみましょう:

eq1=sympy.Eq(a*x+b*y+c*z, j)
eq2=sympy.Eq(d*x+e*y+f*z, k)
eq3=sympy.Eq(g*x+h*y+i*z, l)
sympy.solve ([eq1, eq2, eq3], [x, y, z]) 

より

 \displaystyle{x={{j(ei-fh)-k(bi-ch)+l(bf-ce)}\over{a(ei-fh)-b(di-fg)+c(dh-eg)}} \\ y={{−𝑗(𝑑𝑖−𝑓𝑔)+𝑘(𝑎𝑖−𝑐𝑔)−𝑙(𝑎𝑓−𝑐𝑑)}\over{a(ei-fh)-b(di-fg)+c(dh-eg)}} \\ z={{𝑗(𝑑ℎ−𝑒𝑔)−𝑘(𝑎ℎ−𝑏𝑔)+𝑙(𝑎𝑒−𝑏𝑑)}\over{a(ei-fh)-b(di-fg)+c(dh-eg)} }}

確かに解になっていることがわかります。

行列式を求める

行列  A行列式を求めるには

A.det()

と入力します。

行列  A の余因子行列  \widetilde{A}

A.adjugate()

です:実際

A.inv() == 1/(A.det())*A.adjugate()

より

 \displaystyle{\widetilde{A}=\left(\det A\right) A^{-1}}

が成り立ちます。

単位行列

単位行列

sympy.eye(2)

です。引数の数字は行列のサイズを指定します。

たとえば、逆行列になっていることを確認するときなどに

sympy.simplify(A.inv()*A) == sympy.eye(2)

のようにして使うことができます。

トレースをとる

行列  A のトレースをとるには

A.trace()

とします。

転置をとる

行列  A の転置行列

 \displaystyle{A^t=\begin{bmatrix} a_{11} & a_{21} \\ a_{12} & a_{22} \end{bmatrix}}

A.transpose()

です。

随伴行列をとる

行列  A の随伴行列

 \displaystyle{A^\dagger=\begin{bmatrix} \overline{a_{11}} & \overline{a_{21}} \\ \overline{a_{12}} & \overline{a_{22}} \end{bmatrix}}

をとるには

A.adjoint()

とします。

行列を対角化する

行列の固有値を求める
A.eigenvals()
行列の固有ベクトルを求める
A.eigenvects()
対角行列を求める
A.diagonalize()


キーワードPython、SymPy、行列