2008年4月1日火曜日

誤差の問題

ここ数日ほんとに困ってたある特定の条件に対して数値計算がヘンな値を返すよくわからん問題が解決したぽい。
たぶん、0が正しい値である箇所に1.0E-17とかの超小さいが有限な値が入っていた+対角化に使っていたライブラリのルーチンが適切でなかったのが問題なようだ。すべては数値計算誤差の問題っぽい。物理の考え方が間違っていたわけではなさそうだ。
コンピューターは指示された通りに馬鹿正直に動く。プログラムの意図とか数学を理解しているわけではない。数学的には0となる解もコンピューターの内部的には「メモリのある区画までは0」みたいに扱われるわけで、どうしたってゴミがついてまわる。下手すると(0-0)×1000000で有限の値とかが返ってくる。ってなことは数値解析かなんかで腐るほどきいたけど、まさか自分がそんなんと向き合う必要があるとはねぇ。誤差とか精度についてもういっぺん考える必要があるか。

対角化の数値計算について。行列の対角化なんかは大抵Lapackライブラリを使う。自分でつくるよか精度・速度ともにぶっちぎりでよい。対角化をサポートしている代表的なルーチンはdgeevとdsyevとか。dgeev/zgeevは実数/複素数要素の正方行列ならなんでもOK。dsyev/zheevは実対称/エルミート行列の対角化を行うことができる。量子力学で出てくるハミルトニアンはエルミートなので後者を使ったほうがいいようだ。今の問題ではすべて実数要素なので、これまでdgeevを使っていたが、これだとエラーが出る条件もdsyevを使うようにするとうまく行った。対称性を使ってるので安定性がいいのか?よくわからん。
固有値問題にはいくつか解法があるようだ。ヤコビ法、QR法、ハウスホルダー法、ランチョス法あたりが聞いたことある。名前しか知らん。おもしろそうではあるけど数学にがてな私が近づくとやばそうな領域。固有値ヤバイ。それぞれの手法は対称or非対称で得手不得手があり、丸め誤差に対する頑健さもちがうらしい。2つのルーチンはおそらく違う手法で実装されててその辺が効いてきてるんだろう。まぁとりあえずここのサイトが凄かった。プロの仕業

1 コメント:

匿名 さんのコメント...
このコメントはブログの管理者によって削除されました。