2013/03/27

彗星の尾の計算(ダストテイル)その2

更新がだいぶ間が空いてしまいました・・・

今、C/2011 L4が見えてますが、イマイチ『ショボい!』なんて言われてますが、
STEREO-B衛星が撮影したC/2011L4です。
http://spaceweather.com/gallery/indiv_upload.php?upload_id=78033
シンクロニックバンドがハッキリ見えて綺麗なダストの尾です。
彗星物理屋さんはワクワクしているでしょう。

さて、ダストテイルの計算ですが・・・
運動方程式が解けません!


・ダストの運動方程式

更新してない間、『遊んでない時は真面目に考えてました』という証明のために、まとめ。

太陽系内の天体は、自ら運動しない限り、惑星から塵や埃まですべて

Comet1501

この運動方程式に従っているらしい。
GMは太陽のガウスの引力定数で、今まで使っていた k とは
Comet1504 
という事になる。

さて、彗星核から飛び出た(英語では、release と書かれている)ダストは
光圧と引力の両方の力を受ける。
その力の比を βという係数で表している。

Comet1404

このβからμを計算し、先の運動方程式に入れると

Comet1502_2
Comet1503

これがダストの運動方程式。
μによりGMの値は

μ=1.0 (β=0)  GMのまま=彗星核と同じ運動
μ=.0 (β=1)  GMが0になる=ダストに力は掛からない=等速直線運動
μ=-1.0 (β=2) GMが負になって太陽と反対方向に運動する。

こんな感じに変化する。


・ダストの位置計算

先のダストの運動方程式を解いて、軌道要素を求めようとしたが、これが無理!
よくよく考えたら、軌道要素は「太陽を焦点とする二次曲線運動」が前提、
β≧1.0の時、太陽引力とは無関係な運動になってしまう。よって、

ダストの軌道要素が求まるのはβ<1.0

という事になる。
この結論を出すのに、こんなにも時間が必要だった。
ダストの軌道要素を求めたい人は、

A THEORY OF DUST COMETS. I. MODEL AND EQUATIONS
The application of Finson-Probstein dynamical theory to the dust tail of P/Halley

が参考になる。

軌道計算を求めるのを諦めれば、残る方法は一つ。
「数値解析(数値計算)」で求める方法です。


・数値解析とは

位置推算から離れている気がしますが・・・まぁ、この際だから行くところまで・・・

数値解析とは、数値計算ともいいます。今後は「数値計算」と書きます。
微分式を解かず(式にせず)数値を計算して解を求める方法です。
簡単な例で

Comet1505

これは、自由落下の速度を求める式で、これも立派な微分式です。
この速度(gt)から、ある時間(t)の x 位置を求めるのが微分方程式を解くという事。
微積分を習った人なら、この式は簡単に解ける。時間(t)で積分すればよい。

Comet1506

数値計算の場合は、

Comet1507

(N=刻み幅)
微分式に直接数値を代入して総和を求めます。

実際にやってみましょう。
初期位置の x0 = 0 として、10秒後の x の位置を求めます。
刻み幅(N)を2秒間隔とすれば、t0 = 0, t1 = 2, t2 = 4, t3 = 6, t4 = 8になります。

Comet1508

あまりに単純で原始的で・・・笑われそうですが、
微積分を知らなくても計算できるし、微分式を悩みながら解く必要もない。

でも・・・数値計算には欠点がある
この正しい答え(厳密解という)は、積分で出した式の方で

Comet1509

数値計算で出した 392 と全然違う!
この違いは、刻み幅の2秒が少々粗かったと思うでしょ?
では、刻み幅を変えて計算してみると
刻み幅が1秒なら、441
    0.1→494.9
    0.01→490.49
    0.0001→490.0049
    0.0000001→490.000005
確かに厳密解に近づいてますが・・・490.0ドンピシャにならない。

数値計算法(数値解析法)で出した答えは、必ず誤差がある
誤差は数値計算の宿命という事です。
この誤差を小さくするために色々な方法があります。
先に示した計算例は、「オイラー法」で、誤差が大きいが単純なやりかた。
天文分野では、「4次のルンゲークッタ法」がちょっとした計算に使われています。

次回は、4次のルンゲ-クッタ法でダストの位置を計算してみます。