天文計算と計算誤差

電卓

最近は、あまりコードを書かなくなった。他のことでやることが多いからだ。今書いているコードは殆どが、趣味の天文計算に関するものだ。Python、Octave、GNU Fortranで書いている。Webに関してはさっぱりである。まあ、この歳でエンジニアでもなかろう。

Octaveに関して言えば、昔のコンピュータを思い出す。WindowsもCP/MもMS-DOSもなかった頃だ。スイッチオンでBASICが使え、コードを書いて実行するだけでディスプレイに結果が表示される。Octaveは、端末でも実行できるが、GUIが用意されているので、そちらを利用している。

Octaveは言語仕様から言えば、プログラミングしやすい言語である。基本の演算桁数は5桁だが、format longで16桁の演算桁数になる。本来は科学技術計算向けである。GNUPlotのようにグラフも表示できる。GNU Fortranでコードを書く前に、プリプロセッサ的にコードを書いて試すのが良いかもしれない。

天文計算に限らないが、コードを書いていると、どうしても計算誤差がともなう。計算誤差は、丸め誤差と桁落ちの2つである。実数つまり浮動小数点は仮数と基数で表されるのはご存知であろう。たとえば、0.03141592653589という数値は、基数を10とした場合、3.141592653589×10-2で表される。仮数は一定の桁数の数値として表されるために、前述の丸め誤差と桁落ちが起きる。

3.141592653589×10-2の有効桁数は14桁数である。これを8桁が有効な変数に代入した場合、0.03141593となったりする。この誤差を丸め誤差呼ぶ。とくに丸め誤差が蓄積されるような演算でない限り、必要桁数で切り捨てもしくは四捨五入、切り上げなどの処理を行えばよい。ただ、コンピュータ内部で扱う数値は基数が10ではなく2の場合は、問題はややこしくなる。2進数で扱う浮動小数点は循環小数となり、たとえば、数値が7桁あり、変数も有効桁数が7桁だとしても丸め誤差が入るケースがある。

丸め誤差の影響を軽減するには、単精度を倍精度で扱うといった仮数部の桁数を多くとる。ただし、単精度の数値を倍精度の変数に代入しないように注意する。GNU Fortranでは、倍精度数値とするには、数値の末尾にd0を付加する。

桁落ちは、桁数が近い浮動小数点の減算を行うと、有効桁数が落ちることをいう。軽減するには、丸め誤差と同じように、倍精度として扱う。もしくは、計算のやり方を変えることである。計算誤差は、大きさ数値と小さい数値の加算、減算を行った際に、結果が反映されない情報落ちがあり、数値演算によって厄介なものとなっている。実際は、そう神経質になることはないが、注意するに越したことはない。

天文計算をやっていると、精度にこだわるようになってしまうことがあり、なかなか抜け出せない。実際には演算桁数は14桁もあればよく、必要桁数は10桁でだいたいは足りる。ただし、精密な位置計算ではもう少し必要桁数が多くなる。とはいえ、現在では最大演算桁数が16〜17桁なので、それ以上の演算では10進数演算が必要になるが、そこまでして精度を上げても、位置計算や軌道計算において、あまり意味があるとは思えない。

求める精度についても、計算対象によって分けるべきで、日の出日の入りに秒まで求めたとしても、それは計算結果であって、あまり意味がない。これは確かめる術がないためだ。もし、地球に大気がなく、地平線に何もなければ、計算結果が正しいかどうか観測できるかもしれないが、実際は大気の影響や気象によって、正確に観測はできないと思ってよい。

計算結果を自分で確認するのであれば、概略計算でも十分なのである。小中口径の望遠鏡で低倍率であれば、実視界(見える範囲を角度で表したもの)を広く取ることができるので、計算結果が天文雑誌に掲載されているような位置と多少異なったとしても、視界内に捉えることができるからだ。もちろん、大口径であれば、いくら低倍率とはいえ実視界が狭いために、概略計算は致命的であるといってよいほどで、精度の高い計算を行う必要がある。

いずれにしても、プロでない限り、精度にこだわるとストレスしかない。天文計算では、必要な精度を得られるかどうかが大切なのであり、それ以上の無駄な精度を求めても、それは自己満足なのである。