数値解析法及び演習第十回

[1] 運動方程式

質点mが力Fを受けてx軸上を運動するときの運動方程式は以下のようにかける．

$\begin{displaymath} m{d^2x\over dt^2}=F \end{displaymath}$

(1)

この二階常微分方程式は次の連立微分方程式に変換することができる．

$\begin{displaymath} m{dv\over dt}=F \end{displaymath}$

(2)

$\begin{displaymath} {dx\over dt}=v \end{displaymath}$

(3)

練習問題1

先週の課題を参考に，次の運動方程式（太陽からの重力を受ける惑星の運動：二次元の二階常微分方程式）をルンゲクッタ法で解くプログラムを作成せよ．以下の式では，GMm=1ととってある．

$\displaystyle {d^2x\over dt^2}$	$\textstyle =$	$\displaystyle -{x\over (x^2+y^2)^{3/2}}$
$\displaystyle {d^2y\over dt^2}$	$\textstyle =$	$\displaystyle -{y\over (x^2+y^2)^{3/2}}$	(4)

変数がx,y,vx,vyと四つあることに注意せよ．t=0でx=1, y=0, vx=0, vy=1とすること．

	
	write(1,*) t, x, y

と，時刻，x座標，y座標と並べてファイルに出力する．gnuplotで以下のコマンドを実行することでアニメーションを観察すること．

gnuplot> n=1
gnuplot> load 'anim.txt'

ここで，anim.txtの中身は以下のとおり．

  set xrange [-1:1]
  set yrange [-1:1]

  if (n>ループ回数) exit

  plot "ファイル名" using 2:3 every 1:1:n:0:n:0 with points pt 7 ps 2,\
 "zero" with points pt 6 ps 3 

  pause 0.05

  n=n+1

  reread

"ファイル名"に結果を出力したファイル名を記入すること．時間ループの回数を"ループ回数"に記入すること．また，ファイル"zero"は原点に存在する太陽の位置を表し，内容は

0.D0 0.D0

としておくこと．

[2] ケプラーの法則

ケプラーの法則は、1619年にヨハネス・ケプラーによって解明された惑星の運動に関する法則である。ケプラーは、ティコ・ブラーエの観測記録から、太陽に対する火星の運動を推定し、以下のように定式化した。

第1法則 : 惑星は太陽をひとつの焦点とする楕円軌道上を動く。
第2法則 : 惑星と太陽とを結ぶ線分の描く面積は単位時間あたり常に一定である（面積速度一定）。
第3法則 : 惑星の公転周期の2乗は軌道の半長径の3乗に比例する。

[3] 第二法則を発見しよう

練習問題2

nステップ目の座標(x_n,y_n)，n+1ステップ目の座標(x_n+1,y_n+1)，および原点(0,0)からなる三角形の面積を求めよ．時間とともにどう変化するかgnuplotを用いてグラフを作成せよ．x, yを増やす前に，xとyの値をxold, yold と別の変数に代入しておくと良い．

[4] 第三法則を発見しよう

練習問題3

惑星の初期x座標を5とおり変化させて計算を行う．惑星が太陽の周りを一周するのにかかる時間（周期），および軌道長半径（（初期ｘ座標-軌道を半周したときのx座標）/2）を求めよ．横軸に周期，縦軸に軌道長半径をとりグラフを作成せよ．軌道長半径は周期の何乗に比例しているか？軌道を半周したときのx座標を求めるためには，前述したyoldをつかって，y*yold < 0.D0 となる時のx座標を求めればよい．

[5] 第一法則を発見しよう

練習問題4

惑星の初期ｘ座標を１とせよ．初期速度v_yを変化させよ．どの速度で惑星は無限遠方に飛んでいってしまうか？　またそれはなぜか？

練習問題5

練習問題4で求めた速度よりも小さい初期速度v_yで計算を行え．ケプラーの第一法則により軌道は楕円になっているはずである．

二つの点（焦点）からの距離の和が一定になる点の集合が楕円である．今の場合，一方の焦点は原点（太陽），もう一方の焦点のｘ座標は（初期ｘ座標＋軌道を半周したときのx座標）で与えられる．一度軌道計算を行うと，焦点のｘ座標を求めることができる．まず焦点のｘ座標を求め，同じ条件で再度計算を行い，二つの焦点からの距離の和が一定になっていることを確かめよ．

[6] 今日の宿題 (1.5点：12/22日まで)

プログラムと，練習問題5の結果（距離の和の一定値）を送付すること．

[7] 追加の宿題 (0.5点：12/22日まで) 先を走る宇宙船にたどりつこう

半径が１，速度が１の円軌道上を別の宇宙船が運動している．t=0でその宇宙船はx=-1,y=0の位置にいる．一方で自分はt=0でx=1,y=0の位置にいる．なるべく早くその宇宙船にたどりつくためにはy方向の初速度vyをどのような値にすればよいだろうか？(vx=0とする)

先に走る宇宙船の座標は(piを円周率として）

         write(2,*) t,cos(t+pi),sin(t+pi)

でファイルに出力することができる．このファイル名を"data2.dat"，自分の座標が入ったファイルを"data1.dat"とすると，次のファイルでアニメーションを見ることができる．追い付いたかどうかをアニメを見て確認すること．

  set xrange [-1:1]
  set yrange [-1:1]

  if (n>700) exit

  plot "data1.dat" using 2:3 every 1:1:n:0:n:0 with points pt 7 ps 2,"zero" with points pt 6 ps 3,"data2.dat" u 2:3 every 1:1:n:0:n:0  

  pause 0.05

  n=n+1

  reread

ヒント：今の問題で，軌道長半径aとエネルギーの絶対値との間には以下の関係がある．

$\begin{displaymath} a={1\over \vert E\vert}=\left\vert-{1\over r}+{1\over 2}v^2\right\vert^{-1} \end{displaymath}$

(5)

ただしここでrとvはそれぞれ以下で与えられる．

$\displaystyle r$	$\textstyle =$	$\displaystyle \sqrt{x^2+y^2}$
$\displaystyle v^2$	$\textstyle =$	$\displaystyle \left({dx\over dt}\right)^2+\left({dy\over dt}\right)^2$	(6)

一方で，第3法則から，惑星の公転周期の2乗は軌道の半長径の3乗に比例する。

日程表へ戻る <<

数値解析法及び演習 第十回