02/07/14 16:52
>>59
11章くらいです。
今ちょっと H:matrix([e1,v,0],[v,e2-e1,v],[0,v,-e2]);
これの固有ベクトルを求めることが sci.math.symbolicmathのnewsgroupで
話題になっていて mapleでは 求まって maximaとmupadでは求まらない、
とあったので調べていた。原因は 固有値までは求まるのだが 固有方程式が
解けない。ただの3元1次方程式が 解けないのだから 厨房にまけ?
なお固有値の計算をするとき load("eigen")$としてライブラリをよびだす
のは eigen.mc (新しいバージョンではeigen.mac)をLISPにコンパイルした
ものは 完全に同じ動作をしないからである。このへんは コモンLISPのむずな
ところで 要するに eigen.macをきれいに書けばよいのだが VER7.0くらいまで
は見込がなさそうなので いちいち始めに一回loadする。eigenのライブラリは
ちょっとプログラムスタイルが 悪いというかコワイ部分(ループのなかでルー
プ終端値を書き換えたりする、HP-25のP.ヘンリッチを思い出す)がある。
時代もののせいだろうか。
L1:eigenvalues(H);
L1[1][1];は 固有値のひとつ 。たった３次の正方行列の固有値のひとつが
このようになること自体に 数値処理と数式処理の深いみぞが見れる。
(-SQRT(3)*%I/2-1/2)*(SQRT(-32*v^6-(21*e2^2+6*e1*e2+21*e1^2)*v^4
-(24*e2^4+6*e1*e2^3-36*e1^2*e2^2 +6*e1^3*e2+24*e1^4)
*v^2-4*e2^6+12*e1*e2^5+3*e1^2*e2^4 -26*e1^3*e2^3+3*e1^4*e2^2+12*e1^5*e2
-4*e1^6)/(6*SQRT(3)) +(e2*v^2+e1*(-v^2-e2^2)+e1^2*e2)/2) ^(1/3)
+(SQRT(3)*%I/2-1/2)*(2*v^2+e2^2-e1*e2+e1^2)
/(3*(SQRT(-32*v^6-(21*e2^2+6*e1*e2+21*e1^2)*v^4
-(24*e2^4+6*e1*e2^3-36*e1^2*e2^2+6*e1^3*e2+24*e1^4)
*v^2-4*e2^6+12*e1*e2^5+3*e1^2*e2^4-26*e1^3*e2^3
+3*e1^4*e2^2+12*e1^5*e2-4*e1^6) /(6*SQRT(3))
+(e2*v^2+e1*(-v^2-e2^2)+e1^2*e2)/2)^(1/3))
このくらいになると 検算をしておくのが普通です。(むやみに 数式処理シス
テムを信じないこと)
もとの特性方程式はcharpoly(H,x)；としてもでてくるが,determinant(H
-x*ident(3));としても同じ。ratsimp(%)；で整理する。
-x^3+(2*v^2+e2^2-e1*e2+e1^2)*x+(e2-e1)*v^2-e1*e2^2+e1^2*e2
という３次方程式(=0は省略される)の解になっていることを確かめる。
直接代入すると subst([x=L1[1][1]],この方程式)；ratsimp(%);で
QUOTIENT by ZERO (０割のえらーでっせ)
-- an error. Quitting. To debug this try DEBUGMODE(TRUE);)
えーーそんなあほな この程度の検算もできん数式処理システム降りさせても
らいますとなるところ。でもちょっとまって固有値L1[1][1]をよく見てみると
６行目の/は最後までかかっていて そのなかに別の/がある。つまり数値用と
兼用のライブラリだから この固有値は 整理されていない。計算の途中で0で
割っている用に見えてもみかけだけやないだろうか?ということで
L2:ratsimp(L1[1][1]);subst([x=L2],この方程式);rasimp(%);であたりまえ
だが0が 返って来て 検算ができた。(続く)