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と同じ結果
最後に
長いことはまっていました!
自分用のメモのつもりで書きましたが、
原因分かる方いましたら教えて欲しいです・・・