結局使わなかったけどメモ
背景
http://jikosoft.com/cae/engineering/strmat05.html
みたいな感じ。
ただ胴体が曲がってる場合（とか）を想定して、
応力じゃなくて慣性モーメントの主軸を求めようかなと。
やることは全く同じのようだ。
結論メモ
A をもとの行列とする。
R^{-1} A R = A' とすると、
対角行列 A' の対角項 = 固有値 = 主軸座標系での値
R が固有ベクトル（固有列行列）をまとめた行列 = 座標変換行列（でいいはず）、なので R^{-1} は R^T でOK. まぁ確かめは AR = RA' になってるかでいいので R^{-1} は要らないけど。
んでグリッドを座標変換するときは、グリッド位置を p とすると、
p' = R p
でいいはず。その前に止めたからやってないけど。
というのは、そもそも非対角成分の値が二桁も小さかったので、
A' が A とほとんど変わらなかった。
無視して全然OKっていうか無視すべき。
A っていうか慣性モーメントなので J という記号使ってるけど。
手法
固有値問題なんだけど、一般の n x n のメイトリックスでなくて 3 x 3 なので、
「直説法の方が誤差小さいんじゃね（妄想）」という方針だった。
でも本とかウェブには n x n に対応できるような Jacobi method とか QR method とかの反復法しか出てない。
なので手で解こうとした。
要は det( A -λI) = 0 を解きゃいいはず。
3x3 だから三次方程式。
http://en.wikipedia.org/wiki/Cubic_equation
この解の公式そのままでもできるかもだけど、なんかうまく行かなかった（どうも公式の問題じゃなくて三次方程式の展開をミスってたっぽいorz）
ともあれ、同じこと考えてやってる人がいた（Cで）
http://d.hatena.ne.jp/malibu-bulldog/20081009/1223516459
http://ja.wikipedia.org/wiki/%E4%B8%89%E6%AC%A1%E6%96%B9%E7%A8%8B%E5%BC%8F
http://www004.upp.so-net.ne.jp/s_honma/urawaza/eigenvector.htm
これを Fortran に書き直したらうまくいったっぽいけど、
要素 (1,2) と (2,1) が完全に 0 だとおかしくなる（ちゃんと読んでないけどゼロ割っぽい？
1.0d-5 とか微小な非ゼロいれとけばいける。
あと倍精度より4倍精度の方がよかった。
related articles
http://www22.atwiki.jp/linearalgebra/
あと『プログラミングのための線形代数』が完全にこの問題のためにあるんじゃね的なシンクロ率だった。いまさらだけど。しかも使わないけど。
感想
「n が大きい場合には大変なので～」ってよくあるけど「n > 2 の場合には」じゃね？
No comments:
Post a Comment