[Exp2010]実際に計算してみる (３体問題)

スケジュール表・各回資料 (07/09)

本日の実習でも前回と同様に, すでに途中まで出来上がっているプログラムに手を加えて作業を行います. 前回と同じ要領で情報実験機に以下のファイルをダウンロードしてください.

• 3body.f90
• 3body_sol.f90

ちなみに, 今回用いるプログラムには全て 2 次のルンゲクッタ法が使われています.

実習 その 2-1 : プログラムの穴埋め

ダウンロードした 2 つのファイルのうち, 3body.f90 を開いてみましょう. 中には Fortran のプログラムが書かれた文がならんでいます.

しかし前回と同様に, プログラムが正しく書かれていないところがあります. 前回と同じ要領で穴埋めしていきましょう.

今回のプログラム中には, 穴埋めが必要な箇所が全部で 10 行あります. 全ての箇所に正しい式を書き込み, コンパイル時にエラーが出ないようにしてください.

初期条件について

• 天体の質量

  3 体問題を扱うにあたり, 2 体問題では取り扱わなかった質量の設定も初期条件に含まれています. それはここにあたります.

  !---- 各惑星の質量を決める 
  
     factor = 1.0d0
     mp(1)  = 1.0d0
     mp(2)  = 9.5479d-4*factor
     mp(3)  = 2.859d-4*factor

  これは現在の太陽系での太陽・木星・土星の関係にあたります. 中心星, すなわち太陽の質量を 1 としているところに注意してください. なお factor を変えることで太陽系において木星と土星の質量がそのまま factor 倍された場合を計算できます.

• 惑星の初期位置

  惑星の初期位置は次のように与えられています.

  !---- 初期に惑星を次のように配置する. 
  
     xp(1)=    0.0d0 
     yp(1)=    0.0d0
  
     xp(2)=    1.0d0 
     yp(2)=    0.0d0
  
     xp(3)=    1.837d0
     yp(3)=    0.0d0

  これも現在の太陽系での位置関係にあたります. ここも後で変えて頂いてかまいませんが、変えるのはなるべく, xp(3) だけにしてください.

• 惑星の初期速度

  惑星の初期速度は次のように与えられています.

  !---- おいた惑星を yp の正の方向に初速を与えて放り出す
  
     vx(1) = 0.0d0
     vy(1) = 0.0d0
  
     vx(2) = 0.0d0
     vy(2) = sqrt(mp(1)+mp(2))/sqrt(Abs(xp(2)))
  
     vx(3) = 0.0d0
     vy(3) = sqrt(mp(1)+mp(3))/sqrt(Abs(xp(3)))

  これも現在の太陽系での場合にあたります. 物理的意味としては, 惑星を y 方向にだいたいケプラー速度で飛ばしたということになります. 先の初期位置の y 座標を書き換えるならば, ここも適当に変える必要があります. 興味があればどうぞやってみてください.

実習 その 2-2 : 太陽系は安定か

ファイル 3body.f90 はデフォルトでは現在の太陽系での場合を表しています. これを用いて, 太陽系が本当に安定かどうかを考えてみましょう. その際, プログラム中のパラメータを次のようにしてください.

maxrevo = 5.0d4

上記のようにパラメータを設定して計算することは, 木星が 50000 回公転する期間を計算することに相当します. maxrevo は今回のプログラム内で, 冒頭の変数を定義した直後に登場しています. 各自で確認しておきましょう.

プログラムが完成したら前回と同様に,

$ gfortran 3body.f90
$ ./a.out

と入力してコンパイル行い, 計算を実行してください.

計算を実行に移すと, 端末の 1 行目には m(2) と m(3) の質量, 2 行目には factor, 3 行目以降にそれぞれの天体の初期位置と初速度が表示され, その後は計算回数が 5000 ステップごとに出力されていきます. また, "fort.14" というファイルが新しく作成され, 5000 ステップごとに計算結果が書き込まれていきます.

計算が終了したら, さっそく Gnuplot を起動して計算結果をグラフに出力してみましょう.

$ gnuplot
gnuplot>plot "fort.14" u 2:3

今回の計算で "fort.14"には,

