IT pass HikiWiki - [Exp2010]実際に計算してみる (２体問題) Diff
- Added parts are displayed like this.
- Deleted parts are displayed
like this.
((<"スケジュール表・各回資料 (07/02)"|[Exp2010]スケジュール表・各回資料#07-2F02>))
{{toc}}
= 実習の準備
今回の実習では, プログラムをゼロから作るのではなく, すでに途中まで出来上がっているプログラムを用いて行います. そのためのプログラムのソースファイルはすでに Web 上に置いてあります. なので, まずは情報実験機に適当なディレクトリを作って以下のファイルをダウンロードしてください.
* ((<2body1.f90|URL:http://itpass.scitec.kobe-u.ac.jp/~itpass/exp/fy2010/100702/practice/hiki-sample/2body1.f90>))
* ((<2body1_sol.f90|URL:http://itpass.scitec.kobe-u.ac.jp/~itpass/exp/fy2010/100702/practice/hiki-sample/2body1_sol.f90>))
* ((<2body2.f90|URL:http://itpass.scitec.kobe-u.ac.jp/~itpass/exp/fy2010/100702/practice/hiki-sample/2body2.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 が立ち上がりました. 詳しい使い方は((<こちら|URL:https://itpass.scitec.kobe-u.ac.jp/~itpass/hiki/hiki.cgi?%5BExp2010%5Dgnuplot+%A4%F2%BB%C8%A4%C3%A4%C6%A4%DF%A4%E8%A4%A6>))を見て下さい。では, 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
ここを適当に変えて, 惑星の軌道がどうなるかを確認してください.
本日の実習は以上です. お疲れ様でした.
== 参考資料
* ((<"2009 第11回ITPASS実習資料 実習その１"|URL:https://itpass.scitec.kobe-u.ac.jp/~itpass/hiki/hiki.cgi?%5BExp2009%5D%BC%C2%BA%DD%A4%CB%B7%D7%BB%BB%A4%B7%A4%C6%A4%DF%A4%EB+%28%A3%B2%C2%CE%CC%E4%C2%EA%29">))
* 本実習で用いているプログラム 2body1.f90, 2body1_sol.f90, 2body2.f90 は 2006 - 2008 年度惑星科学実習 (計算機) 用として、今枝佑輔氏が作成されたものを基にして作成されました.
((<"スケジュール表・各回資料 (07/02)"|[Exp2010]スケジュール表・各回資料#07-2F02>))
{{toc}}
= 実習の準備
今回の実習では, プログラムをゼロから作るのではなく, すでに途中まで出来上がっているプログラムを用いて行います. そのためのプログラムのソースファイルはすでに Web 上に置いてあります. なので, まずは情報実験機に適当なディレクトリを作って以下のファイルをダウンロードしてください.
* ((<2body1.f90|URL:http://itpass.scitec.kobe-u.ac.jp/~itpass/exp/fy2010/100702/practice/hiki-sample/2body1.f90>))
* ((<2body1_sol.f90|URL:http://itpass.scitec.kobe-u.ac.jp/~itpass/exp/fy2010/100702/practice/hiki-sample/2body1_sol.f90>))
* ((<2body2.f90|URL:http://itpass.scitec.kobe-u.ac.jp/~itpass/exp/fy2010/100702/practice/hiki-sample/2body2.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 が立ち上がりました. 詳しい使い方は((<こちら|URL:https://itpass.scitec.kobe-u.ac.jp/~itpass/hiki/hiki.cgi?%5BExp2010%5Dgnuplot+%A4%F2%BB%C8%A4%C3%A4%C6%A4%DF%A4%E8%A4%A6>))を見て下さい。では, 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
ここを適当に変えて, 惑星の軌道がどうなるかを確認してください.
本日の実習は以上です. お疲れ様でした.
== 参考資料
* ((<"2009 第11回ITPASS実習資料 実習その１"|URL:https://itpass.scitec.kobe-u.ac.jp/~itpass/hiki/hiki.cgi?%5BExp2009%5D%BC%C2%BA%DD%A4%CB%B7%D7%BB%BB%A4%B7%A4%C6%A4%DF%A4%EB+%28%A3%B2%C2%CE%CC%E4%C2%EA%29">))
* 本実習で用いているプログラム 2body1.f90, 2body1_sol.f90, 2body2.f90 は 2006 - 2008 年度惑星科学実習 (計算機) 用として、今枝佑輔氏が作成されたものを基にして作成されました.
((<"スケジュール表・各回資料 (07/02)"|[Exp2010]スケジュール表・各回資料#07-2F02>))