の中で固有値を求める関数eivals()を使うとエラーが発生する場合がある、という話を書きました。具体的には以下のような行列mで、
(%i1) m:matrix([3,0,0,0,-1,0,0,-1,0,0,-1,0,0],[ 0,1,-1,0,0,0,0,0,0,0,0,0,0],[ 0,-1,1,0,0,0,0,0,0,0,0,0,0],[ 0,0,0,0,0,0,0,0,0,0,0,0,0],[ -1,0,0,0,2,0,0,0,-1,0,0,0,0],[ 0,0,0,0,0,1,0,0,0,0,-1,0,0],[ 0,0,0,0,0,0,0,0,0,0,0,0,0],[ -1,0,0,0,0,0,0,2,0,0,-1,0,0],[ 0,0,0,0,-1,0,0,0,1,0,0,0,0],[ 0,0,0,0,0,0,0,0,0,1,-1,0,0],[ -1,0,0,0,0,-1,0,-1,0,-1,4,0,0],[ 0,0,0,0,0,0,0,0,0,0,0,0,0],[ 0,0,0,0,0,0,0,0,0,0,0,0,0 ]);
$$ \tag{%o1} \begin{pmatrix}3&0&0&0&-1&0&0&-1&0&0&-1&0&0\\ 0&1&-1&0&0&0&0&0&0&0&0&0&0\\ 0&-1&1&0&0&0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0&0&0&0&0&0\\ -1&0&0&0&2&0&0&0&-1&0&0&0&0\\ 0&0&0&0&0&1&0&0&0&0&-1&0&0\\ 0&0&0&0&0&0&0&0&0&0&0&0&0\\ -1&0&0&0&0&0&0&2&0&0&-1&0&0\\ 0&0&0&0&-1&0&0&0&1&0&0&0&0\\ 0&0&0&0&0&0&0&0&0&1&-1&0&0\\ -1&0&0&0&0&-1&0&-1&0&-1&4&0&0\\ 0&0&0&0&0&0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0&0&0&0&0&0 \end{pmatrix} $$
固有値を求めると、
(%i2) eivals(m);
part: fell off the end.
#0: eigenvalues(mat=matrix([3,0,0,0,-1,0,0,-1,0,0,-1,0,0],[0,1,-1,0,0,0,0,0,0,0,0,0,0],[0,-1,1,0,0,0,0,0,0,0,0,0,0],[0,0...)(eigen.mac line 94)
#1: eivals(mat=matrix([3,0,0,0,-1,0,0,-1,0,0,-1,0,0],[0,1,-1,0,0,0,0,0,0,0,0,0,0],[0,-1,1,0,0,0,0,0,0,0,0,0,0],[0,0...)(eigen.mac line 189)
#2: eivals(_l=[matrix([3,0,0,0,-1,0,0,-1,0,0,-1,0,0],[0,1,-1,0,0,0,0,0,0,0,0,0,0],[0,-1,1,0,0,0,0,0,0,0,0,0,0],[0,...)
-- an error. To debug this try: debugmode(true);
このようにエラーが発生します。これは明らかにmaximaのバグです。なんでこんなエラーが発生するのでしょうか。
細かい理由を調べる前に、まずこの行列の特性方程式を求めて、そいつを因数分解してみます。
(%i3) factor(charpoly(m,x));
$$ \tag{%o3} -\left(x-2\right)\,\left(x-1\right)\,x^6\,\left(x^5-13\,x^4+60\,x^3-118\,x^2+93\,x-21\right) $$
このように因数分解出来ました。solve()を使って解くと以下のように表示されます。
(%i4) solve(%,x);
$$ \tag{%o4} \left[ x=0 , x=1 , x=2 , 0=x^5-13\,x^4+60\,x^3-118\,x^2+93\,x-21 \right] $$
5次の因数の部分が解けなくて、その部分だけ方程式がそのまま解のリストに入っています。
これは、、、いかにもダメっぽそうです。
eigen.macのコードを調べるとこの部分の処理が間違っていることがわかり、そこを修正しました。修正版のeigen.macを読み込みます。
(%i5) load("eigen.mac")$
修正版のeivals()を使って固有値を求めてみます。
(%i6) eivals(m);
(*) eigenvalues: solve is unable to find some of the roots of the characteristic polynomial.
$$ \tag{%o6} \left[ \left[ 2 , 0 , 1 \right] , \left[ 1 , 6 , 1 \right] \right] $$
solve()の解として求まっている2, 0, 1についてはその重複度まで正しく返しています。また5次の因数は解けずそこに対応する固有値は求まっていないので、そのことをメッセージとして表示しています。
この修正コードはMaxima本家のgitレポジトリにコミットしましたので、5.36.1以降のリリース(おそらく5.37.0)に含まれることになります。
先日Google PlayにリリースしたMaxima on Android 2.8のMaximaは5.36.1をインテグレーションしてありますが、上記の修正は取り込んで対応済みです。
