IT pass HikiWiki - [Exp2011]実際に計算してみる (円制限 3 体問題) Diff
- Added parts are displayed like this.
- Deleted parts are displayed
like this.
((<"スケジュール表・各回資料 (07/15)"|[Exp2011]スケジュール表・各回資料#07-2F15>))
{{toc}}
= 実習の準備
本日の実習でも前回と同様に, すでに途中まで出来上がっているプログラムに手を加えて作業を行います. 前回と同じ要領で情報実験機に以下のファイルをダウンロードしてください.
* ((<res-3body_4runge.f90|URL:http://itpass.scitec.kobe-u.ac.jp/~itpass/exp/fy2011/110715/practice/hiki-sample/res-3body_4runge.f90>))
* ((<res-3body_4runge_sol.f90|URL:http://itpass.scitec.kobe-u.ac.jp/~itpass/exp/fy2011/110715/practice/hiki-sample/res-3body_4runge_sol.f90>))
* ((<res-3body_2runge.f90|URL:http://itpass.scitec.kobe-u.ac.jp/~itpass/exp/fy2011/110715/practice/hiki-sample/res-3body_2runge.f90>))
なお, 今回用いるプログラムには 2 次のルンゲクッタ法と 4 次のルンゲクッタ法が使われています.
= 実習 その 4
== プログラムの穴埋め
ダウンロードした 3 つのファイルのうち, res-3body_4runge.f90 を開いてみましょう. 中には Fortran のプログラムが書かれた文がならんでいます.
しかし前回と同様に, プログラムが正しく書かれていないところがあります. 前回と同じ要領で穴埋めしていきましょう.
このプログラムは 4 次のルンゲクッタ法で書かれていますが, 穴埋めが必要な箇所が全部で 16 行あります. レクチャー編の資料を参考にしながら全ての箇所に正しい式を書き込み, コンパイル時にエラーが出ないようにしてください. Fortran を書くうえで分からないことがあれば((<Fortran の簡単な解説|[Exp2011]Fortran の簡単な解説>))を参考にしてください. コンパイルができると, 次のような表示が端末に出力されます.
$ gfortran res-3body_4runge.f90
$ ./a.out
factor
mu1= 0.99904640816924728 mu2= 9.53591830752668532E-004
x1= -9.53591830752668532E-004 y1= 0.0000000000000000
x2= 0.99904640816924728 y2= 0.0000000000000000
initial condition
x = 0.28754853342559489 y = 0.0000000000000000
vx= 0.0000000000000000 vy= 1.5700000000000001
Cj= 2.2731059913538694
factor は中心星と惑星の質量比と位置を, initial condition は 3 つ目の天体 (粒子) の初期条件を表しています. ここでは, 中心星を太陽, 惑星を木星としています.
粒子の軌道は 4orbit.dat に, 中心星と惑星の位置は locate.datに, ヤコビ定数から定まる運動の禁止領域は forbidden.dat に出力されています.
前回と同様に Gnuplot を立ち上げて, 軌道を描いてみましょう.
$ gnuplot
gnuplot> set size square
gnuplot> plot "forbidden.dat", "4orbit.dat" w l, "locate.dat" pt 7 ps 2
== 初期条件について
* 天体の質量
中心星と惑星の質量は次のように与えられています.
!--中心星と惑星の質量比を求める-----------------------------------
m1 = 1.9891d33
m2 = 1.8986d30
mu = m2/(m1+m2)
mu1 = 1.0d0-mu
mu2 = mu
mu1 と mu2 が中心星と惑星の質量比を表しており, ここでは太陽と木星の関係になっています. 円制限三体問題では, 中心星と惑星に比べて粒子の質量は無視できるほど小さいとしているので, 条件には入ってきません.
* 粒子の初期位置と初速度
粒子の初期位置と初速度は次のように与えられています.
!--初期条件---------------------------------------------------
!--粒子の初期位置および初速度を与える----------------------
x = 1.4960d0 / 5.2026d0
y = 0.0d0
vx = 0.0d0
vy = 1.57d0
ここでは, 粒子がほぼ地球軌道をとるように初期条件を与えています.
※ どうしてもコンパイルがうまくいかないとき, もしくは軌道がうまく描けないときは, res-3body_4runge_sol.f90 を使ってください.
= 実習 その 5
== 計算精度の比較
これまで皆さんに扱っていただいていた res-3body_4runge.f90 には, 4 次のルンゲクッタ法が使われています. ここでは, 前回の実習で使われていた 2 次のルンゲクッタ法で計算したときと比較をしてみましょう. それぞれコンパイルをするのですが, 初期条件は同じにして, 終了する時間を
maxrevo = 2.0d0*PI * 1.0d2
として計算を行ってください. これは, 惑星が中心星の周りを100回公転する時間に相当します.
res-3body_2runge.f90を実行した結果できるファイルは 2runge.dat になります. 同じように Gnuplot で軌道を描くとどうなるでしょうか.
また, ヤコビ定数が本当に保存されているかを見て, それぞれの計算精度を確認してみましょう. 出力されるファイル名は, 2 次は 2jacobi.dat, 4 次は 4jacobi.dat です.
gnuplot>plot "2jacobi.dat", "4jacobi.dat"
として, プロットしてみましょう. また,
gnuplot>set log y
gnuplot>plot "2jacobi.dat" u 1:3,"4jacobi.dat" u 1:3
として, 誤差をみてみましょう.
== 解析解と数値解
レクチャーの中で, 回転系で粒子が静止しているような解析解であるラグランジュ平衡点の話をしました. ここでは, 数値計算を行って本当に静止しているか確認してみましょう.
中心星−惑星−粒子の配置が正三角形で, 速度がゼロであるように初期条件を与えて計算を行ってみてください.
粒子が静止している場合は 1 点だけしかプロットされないので, 線ではなく点でプロットしましょう.
$ gnuplot
gnuplot> set size square
gnuplot> plot "forbidden.dat", "4orbit.dat", "locate.dat" pt 7 ps 2
= 実習 その 6
== 初期条件の変更
* 中心星と惑星の質量を変化させるとどうなるでしょうか.
* 粒子の初期位置や初速度を変化させるとどうなるでしょうか.
円制限 3 体問題では 3 つ目の天体 (粒子) の質量は変えませんが, 中心星と惑星の質量や粒子の初期条件を変えることで様々な軌道が描かれます. パラメータを変えて軌道を 3 種類以上描き, それぞれの軌道がどのようなものか考えてみましょう. 必要に応じて, 計算時間 maxrevo の値も変えてください.
さらに, この結果についてレポートにまとめ, 7/01(金) に出題した((<「課題 その1」|[Exp2011]数値計算実習課題その1>))と一緒に提出してください. 課題その 2 の詳細は((<こちら|[Exp2011]数値計算実習課題その2>))です.
本日の実習は以上です. お疲れ様でした.
= 参考資料
* ((<2010 年度 ITPASS実習 「数値計算実習 その2」 実習資料|[Exp2010]実際に計算してみる (３体問題)>))
* 木下宙 「天体と軌道の力学」 東京大学出版会
* C.D. Murray and S.F. Dermott 「Solar System Dynamics」 CAMBRIDGE UNIVERSITY PRESS
((<"スケジュール表・各回資料 (07/15)"|[Exp2011]スケジュール表・各回資料#07-2F15>))
{{toc}}
= 実習の準備
本日の実習でも前回と同様に, すでに途中まで出来上がっているプログラムに手を加えて作業を行います. 前回と同じ要領で情報実験機に以下のファイルをダウンロードしてください.
* ((<res-3body_4runge.f90|URL:http://itpass.scitec.kobe-u.ac.jp/~itpass/exp/fy2011/110715/practice/hiki-sample/res-3body_4runge.f90>))
* ((<res-3body_4runge_sol.f90|URL:http://itpass.scitec.kobe-u.ac.jp/~itpass/exp/fy2011/110715/practice/hiki-sample/res-3body_4runge_sol.f90>))
* ((<res-3body_2runge.f90|URL:http://itpass.scitec.kobe-u.ac.jp/~itpass/exp/fy2011/110715/practice/hiki-sample/res-3body_2runge.f90>))
なお, 今回用いるプログラムには 2 次のルンゲクッタ法と 4 次のルンゲクッタ法が使われています.
= 実習 その 4
== プログラムの穴埋め
ダウンロードした 3 つのファイルのうち, res-3body_4runge.f90 を開いてみましょう. 中には Fortran のプログラムが書かれた文がならんでいます.
しかし前回と同様に, プログラムが正しく書かれていないところがあります. 前回と同じ要領で穴埋めしていきましょう.
このプログラムは 4 次のルンゲクッタ法で書かれていますが, 穴埋めが必要な箇所が全部で 16 行あります. レクチャー編の資料を参考にしながら全ての箇所に正しい式を書き込み, コンパイル時にエラーが出ないようにしてください. Fortran を書くうえで分からないことがあれば((<Fortran の簡単な解説|[Exp2011]Fortran の簡単な解説>))を参考にしてください. コンパイルができると, 次のような表示が端末に出力されます.
$ gfortran res-3body_4runge.f90
$ ./a.out
factor
mu1= 0.99904640816924728 mu2= 9.53591830752668532E-004
x1= -9.53591830752668532E-004 y1= 0.0000000000000000
x2= 0.99904640816924728 y2= 0.0000000000000000
initial condition
x = 0.28754853342559489 y = 0.0000000000000000
vx= 0.0000000000000000 vy= 1.5700000000000001
Cj= 2.2731059913538694
factor は中心星と惑星の質量比と位置を, initial condition は 3 つ目の天体 (粒子) の初期条件を表しています. ここでは, 中心星を太陽, 惑星を木星としています.
粒子の軌道は 4orbit.dat に, 中心星と惑星の位置は locate.datに, ヤコビ定数から定まる運動の禁止領域は forbidden.dat に出力されています.
前回と同様に Gnuplot を立ち上げて, 軌道を描いてみましょう.
$ gnuplot
gnuplot> set size square
gnuplot> plot "forbidden.dat", "4orbit.dat" w l, "locate.dat" pt 7 ps 2
== 初期条件について
* 天体の質量
中心星と惑星の質量は次のように与えられています.
!--中心星と惑星の質量比を求める-----------------------------------
m1 = 1.9891d33
m2 = 1.8986d30
mu = m2/(m1+m2)
mu1 = 1.0d0-mu
mu2 = mu
mu1 と mu2 が中心星と惑星の質量比を表しており, ここでは太陽と木星の関係になっています. 円制限三体問題では, 中心星と惑星に比べて粒子の質量は無視できるほど小さいとしているので, 条件には入ってきません.
* 粒子の初期位置と初速度
粒子の初期位置と初速度は次のように与えられています.
!--初期条件---------------------------------------------------
!--粒子の初期位置および初速度を与える----------------------
x = 1.4960d0 / 5.2026d0
y = 0.0d0
vx = 0.0d0
vy = 1.57d0
ここでは, 粒子がほぼ地球軌道をとるように初期条件を与えています.
※ どうしてもコンパイルがうまくいかないとき, もしくは軌道がうまく描けないときは, res-3body_4runge_sol.f90 を使ってください.
= 実習 その 5
== 計算精度の比較
これまで皆さんに扱っていただいていた res-3body_4runge.f90 には, 4 次のルンゲクッタ法が使われています. ここでは, 前回の実習で使われていた 2 次のルンゲクッタ法で計算したときと比較をしてみましょう. それぞれコンパイルをするのですが, 初期条件は同じにして, 終了する時間を
maxrevo = 2.0d0*PI * 1.0d2
として計算を行ってください. これは, 惑星が中心星の周りを100回公転する時間に相当します.
res-3body_2runge.f90を実行した結果できるファイルは 2runge.dat になります. 同じように Gnuplot で軌道を描くとどうなるでしょうか.
また, ヤコビ定数が本当に保存されているかを見て, それぞれの計算精度を確認してみましょう. 出力されるファイル名は, 2 次は 2jacobi.dat, 4 次は 4jacobi.dat です.
gnuplot>plot "2jacobi.dat", "4jacobi.dat"
として, プロットしてみましょう. また,
gnuplot>set log y
gnuplot>plot "2jacobi.dat" u 1:3,"4jacobi.dat" u 1:3
として, 誤差をみてみましょう.
== 解析解と数値解
レクチャーの中で, 回転系で粒子が静止しているような解析解であるラグランジュ平衡点の話をしました. ここでは, 数値計算を行って本当に静止しているか確認してみましょう.
中心星−惑星−粒子の配置が正三角形で, 速度がゼロであるように初期条件を与えて計算を行ってみてください.
粒子が静止している場合は 1 点だけしかプロットされないので, 線ではなく点でプロットしましょう.
$ gnuplot
gnuplot> set size square
gnuplot> plot "forbidden.dat", "4orbit.dat", "locate.dat" pt 7 ps 2
= 実習 その 6
== 初期条件の変更
* 中心星と惑星の質量を変化させるとどうなるでしょうか.
* 粒子の初期位置や初速度を変化させるとどうなるでしょうか.
円制限 3 体問題では 3 つ目の天体 (粒子) の質量は変えませんが, 中心星と惑星の質量や粒子の初期条件を変えることで様々な軌道が描かれます. パラメータを変えて軌道を 3 種類以上描き, それぞれの軌道がどのようなものか考えてみましょう. 必要に応じて, 計算時間 maxrevo の値も変えてください.
さらに, この結果についてレポートにまとめ, 7/01(金) に出題した((<「課題 その1」|[Exp2011]数値計算実習課題その1>))と一緒に提出してください. 課題その 2 の詳細は((<こちら|[Exp2011]数値計算実習課題その2>))です.
本日の実習は以上です. お疲れ様でした.
= 参考資料
* ((<2010 年度 ITPASS実習 「数値計算実習 その2」 実習資料|[Exp2010]実際に計算してみる (３体問題)>))
* 木下宙 「天体と軌道の力学」 東京大学出版会
* C.D. Murray and S.F. Dermott 「Solar System Dynamics」 CAMBRIDGE UNIVERSITY PRESS
((<"スケジュール表・各回資料 (07/15)"|[Exp2011]スケジュール表・各回資料#07-2F15>))