拝読中の溝口純敏様著「Maxima を使った物理数学基礎演習ノート」（以下「演習ノート」と略）の今回は「対角化」の練習です。練習といってもひたすら他力本願、Maximaさまに御すがりするバカリ。さて繰り返す対角化を関数化、そこにまだ問題ありーの。中途半端なところだのう。タイパなんか考えるなよ。なんまんだぶ。
※ 「忘却の微分方程式」投稿順 index はこちら
※ MaximaおよびそのGUIであるwxMaximaの以下バージョンを使用させていただいております。
-
- wxMaxima 22.04.0
- Maxima 5.46.0(x86_64-w64-mingw32)
- SBCL 2.2.2 (SBCL = Steel Bank Common Lisp )
※Maxima を使った物理数学基礎演習ノート は以下のバージョンをダウンロードさせていただきました。
令和4 年3 月 第八回改訂
今回「演習ノート」を学ばせていただくのは、「4.2 行列とテンソル」の「4.2.10 対称行列の対角化」の後半あたりデス。
行列対角化用の関数を定義
さて演習ノートの「4.2.10 対称行列の対角化」の後半ではいくつか具体的な数値が詰まった行列を持ってきて対角化をやってます。でもね、そういえば過去回でも以下の回で同じような練習やってます。
そのとき、「定型的だし、関数化できそうだなあ」と感想を述べつつ、関数化はしてませんでした。そこで今回、心を入れ替えて（ホントか？）関数化してみました。しかし中途半端な関数でしかないのでありますが。中途半端なヤツが以下に。
/* diagonalize2 */
/* ２ｘ２行列 A の対角化 */
/* 返却：対角化されたA */
diagonalize2(A):=block(
[EA, X1, X2, P, PI],
EA: eigenvectors(A),
X1: transpose(EA[2][1][1]),
X2: transpose(EA[2][2][1]),
P: addcol(4*X1, X2),
PI: invert(P),
(P^^-1) . A . P
)$
/* diagonalize3 */
/* ３ｘ３行列 B の対角化 */
/* 返却：対角化されたB */
diagonalize3(B):=block(
[EB, X1, X2, X3, P, PI],
EB: eigenvectors(B),
X1: transpose(EB[2][1][1]),
X2: transpose(EB[2][2][1]),
X3: transpose(EB[2][3][1]),
P: addcol(X3, X2, 3*X1),
PI: invert(P),
(P^^-1) . B . P
)$
/* diagonalize3e */
/* ３ｘ３行列 A の対角化 演習ノート式 */
/* 返却：対角化されたA */
diagonalize3e(A):=block(
[n, B, C, D, B1, B2, B3, Q, m, P, P1],
eigenvectors(A),
n : length(A),
B : eigenvectors(A),
C : zeromatrix(n, n),
B1:transpose(B[2][1][1]),
B2:transpose(B[2][2][1]),
B3:transpose(B[2][2][2]),
Q:addcol(B1,B2,B3),
m : length(Q),
B : transpose(Q),
n : length(B),
C : gramschmidt(B),
D : zeromatrix(n, m),
for i:1 thru n do
(D[i] : unitvector(C[i])),
P:transpose(D),
P1:ratsimp(invert(P)),
ratsimp(P1.A.P)
)$
最初の２つ、
-
- diagonalize2
- diagonalize3
は、上記の過去回で練習につかったシーケンスをそのまま関数化したものでやんす。行列サイズを判断して切り替えるようなことはまったくしていないので、「２」の方が２ｘ２行列用、「３」の方が３ｘ３行列用です。一方、
-
- diagonalize3e
は「演習ノート」のP.97の左側に掲載されている例を関数化してみたものです。この例は「ちょいとクセあり」で上記の単純なdiagonalize3関数では失敗してしまうものためです。
「中途半端」の理由の一つは、渡される行列が対角化できるかどうかなどまったくチェックせずに、ただただ決められたシーケンスで処理しているところです。よって、上記の関数どもを呼び出す前に、対角化できる行列なのか否かは自分でチェックしてからにする必要があります。まあ、今回みたいに対角化の練習問題なら対角化できるに決まってますけど。
まずは２ｘ２行列
上記の関数定義をバッチ実行すれば、関数どもが使えるようになります。まずは過去回の２ｘ２行列を対角化する例。
余勢を駆って、演習ノートP.96の右側の例。
対角化できてるみたい。
お次は３ｘ３行列なんだが
関数切り替えて、３ｘ３行列を処理してみます。まずは過去回でもやった例。
この例はOKね。
しかし、「演習ノート」のP.97の左側に掲載されている例ではdiagonalize3はエラーを吐きます。こんな感じ。
なんでかと eigenvectors(A) を求めてみると、重解っす。diagonalize3は重解を想定していないのでエラーになってます。
そこで「演習ノート」のP.97左側のシーケンスをそのまま関数化した diagonalize3eを使ってみると、上記のように対角化できているみたい。
この辺、勝手に判定して適切な処理を行ってくれるといい感じなんだが、上記は中途半端。まあ、お惚け老人が介護の合間にやっているので許してちょ。