IT pass HikiWiki - [Exp2009]実際に計算してみる (３体問題) Diff
- Added parts are displayed like this.
- Deleted parts are displayed
like this.
((<"スケジュール表・各回資料 (07/17)"|[Exp2009]スケジュール表・各回資料#07-2F17>))
{{toc}}
= 実習の準備
本日の実習でも前回と同様に, すでに途中まで出来上がっているプログラムに手を加えて作業を行います. なので, 前回と同じ要領で情報実験機に以下のファイルをダウンロードしてください.
* ((<3body.f90|URL:http://itpass.scitec.kobe-u.ac.jp/~itpass/exp/fy2009/090717/practice/hiki-sample/3body.f90>))
* ((<3body_sol.f90|URL:http://itpass.scitec.kobe-u.ac.jp/~itpass/exp/fy2009/090717/practice/hiki-sample/3body_sol.f90>))
ちなみに, 今回用いるプログラムには全て２次のルンゲクッタ法が使われています.
= 実習その２
== 実習その２−１ プログラムの穴埋めをしましょう
ダウンロードしたファイルのうちから 3body.f90 を使ってみましょう. 開いてみると, Fortran のプログラムが書かれた文がならんでいます. これにも前回と同様に, プログラムが正しく書かれていないところがあります. 前回と同じ要領で穴埋めしていきましょう.
== 初期設定について
* 天体の質量
３体問題を扱うにあたり, ２体問題では取り扱わなかった質量の設定も初期の
設定に含まれています. それはここにあたります.
!---- 各惑星の質量を決める
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)= 9.5546d0/5.2026d0
yp(3)= 0.0d0
これも太陽系の場合にあたります. ここも後で変えて頂いてかまいませんが、
変えるのはなるべく, xp(3) だけにしてください.
* 惑星の初期速度
惑星の初期速度は次のように与えられています.
!---- 初期に惑星を次のように配置する.
xp(1)= 0.0d0
yp(1)= 0.0d0
xp(2)= 0.0d0
yp(2)= sqrt(mp(1)+mp(2))/sqrt(abs(xp(2)))
xp(3)= 0.0d0
yp(3)= sqrt(mp(1)+mp(3))/sqrt(abs(xp(3)))
これも太陽系の場合にあたります. 物理的意味としては, 惑星を y 方向に
だいたいケプラー速度で飛ばしたということになります. 先の初期位置の y
座標を書き換えるならば, ここも適当に変える必要があります.
興味があればどうぞやってみてください.
== 実習その２−２ 太陽系は安定かどうか確認してみましょう
ファイル 3body.f90 はデフォルトでは太陽系の場合を表しています. これを用いて, 太陽系が本当に安定かどうかを考えてみましょう. その際, プログラム中のパラメータを次のようにしてください.
nstep = 1.0d8
maxtime = 3.0d5
上記のようにパラメータを設定して計算することは, 木星が 50000 回公転する期間を計算することに相当します.
ちなみに, ここで言う安定とは, 「計算している範囲で惑星の軌道が楕円形を保ち続ける」という意味だと思ってください.
※ どうしてもコンパイルがうまくいかない, もしくは軌道の図がうまく描けないときは, 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 を丁度, ８個用意したので, 各実験機につき, factor 1つ分をこなしてください.
この結果をレポートにまとめて提出してください. 詳細は((<こちら|URL:https://epa.scitec.kobe-u.ac.jp/~itpass/hiki/hiki.cgi?%5BExp2009%5D%BF%F4%C3%CD%B7%D7%BB%BB%BC%C2%BD%AC%B2%DD%C2%EA%A4%BD%A4%CE2>))です.
= Extra
今回の実習が物足りない, と思われる方は次のテーマに取り組んでみましょう.
* プログラムを 4 次のルンゲクッタ法で書き換えてみましょう.
単純な数値計算の問題ですがぜひやってみましょう. 公式は参考書などで調べ
てみてください.
* 保存則が成立しているのかを確認してみましょう.
数値計算をする上で, 系のエネルギー, 運動量, 角運動量がどうなっているか
を検討するのは大切なことです. 理論の上では, 惑星の軌道では系の全エネル
ギー, 全運動量, 角運動量は保存するはずですが, 実際に計算をしてみるとど
うなるでしょうか. それぞれの保存量を求めるプログラムを書き加えて計算し
て見ましょう.
* 3 体以上の惑星を扱ってみましょう.
3body.f90 は配列を用いて書かれています. なので比較的簡単に, 3 体以上の
惑星を扱うように書き換えることができるでしょう.
* 天体力学では, 3 体以上の惑星があると惑星の軌道は不安定になるというこ
とが知られています. ( ((<下記参考資料参照|URL:#参考資料>)) )nstep や
maxtime を変更したり, 書き出しの条件の下記の if 文
if(mod(step, 50000).eq.1) then
の 50000 や
if(const.eq.5000) Goto 2001
の 5000 を書き換えたりして, 調べてみましょう.
* 現在の太陽系には 8 個の惑星があります. 今回もプログラムをうまく使い,
惑星の質量, 初期速度を適当な値に取り, 現在の太陽系の姿を再現できるか
を確かめて見ましょう. また, もし再現ができないとすると, 何が原因かも
考えてみましょう.
本日の実習は以上です. お疲れ様でした.
== 参考資料
* 2006 - 2008 年度 惑星科学実習資料 (計算機)
* 本実習の実習資料は 2006 - 2008 年度 惑星科学実習資料 (計算機)を基に作成されました.
* 本実習で用いているプログラム 3body.f90, 3body_sol.f90 は 2006 - 2008 年度惑星科学実習 (計算機) 用として、今枝佑輔氏が作成されたものを基にして作成されました.
* 異形の惑星 系外惑星形成論から (井田 茂 著 / NHK BOOKS)
* The Stability of Multi-Planet Systems
(J. E. Chambers, G. W. Wetherill, A. P. Boss, /
Icarus, Volume 119, Issue 2, February 1996, Pages 261-268)
((<"スケジュール表・各回資料 (07/17)"|[Exp2009]スケジュール表・各回資料#07-2F17>))
{{toc}}
= 実習の準備
本日の実習でも前回と同様に, すでに途中まで出来上がっているプログラムに手を加えて作業を行います. なので, 前回と同じ要領で情報実験機に以下のファイルをダウンロードしてください.
* ((<3body.f90|URL:http://itpass.scitec.kobe-u.ac.jp/~itpass/exp/fy2009/090717/practice/hiki-sample/3body.f90>))
* ((<3body_sol.f90|URL:http://itpass.scitec.kobe-u.ac.jp/~itpass/exp/fy2009/090717/practice/hiki-sample/3body_sol.f90>))
ちなみに, 今回用いるプログラムには全て２次のルンゲクッタ法が使われています.
= 実習その２
== 実習その２−１ プログラムの穴埋めをしましょう
ダウンロードしたファイルのうちから 3body.f90 を使ってみましょう. 開いてみると, Fortran のプログラムが書かれた文がならんでいます. これにも前回と同様に, プログラムが正しく書かれていないところがあります. 前回と同じ要領で穴埋めしていきましょう.
== 初期設定について
* 天体の質量
３体問題を扱うにあたり, ２体問題では取り扱わなかった質量の設定も初期の
設定に含まれています. それはここにあたります.
!---- 各惑星の質量を決める
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)= 9.5546d0/5.2026d0
yp(3)= 0.0d0
これも太陽系の場合にあたります. ここも後で変えて頂いてかまいませんが、
変えるのはなるべく, xp(3) だけにしてください.
* 惑星の初期速度
惑星の初期速度は次のように与えられています.
!---- 初期に惑星を次のように配置する.
xp(1)= 0.0d0
yp(1)= 0.0d0
xp(2)= 0.0d0
yp(2)= sqrt(mp(1)+mp(2))/sqrt(abs(xp(2)))
xp(3)= 0.0d0
yp(3)= sqrt(mp(1)+mp(3))/sqrt(abs(xp(3)))
これも太陽系の場合にあたります. 物理的意味としては, 惑星を y 方向に
だいたいケプラー速度で飛ばしたということになります. 先の初期位置の y
座標を書き換えるならば, ここも適当に変える必要があります.
興味があればどうぞやってみてください.
== 実習その２−２ 太陽系は安定かどうか確認してみましょう
ファイル 3body.f90 はデフォルトでは太陽系の場合を表しています. これを用いて, 太陽系が本当に安定かどうかを考えてみましょう. その際, プログラム中のパラメータを次のようにしてください.
nstep = 1.0d8
maxtime = 3.0d5
上記のようにパラメータを設定して計算することは, 木星が 50000 回公転する期間を計算することに相当します.
ちなみに, ここで言う安定とは, 「計算している範囲で惑星の軌道が楕円形を保ち続ける」という意味だと思ってください.
※ どうしてもコンパイルがうまくいかない, もしくは軌道の図がうまく描けないときは, 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 を丁度, ８個用意したので, 各実験機につき, factor 1つ分をこなしてください.
この結果をレポートにまとめて提出してください. 詳細は((<こちら|URL:https://epa.scitec.kobe-u.ac.jp/~itpass/hiki/hiki.cgi?%5BExp2009%5D%BF%F4%C3%CD%B7%D7%BB%BB%BC%C2%BD%AC%B2%DD%C2%EA%A4%BD%A4%CE2>))です.
= Extra
今回の実習が物足りない, と思われる方は次のテーマに取り組んでみましょう.
* プログラムを 4 次のルンゲクッタ法で書き換えてみましょう.
単純な数値計算の問題ですがぜひやってみましょう. 公式は参考書などで調べ
てみてください.
* 保存則が成立しているのかを確認してみましょう.
数値計算をする上で, 系のエネルギー, 運動量, 角運動量がどうなっているか
を検討するのは大切なことです. 理論の上では, 惑星の軌道では系の全エネル
ギー, 全運動量, 角運動量は保存するはずですが, 実際に計算をしてみるとど
うなるでしょうか. それぞれの保存量を求めるプログラムを書き加えて計算し
て見ましょう.
* 3 体以上の惑星を扱ってみましょう.
3body.f90 は配列を用いて書かれています. なので比較的簡単に, 3 体以上の
惑星を扱うように書き換えることができるでしょう.
* 天体力学では, 3 体以上の惑星があると惑星の軌道は不安定になるというこ
とが知られています. ( ((<下記参考資料参照|URL:#参考資料>)) )nstep や
maxtime を変更したり, 書き出しの条件の下記の if 文
if(mod(step, 50000).eq.1) then
の 50000 や
if(const.eq.5000) Goto 2001
の 5000 を書き換えたりして, 調べてみましょう.
* 現在の太陽系には 8 個の惑星があります. 今回もプログラムをうまく使い,
惑星の質量, 初期速度を適当な値に取り, 現在の太陽系の姿を再現できるか
を確かめて見ましょう. また, もし再現ができないとすると, 何が原因かも
考えてみましょう.
本日の実習は以上です. お疲れ様でした.
== 参考資料
* 2006 - 2008 年度 惑星科学実習資料 (計算機)
* 本実習の実習資料は 2006 - 2008 年度 惑星科学実習資料 (計算機)を基に作成されました.
* 本実習で用いているプログラム 3body.f90, 3body_sol.f90 は 2006 - 2008 年度惑星科学実習 (計算機) 用として、今枝佑輔氏が作成されたものを基にして作成されました.
* 異形の惑星 系外惑星形成論から (井田 茂 著 / NHK BOOKS)
* The Stability of Multi-Planet Systems
(J. E. Chambers, G. W. Wetherill, A. P. Boss, /
Icarus, Volume 119, Issue 2, February 1996, Pages 261-268)
((<"スケジュール表・各回資料 (07/17)"|[Exp2009]スケジュール表・各回資料#07-2F17>))