Python
sympy
Julia

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

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と同じ結果

最後に

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