LoginSignup
2
2

More than 5 years have passed since last update.

Julia+SymPyで行列式が正しく求められなかった件

Last updated at Posted at 2017-11-10

Julia+SymPyで行列式が正しく求められない場合???

JuliaのSymPyパッケージで文字の入った行列から行列式を求めようとしていました。うまく行く場合とうまく行かない場合があり、原因ははっきりとさせることできませんでしたが、調べた結果を記します。

実行環境

Windows 7 Professional Service Pack 1
Julia 0.6.1
SymPy 0.5.4
Conda 0.7.0

例1【Julia版】(エラー出たコード!)

Julia
using PyCall
using SymPy

@syms L a

M=[1        0        0 1
   1        0        0 0
   cos(L*a) sin(L*a) L 1
   cos(L*a) sin(L*a) 0 0]

det(M)       #←なぜかエラー出る
SymPy.det(M) #←これもエラー出る
Base.det(M)  #←これもエラー出る

例1'【Python版】(例1と同じ内容のつもりです)

本家Pythonで行うとやはりエラーなど出ませんでした。

Python
#以下3行はPythonでmatplotlibを使うためのおまじないみたい。よく意味は理解していない・・・
%matplotlib inline
from sympy import *
init_printing()

L,a = symbols("L a")

M=Matrix([[1,       0,       0,1],
          [1,       0,       0,0],
          [cos(L*a),sin(L*a),L,1],
          [cos(L*a),sin(L*a),0,0]])

det(M)

※検証のため使ったことのないPythonをなんとか操作しました・・・

例2【Julia版】(前の例に続けて)

下記のような行列Nでは正しく計算できました。何でだろう???

Julia
N=[1           0          0  1
   0           a          1  0
   cos(L*a)    sin(L*a)   L  1
   -a*sin(L*a) a*cos(L*a) 1  0]
det(N) |> simplify #←エラー出ない

例3【Julia版】(前の例に続けて)

仕方が無いので、Julia上では行列式を求める関数を自作して、それを使うようにしてみました。エラーは出ません。

Julia
#■小行列を求める関数
function submat(M,R::Int,C::Int)
  (Rend,Cend)=size(M)
  R>Rend  && error("行番号の指定が行列の最大行数を超えています")
  C>Cend  && error("列番号の指定が行列の最大列数を超えています")
  R==1    ? M=M[2:end,  :] :
  R==Rend ? M=M[1:end-1,:] :
            M=[M[1:R-1, :];M[R+1:end, :]]
  C==1    ? M=M[:,2:end  ] :
  C==Cend ? M=M[:,1:end-1] :
            M=[M[:,1:C-1] M[:,C+1:end]]
  (1,1)==size(M) ? M[1] :M
end
#↑の自作関数はSymPyのminorMatrix(M,R-1,C-1)と同じような働きになっているはず

#■行列式を求める関数
function determ(M)
  ndims(M)>=3 ? error("行列になっていないので計算できません") : begin
    s=size(M)
    s==()      && return M    #Mが非Arrayだった場合
    s==(1,)    && return M[1] #Mが1-element Arrayだった場合
    s==(1,1)   && return M[1] #Mが1×1 Arrayだった場合
    s==(2,2)   && return M[1,1]*M[2,2]-        #Mが2×2 Arrayだった場合
                         M[1,2]*M[2,1]
    s==(3,3)   && return M[1,1]*M[2,2]*M[3,3]+ #Mが2×2 Arrayだった場合
                         M[1,2]*M[2,3]*M[3,1]+
                         M[1,3]*M[2,1]*M[3,2]-
                         M[1,1]*M[2,3]*M[3,2]-
                         M[1,2]*M[2,1]*M[3,3]-
                         M[1,3]*M[2,2]*M[3,1]
    #↑s==(2,2)とs==(3,3)の場合はたすき掛けの法で求めるようにしている
    # この場合を設けなくても関数は成立するが、計算時間が早そうなので入れている
    s[1]!=s[2] && error("正方行列ではないので計算できません")
    return sum((-1)^(1+i)*M[1,i]*determ(submat(M,1,i)) for i in 1:s[2])
  end
end

determ(M)             #エラー出なくなった
determ(N) |> simplify #これもエラー出ない。例2と同じ結果

最後に

長いことはまっていました!
自分用のメモのつもりで書きましたが、
原因分かる方いましたら教えて欲しいです・・・

2
2
5

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
2
2