どうも、木村(@kimu3_slime)です。
今回は、Julia(SymPy)で行列式、逆行列を求める方法を紹介します。
準備
SymPyを使うので、持っていなければインストールしておきましょう。
1 2 | using Pkg Pkg.add("SymPy") |
準備として、以下のコードを実行しておきます。
1 | using SymPy |
行列式、逆行列を求める方法
行列式は「行列.det()」で、逆行列は「行列^-1」または「行列.inv()」で求められます。
1 2 3 4 | B = Sym[1 2; 3 4] B.det() B^-1 B.inv() |
\[ \begin{aligned}\left[ \begin{array}{rr}1&2\\3&4\end{array}\right]\end{aligned} \]
\[ \begin{aligned}-2\end{aligned} \]
\[ \begin{aligned}\left[ \begin{array}{rr}-2&1\\\frac{3}{2}&- \frac{1}{2}\end{array}\right]\end{aligned} \]
\[ \begin{aligned}\left[ \begin{array}{rr}-2&1\\\frac{3}{2}&- \frac{1}{2}\end{array}\right]\end{aligned} \]
行列式の余因子展開に伴って現れる、行列の余因子行列\(\mathrm{adj\,}(B)\)を求めることもできます。
1 | B.adjugate() |
\[ \begin{aligned}\left[ \begin{array}{rr}4&-2\\-3&1\end{array}\right]\end{aligned} \]
余因子行列は、
\[ \begin{aligned}B \mathrm{adj\,}(B)= (\det B) I\end{aligned} \]
\[ \begin{aligned}B^{-1} = \frac{1}{\det B} \mathrm{adj\,}(B)\end{aligned} \]
といった性質を満たします。これを計算して確かめてみましょう。「sympy.eye(2)」はサイズ2の単位行列です。
1 2 | B*B.adjugate() == B.det() *sympy.eye(2) B^-1 == 1/B.det()*B.adjugate() |
1 2 | true true |
ここからは、成分に記号を使います。
1 | @vars a11 a12 a13 a21 a22 a23 a31 a32 a33 x y z n α β |
回転行列の行列式、逆行列を求めてみましょう。
1 2 3 | R_α= Sym[cos(α) -sin(α); sin(α) cos(α)] R_α.det() R_α.inv() |
\[ \begin{aligned}\left[ \begin{array}{rr}\cos{\left(α \right)}&- \sin{\left(α \right)}\\\sin{\left(α \right)}&\cos{\left(α \right)}\end{array}\right]\end{aligned} \]
\[ \begin{aligned}\sin^{2}{\left(α \right)} + \cos^{2}{\left(α \right)}\end{aligned} \]
\[ \begin{aligned}\left[ \begin{array}{rr}- \frac{\sin^{2}{\left(α \right)}}{\cos{\left(α \right)}} + \frac{1}{\cos{\left(α \right)}}&\sin{\left(α \right)}\\- \sin{\left(α \right)}&\cos{\left(α \right)}\end{array}\right]\end{aligned} \]
結果がシンプルでないので、「simplify」で単純化します。
1 2 | simplify(R_α.det()) sympy.simplify(R_α.inv()) |
\[ \begin{aligned}1\end{aligned} \]
\[ \begin{aligned}\left[\begin{smallmatrix}\cos{\left(α \right)} & \sin{\left(α \right)}\\- \sin{\left(α \right)} & \cos{\left(α \right)}\end{smallmatrix}\right]\end{aligned} \]
角度\(\alpha\)によらず行列式が1で、\((R_{\alpha})^{-1}= R_{-\alpha}\)であることがわかりました。
成分が一般的な3次正方行列の行列式、逆行列を求めてみましょう。
1 2 3 | A = [Sym("a$i$j") for i=1:3 for j=1:3].reshape(3,3) A.det() A^-1 |
\[ \begin{aligned}\left[ \begin{array}{rrr}a_{11}&a_{12}&a_{13}\\a_{21}&a_{22}&a_{23}\\a_{31}&a_{32}&a_{33}\end{array}\right]\end{aligned} \]
\[ \begin{aligned}a_{11} a_{22} a_{33} – a_{11} a_{23} a_{32} – a_{12} a_{21} a_{33} \\+ a_{12} a_{23} a_{31} + a_{13} a_{21} a_{32} – a_{13} a_{22} a_{31}\end{aligned} \]
\[ \begin{aligned}\left[ \begin{array}{rrr}\frac{a_{22} a_{33} – a_{23} a_{32}}{a_{11} a_{22} a_{33} – a_{11} a_{23} a_{32} – a_{12} a_{21} a_{33} + a_{12} a_{23} a_{31} + a_{13} a_{21} a_{32} – a_{13} a_{22} a_{31}}&\frac{- a_{12} a_{33} + a_{13} a_{32}}{a_{11} a_{22} a_{33} – a_{11} a_{23} a_{32} – a_{12} a_{21} a_{33} + a_{12} a_{23} a_{31} + a_{13} a_{21} a_{32} – a_{13} a_{22} a_{31}}&\frac{a_{12} a_{23} – a_{13} a_{22}}{a_{11} a_{22} a_{33} – a_{11} a_{23} a_{32} – a_{12} a_{21} a_{33} + a_{12} a_{23} a_{31} + a_{13} a_{21} a_{32} – a_{13} a_{22} a_{31}}\\\frac{- a_{21} a_{33} + a_{23} a_{31}}{a_{11} a_{22} a_{33} – a_{11} a_{23} a_{32} – a_{12} a_{21} a_{33} + a_{12} a_{23} a_{31} + a_{13} a_{21} a_{32} – a_{13} a_{22} a_{31}}&\frac{a_{11} a_{33} – a_{13} a_{31}}{a_{11} a_{22} a_{33} – a_{11} a_{23} a_{32} – a_{12} a_{21} a_{33} + a_{12} a_{23} a_{31} + a_{13} a_{21} a_{32} – a_{13} a_{22} a_{31}}&\frac{- a_{11} a_{23} + a_{13} a_{21}}{a_{11} a_{22} a_{33} – a_{11} a_{23} a_{32} – a_{12} a_{21} a_{33} + a_{12} a_{23} a_{31} + a_{13} a_{21} a_{32} – a_{13} a_{22} a_{31}}\\\frac{a_{21} a_{32} – a_{22} a_{31}}{a_{11} a_{22} a_{33} – a_{11} a_{23} a_{32} – a_{12} a_{21} a_{33} + a_{12} a_{23} a_{31} + a_{13} a_{21} a_{32} – a_{13} a_{22} a_{31}}&\frac{- a_{11} a_{32} + a_{12} a_{31}}{a_{11} a_{22} a_{33} – a_{11} a_{23} a_{32} – a_{12} a_{21} a_{33} + a_{12} a_{23} a_{31} + a_{13} a_{21} a_{32} – a_{13} a_{22} a_{31}}&\frac{a_{11} a_{22} – a_{12} a_{21}}{a_{11} a_{22} a_{33} – a_{11} a_{23} a_{32} – a_{12} a_{21} a_{33} + a_{12} a_{23} a_{31} + a_{13} a_{21} a_{32} – a_{13} a_{22} a_{31}}\end{array}\right]\end{aligned} \]
それぞれ、サラスの方法、クラメルの公式と呼ばれるものです。この方法での手計算は大変なので、コンピュータで行うのが楽ですね。
参考:サラスの方法、n次行列式の覚え方:置換を使った定義を理解する
最後に、ヴァンデルモンドの行列式
\[ \begin{aligned}V_n :=\det \begin{pmatrix} 1& \dots &1\\ \lambda_1 &\cdots&\lambda_n \\ \vdots & \vdots \\ \lambda_1^{n-1}&\cdots& \lambda_n^{n-1}\end{pmatrix}\end{aligned} \]
を求めてみましょう。これは差積
\[ \begin{aligned}V_n = \prod_{1 \leq i <j \leq n} (\lambda_j -\lambda_i)\end{aligned} \]
と等しいことが知られています。
行列の中身としては、j列目に変数\(\lambda_j\)があり、i行目にそのべき乗\(\lambda_j ^{i-1}\)が並んだものです。
1 2 3 | V3 = [Sym("λ$j")^(i-1) for i=1:3 for j=1:3].reshape(3,3) V3.det() factor(V3.det()) |
\[ \begin{aligned}\left[ \begin{array}{rrr}1&1&1\\λ_{1}&λ_{2}&λ_{3}\\λ_{1}^{2}&λ_{2}^{2}&λ_{3}^{2}\end{array}\right]\end{aligned} \]
\[ \begin{aligned}- λ_{1}^{2} λ_{2} + λ_{1}^{2} λ_{3} + λ_{1} λ_{2}^{2} – λ_{1} λ_{3}^{2} – λ_{2}^{2} λ_{3} + λ_{2} λ_{3}^{2}\end{aligned} \]
\[ \begin{aligned}- \left(λ_{1} – λ_{2}\right) \left(λ_{1} – λ_{3}\right) \left(λ_{2} – λ_{3}\right)\end{aligned} \]
因数分解すると、確かに差積であることがわかりますね。
5次にして計算してみます。
1 2 | V5 = [Sym("λ$j")^(i-1) for i=1:5 for j=1:5].reshape(5,5) @time factor(V5.det()) |
\[ \begin{aligned}\left[ \begin{array}{rrrrr}1&1&1&1&1\\λ_{1}&λ_{2}&λ_{3}&λ_{4}&λ_{5}\\λ_{1}^{2}&λ_{2}^{2}&λ_{3}^{2}&λ_{4}^{2}&λ_{5}^{2}\\λ_{1}^{3}&λ_{2}^{3}&λ_{3}^{3}&λ_{4}^{3}&λ_{5}^{3}\\λ_{1}^{4}&λ_{2}^{4}&λ_{3}^{4}&λ_{4}^{4}&λ_{5}^{4}\end{array}\right]\end{aligned} \]
\[ \begin{aligned}\left(λ_{1} – λ_{2}\right) \left(λ_{1} – λ_{3}\right) \left(λ_{1} – λ_{4}\right) \left(λ_{1} – λ_{5}\right) \left(λ_{2} – λ_{3}\right)\\ \left(λ_{2} – λ_{4}\right) \left(λ_{2} – λ_{5}\right) \left(λ_{3} – λ_{4}\right) \left(λ_{3} – λ_{5}\right) \left(λ_{4} – λ_{5}\right)\end{aligned} \]
約2秒とそれなりの時間がかかり、工夫をしないとこれ以上のサイズは大変そうです。
以上、Julia(SymPy)で行列式、逆行列を求める方法を紹介してきました。
数値計算的には、線形方程式\(Ax=b\)を解くためには、逆行列の掛け算\(x= A^{-1}b\)は誤差が増えるので、LU分解などの方法が推奨されています。
しかし、記号を使って理論的に行列式や逆行列を求めたいときに、すべて手計算するのも大変です。SymPyならある程度は求められると知っておくのと良いでしょう。
木村すらいむ(@kimu3_slime)でした。ではでは。
コロナ社 (2020-03-26T00:00:01Z)
¥7,353 (コレクター商品)
こちらもおすすめ
サラスの方法、n次行列式の覚え方:置換を使った定義を理解する
Julia(SymPy)で行列の和・積・転置、線形方程式を解く方法