• 1 列目に時刻
• 2, 3 列目に m(1) の x 座標と y 座標
• 4, 5 列目に m(2) の x 座標と y 座標
• 6, 7 列目に m(3) の x 座標と y 座標

という順で計算結果が出力されています. よって gnuplot で上のように入力すると, m(1) の軌道がプロットされるはずです.これと同様にして, m(2) と m(3) の計算結果も一緒にプロットしてみましょう.

ちなみに, ここで言う安定とは「計算している範囲で惑星の軌道が楕円形を保ち続ける」という意味だと思ってください.

※　どうしてもコンパイルがうまくいかない, もしくは軌道の図がうまく描けないときは, 3body_sol.f90 を使ってください.

初期条件の変更

• 惑星の質量を変化させるとどうなるでしょうか.
• 惑星の初期位置を変化させるとどうなるでしょうか.

初期条件を次のように変化させてみましょう.

• factor = 0.5, 1.0, 2.0, 4.0, 8.0, 16.0, 32.0, 64.0
• xp(3) = 1.1, 1.3, 1.5, 太陽系での値, 2.0, 3.0

惑星の軌道は惑星の質量や初期位置によってどのように変わりますか.

• これらの計算を全て行うのは, さすがに時間がかかってしまいますので, 各実験機で分担することにします.
• factor を 8 個用意したので, 各実験機につき factor 1 つ分をこなしてください.

レクチャー資料に登場したものとと同じ表を, 前のホワイトボードに用意してあります. 自分たちが担当している factor について, 6 通り全ての xp(3) を計算し終えた人から, 計算結果を記入しに来て下さい. 計算した結果, 全ての天体の軌道が安定であった場合は○を, 不安定なものがあった場合は×を表の中に記入してください.

さらに, この結果についてレポートにまとめ, 6/25(金) に出題した「課題 その1」と一緒に提出してください. 詳細はこちらです.

今回の実習が物足りない, と思われる方は次のテーマに取り組んでみましょう.

• より計算時間を長くしても, このまま安定性を保てるかどうか確認してみましょう. 計算時間を変えるためには, 先ほど出てきたように maxrevo の値を変えればよいのですが, 実は計算を終了させる条件がプログラム中に複数存在し, maxrevo の値を増やしただけでは計算時間が変わらない場合もあります. 計算時間を変えるためには, maxrevo 以外にプログラム中のどのパラメータを書き換えればよいでしょうか. 計算を終了させる条件について着目し, もう一度プログラムの中身を確認してみましょう.
• プログラムを 4 次のルンゲクッタ法で書き換えてみましょう. 単純な数値計算の問題ですがぜひやってみてください. 公式は参考書などで調べてみましょう.
• 保存則が成立しているのかを確認してみましょう. 数値計算をする上で, 系のエネルギー, 運動量, 角運動量がどうなっているか を検討するのは大切なことです. 理論の上では, 惑星の軌道では系の全エネル ギー, 全運動量, 角運動量は保存するはずですが, 実際に計算をしてみるとど うなるでしょうか. それぞれの保存量を求めるプログラムを書き加えて計算し て見ましょう.
• さらに惑星の数を増やした計算を扱ってみましょう. 3body.f90 は配列を用いて書かれています. なので比較的簡単に, 3 体以上の 惑星を扱うように書き換えることができるでしょう.
  • 現在の太陽系には 8 個の惑星があります. 今回もプログラムをうまく使い, 惑星の質量, 初期速度を適当な値に取り, 現在の太陽系の姿を再現できるか を確かめて見ましょう. また, もし再現ができないとすると, 何が原因かも 考えてみましょう.

本日の実習は以上です. お疲れ様でした.

• 2006 - 2008 年度 惑星科学実習資料 (計算機)
  • 本実習で用いているプログラム 3body.f90, 3body_sol.f90 は 2006 - 2008 年度惑星科学実習 (計算機) 用として、今枝佑輔氏が作成されたものを基にして作成されました.
• 2009 年度 地球および惑星大気科学実習 「数値計算実習 その2」 実習資料
• 木下宙 「天体と軌道の力学」 東京大学出版会

スケジュール表・各回資料 (07/09)