IT pass HikiWiki - [Exp2011]実際に計算してみる Diff
- Added parts are displayed like this.
- Deleted parts are displayed
like this.
((<"スケジュール表・各回資料 (07/08)"|[Exp2011]スケジュール表・各回資料#07-2F08>))
{{toc}}
= 実習の準備
今回の実習では, プログラムをゼロから作るのではなく, すでに途中まで出来上がっているプログラムを用いて行います. そのためのプログラムのソースファイルはすでに Web 上に置いてあります. なので, まずは情報実験機に適当なディレクトリを作って以下のファイルをダウンロードしてください.
* ((<2body1.f90|URL:http://itpass.scitec.kobe-u.ac.jp/~itpass/exp/fy2011/110708/practice/hiki-sample/2body1.f90>))
* ((<2body1_sol.f90|URL:http://itpass.scitec.kobe-u.ac.jp/~itpass/exp/fy2011/110708/practice/hiki-sample/2body1_sol.f90>))
* ((<2body2.f90|URL:http://itpass.scitec.kobe-u.ac.jp/~itpass/exp/fy2011/110708/practice/hiki-sample/2body2.f90>))
* ((<3body.f90|URL:http://itpass.scitec.kobe-u.ac.jp/~itpass/exp/fy2011/110708/practice2/hiki-sample/3body.f90>))
* ((<3body_sol.f90|URL:http://itpass.scitec.kobe-u.ac.jp/~itpass/exp/fy2011/110708/practice2/hiki-sample/3body_sol.f90>))
= 実習その１(二体問題)
== 実習その１−１ プログラムの穴埋めをしよう
ではダウンロードしたファイルのうちから 2body1.f90 を使ってみましょう. 開いてみると, Fortran のプログラムが書かれた文がならんでいます. 大まかに見ると, これらは次のような順番で文が書かれています.
* 使う変数の宣言する.
* いくつかの変数に値を与える.
* 実行する計算を設定する.
* 値の書き出しは随時行う.
よろしいでしょうか. では早速これらをコンパイルしてみましょう. ターミナルを開き, ファイルをダウンロードしたディレクトリに移動してから, 次のように入力してください.
$ gfortran 2body1.f90
すると, 以下のような内容がターミナル上に出力されるでしょう.
2body1.f90:55.16:
rp = sqrt(??????????)
1
Error: Syntax error in argument list at (1)
これは,
2body1.f90 というファイルの 55 行目の 16 番目の文字から後ろ, すなわち 1 で指す部分から先に, 本来 Fortarn で書かれるべき形式で書かれていない所がある.
という意味です. なので, この「??????????」に適当な文を入れてやる必要があります. ここ以外にも同様のエラー合計6つが出ているので, 先程のレクチャーを思い出しながら, 適当な文を入れていき, エラーがでないようにして下さい.
== 実習その１−２ プログラムを実行し, 惑星の軌道を図示しましょう.
プログラムをエラーなくコンパイルできたら, コンパイルされたプログラムを実行してみましょう.問題なくコンパイルが終わった後、ディレクトリに 「 a.out 」というファイルができているでしょう. それを確認した後, 次のように入力してください.
$ ./a.out
すると, ターミナル上に
xp= 1.0000 yp= 0.0000
vx= 0.0000 vy= 1.1000
という出力があった後, ディレクトリの中に「 fort.12 」というファイルができているでしょう. ちなみに, なぜターミナル上に xp, yp, vx, vy が出力され, ディレクトリに fort.12 ができたのでしょうか. 無論それを作るようなプログラムが 2body1.f90 に書かれているのですが, どの部分にかかれているかを確認しておきましょう.
さて次に, 作った fort.12 を用いて, 実際に絵を書いてみましょう. ここで Gnuplot を立ち上げます. そのために今お使いのターミナルとは別のターミナルを立ち上げて下さい. その後, fort.12 があるディレクトリに移動し, 次のように入力しましょう.
$ gnuplot
すると, Gnuplot が立ち上がりました. 詳しい使い方は((<こちら|[Exp2011]gnuplot を使ってみよう>))を見て下さい。では, Gnuplot を立ち上げたターミナルで次のように入力してください.
gnuplot>plot "fort.12" u 2:3
このとき楕円軌道が描けているでしょうか. あるいは
gnuplot>plot "fort.12" u 2:3 w l
としたらどうなるでしょうか. 試してみてください.
※ どうしてもコンパイルがうまくいかないとき, もしくは楕円軌道がうまく描けないときは, 2body1_sol.f90 を使ってください.
== 実習その１−３ 2body1.f90 と 2body2.f90 の実行結果を比較し, 計算精度を確認してみましょう.
これまで皆さんに扱っていただいていた 2body1.f90 には, オイラー法が使わ
れています. この計算方法は, 実はそれほど精度がよい計算方法ではありません. 例えば 2body1.f90 をエディタで開き
do step = 1, 100000
として実行してください. すると, 軌道がどんどん膨らんでいっているのが分かります. では次に 2body2.f90 を使ってみましょう. これには２次のルンゲクッタ法が使われています. これについても同様に,
do step = 1, 100000
として実行してください. このとき, 実行した結果できるファイルは fort.13 になります. 以下同じように Gnuplot で絵を描き出すとどうなるでしょうか.
== 実習その１−４ 惑星の初期位置や初期速度を変えたとき, 惑星の運動はどうなるかを確認しましょう.
ここまででプログラムの仕組みはなんとなく分かりましたか？ では初期条件をいろいろ変えてみましょう. 初期条件を与えているのは, プログラム中の次の部分です.
!---- 初期に惑星を(1,0)に配置する
xp = 1.0d0
yp = 0.0d0
!---- おいた惑星を ypの正の方向に初速を与えて放り出す
vx = 0.0d0
vy = 1.1d0
ここを適当に変えて, 惑星の軌道がどうなるかを確認してください.
= 実習 その 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 を使ってください.
= 実習 その 3
== 初期条件の変更
* 惑星の質量を変化させるとどうなるでしょうか.
* 惑星の初期位置を変化させるとどうなるでしょうか.
初期条件を次のように変化させてみましょう.
* 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) を計算し終えた人から, 計算結果を記入しに来て下さい.
計算した結果, 全ての天体の軌道が安定であった場合は○を, 不安定なものがあった場合は×を表の中に記入してください.
本日の実習は以上です. お疲れ様でした.
= 参考資料
* 2006 - 2008 年度 惑星科学実習資料 (計算機)
* ((<2010 第11回ITPASS実習資料 「数値計算実習 その1」|[Exp2010]実際に計算してみる (２体問題)>))
* 本実習で用いているプログラム 2body1.f90, 2body1_sol.f90, 2body2.f90 は 2006 - 2008 年度惑星科学実習 (計算機) 用として、今枝佑輔氏が作成されたものを基にして作成されました.
* ((<2010 第12回ITPASS実習資料 「数値計算実習 その2」|[Exp2010]実際に計算してみる (３体問題)>))
* 本実習で用いているプログラム 3body.f90, 3body_sol.f90 は 2006 - 2008 年度惑星科学実習 (計算機) 用として、今枝佑輔氏が作成されたものを基にして作成されました.
* 木下宙 「天体と軌道の力学」 東京大学出版会
((<"スケジュール表・各回資料 (07/08)"|[Exp2011]スケジュール表・各回資料#07-2F08>))
{{toc}}
= 実習の準備
今回の実習では, プログラムをゼロから作るのではなく, すでに途中まで出来上がっているプログラムを用いて行います. そのためのプログラムのソースファイルはすでに Web 上に置いてあります. なので, まずは情報実験機に適当なディレクトリを作って以下のファイルをダウンロードしてください.
* ((<2body1.f90|URL:http://itpass.scitec.kobe-u.ac.jp/~itpass/exp/fy2011/110708/practice/hiki-sample/2body1.f90>))
* ((<2body1_sol.f90|URL:http://itpass.scitec.kobe-u.ac.jp/~itpass/exp/fy2011/110708/practice/hiki-sample/2body1_sol.f90>))
* ((<2body2.f90|URL:http://itpass.scitec.kobe-u.ac.jp/~itpass/exp/fy2011/110708/practice/hiki-sample/2body2.f90>))
* ((<3body.f90|URL:http://itpass.scitec.kobe-u.ac.jp/~itpass/exp/fy2011/110708/practice2/hiki-sample/3body.f90>))
* ((<3body_sol.f90|URL:http://itpass.scitec.kobe-u.ac.jp/~itpass/exp/fy2011/110708/practice2/hiki-sample/3body_sol.f90>))
= 実習その１(二体問題)
== 実習その１−１ プログラムの穴埋めをしよう
ではダウンロードしたファイルのうちから 2body1.f90 を使ってみましょう. 開いてみると, Fortran のプログラムが書かれた文がならんでいます. 大まかに見ると, これらは次のような順番で文が書かれています.
* 使う変数の宣言する.
* いくつかの変数に値を与える.
* 実行する計算を設定する.
* 値の書き出しは随時行う.
よろしいでしょうか. では早速これらをコンパイルしてみましょう. ターミナルを開き, ファイルをダウンロードしたディレクトリに移動してから, 次のように入力してください.
$ gfortran 2body1.f90
すると, 以下のような内容がターミナル上に出力されるでしょう.
2body1.f90:55.16:
rp = sqrt(??????????)
1
Error: Syntax error in argument list at (1)
これは,
2body1.f90 というファイルの 55 行目の 16 番目の文字から後ろ, すなわち 1 で指す部分から先に, 本来 Fortarn で書かれるべき形式で書かれていない所がある.
という意味です. なので, この「??????????」に適当な文を入れてやる必要があります. ここ以外にも同様のエラー合計6つが出ているので, 先程のレクチャーを思い出しながら, 適当な文を入れていき, エラーがでないようにして下さい.
== 実習その１−２ プログラムを実行し, 惑星の軌道を図示しましょう.
プログラムをエラーなくコンパイルできたら, コンパイルされたプログラムを実行してみましょう.問題なくコンパイルが終わった後、ディレクトリに 「 a.out 」というファイルができているでしょう. それを確認した後, 次のように入力してください.
$ ./a.out
すると, ターミナル上に
xp= 1.0000 yp= 0.0000
vx= 0.0000 vy= 1.1000
という出力があった後, ディレクトリの中に「 fort.12 」というファイルができているでしょう. ちなみに, なぜターミナル上に xp, yp, vx, vy が出力され, ディレクトリに fort.12 ができたのでしょうか. 無論それを作るようなプログラムが 2body1.f90 に書かれているのですが, どの部分にかかれているかを確認しておきましょう.
さて次に, 作った fort.12 を用いて, 実際に絵を書いてみましょう. ここで Gnuplot を立ち上げます. そのために今お使いのターミナルとは別のターミナルを立ち上げて下さい. その後, fort.12 があるディレクトリに移動し, 次のように入力しましょう.
$ gnuplot
すると, Gnuplot が立ち上がりました. 詳しい使い方は((<こちら|[Exp2011]gnuplot を使ってみよう>))を見て下さい。では, Gnuplot を立ち上げたターミナルで次のように入力してください.
gnuplot>plot "fort.12" u 2:3
このとき楕円軌道が描けているでしょうか. あるいは
gnuplot>plot "fort.12" u 2:3 w l
としたらどうなるでしょうか. 試してみてください.
※ どうしてもコンパイルがうまくいかないとき, もしくは楕円軌道がうまく描けないときは, 2body1_sol.f90 を使ってください.
== 実習その１−３ 2body1.f90 と 2body2.f90 の実行結果を比較し, 計算精度を確認してみましょう.
これまで皆さんに扱っていただいていた 2body1.f90 には, オイラー法が使わ
れています. この計算方法は, 実はそれほど精度がよい計算方法ではありません. 例えば 2body1.f90 をエディタで開き
do step = 1, 100000
として実行してください. すると, 軌道がどんどん膨らんでいっているのが分かります. では次に 2body2.f90 を使ってみましょう. これには２次のルンゲクッタ法が使われています. これについても同様に,
do step = 1, 100000
として実行してください. このとき, 実行した結果できるファイルは fort.13 になります. 以下同じように Gnuplot で絵を描き出すとどうなるでしょうか.
== 実習その１−４ 惑星の初期位置や初期速度を変えたとき, 惑星の運動はどうなるかを確認しましょう.
ここまででプログラムの仕組みはなんとなく分かりましたか？ では初期条件をいろいろ変えてみましょう. 初期条件を与えているのは, プログラム中の次の部分です.
!---- 初期に惑星を(1,0)に配置する
xp = 1.0d0
yp = 0.0d0
!---- おいた惑星を ypの正の方向に初速を与えて放り出す
vx = 0.0d0
vy = 1.1d0
ここを適当に変えて, 惑星の軌道がどうなるかを確認してください.
= 実習 その 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 を使ってください.
= 実習 その 3
== 初期条件の変更
* 惑星の質量を変化させるとどうなるでしょうか.
* 惑星の初期位置を変化させるとどうなるでしょうか.
初期条件を次のように変化させてみましょう.
* 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) を計算し終えた人から, 計算結果を記入しに来て下さい.
計算した結果, 全ての天体の軌道が安定であった場合は○を, 不安定なものがあった場合は×を表の中に記入してください.
本日の実習は以上です. お疲れ様でした.
= 参考資料
* 2006 - 2008 年度 惑星科学実習資料 (計算機)
* ((<2010 第11回ITPASS実習資料 「数値計算実習 その1」|[Exp2010]実際に計算してみる (２体問題)>))
* 本実習で用いているプログラム 2body1.f90, 2body1_sol.f90, 2body2.f90 は 2006 - 2008 年度惑星科学実習 (計算機) 用として、今枝佑輔氏が作成されたものを基にして作成されました.
* ((<2010 第12回ITPASS実習資料 「数値計算実習 その2」|[Exp2010]実際に計算してみる (３体問題)>))
* 本実習で用いているプログラム 3body.f90, 3body_sol.f90 は 2006 - 2008 年度惑星科学実習 (計算機) 用として、今枝佑輔氏が作成されたものを基にして作成されました.
* 木下宙 「天体と軌道の力学」 東京大学出版会
((<"スケジュール表・各回資料 (07/08)"|[Exp2011]スケジュール表・各回資料#07-2F08>))