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

[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とすること．座標の変化分（つまり速度）を与える関数をdxdt, dydt，速度の変化分（加速度）を与える関数をdvxdt, dvydt とすること．なので関数f, gと増加分のf1, g1などは概ね次のようになる．dx1は先週のf1，dvx1は先週のg1に対応する．

    double precision function dxdt(t, x, y, vx, vy)
    implicit none
    .....


    double pcecision function dvxdt(t, x, y, vx, vy)
    implicit none
    .....


    dx1 = dt * dxdt( t, x, y, vx, vy )
    dvx1 = dt * dvxdt( t, x, y, vx, vy )
    dy1 = dt * dydt( t, x, y, vx, vy )
    dvy1 = dt * dvydt( t, x, y, vx, vy )

    .....

    x = x + (dx1 + 2*dx2 + 2*dx3 + dx4)/6.
    .....

初期条件は t=0.D0 で x=1.D0, y=0.D0, vx=0.D0, vy=1.D0とすること．

	
	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,\
 "sun.dat" with points pt 6 ps 3 

  pause 0.05

  n=n+1

  reread

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

0.D0 0.D0

としておくこと．xrange, yrangeで指定されている値を変化させると，画面で表示される座標の範囲が変わる．適宜変更すること．

[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座標を3とおり変化させて計算を行う(0.25D0, 0.5D0, 1.D0)．惑星が太陽の周りを一周するのにかかる時間（周期），および軌道長半径（（初期ｘ座標-軌道を半周したときのx座標）/2）を求めよ（位置関係は下図参照）．横軸に周期，縦軸に軌道長半径をとりグラフを作成せよ．軌道長半径は周期の何乗に比例しているか？軌道を半周したときのx座標を求めるためには，前述したyoldをつかって，y*yold < 0.D0 となる時のx座標を求めればよい．

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

練習問題4

練習問題3で用いた，初期x座標=0.5D0で計算を行え．ケプラーの第一法則により軌道は楕円になっているはずである．

二つの点（焦点）からの距離の和が一定になる点の集合が楕円である．今の場合，一方の焦点は原点（0.D0, 0.D0），もう一方の焦点のy座標は0.D0, ｘ座標は（初期ｘ座標＋軌道を半周したときのx座標）で与えられる（上図参照）．軌道を半周したときのx座標は練習問題３で求めているはずである．二つの焦点から惑星までの距離の和が一定になっていることを確かめよ．

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

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

[7] 追加の宿題 (1点：12/20日まで) 散乱角度：斥力の場合との比較

dt=1.D-2, 初期位置をx=20.D0, y=0.2D0，初期速度をvx=-5.D0, vy=0.D0とする．

この条件でシミュレーションすると，質点は右から原点に近づき，原点近くで引力により曲げられ，原点から遠ざかる．まったく同じシミュレーションを，力の符号をマイナスからプラスに変えて行えば，斥力が働く場合のシミュレーションとなる（同一電荷を持った粒子が接近する場合）．これらふた通りのシミュレーションを行い，結果のアニメーションを以下のファイルで観察すること．

  set xrange [-20:20]
  set yrange [-20:20]

  
  if (n>700) 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,"斥力でのファイル名" u 2:3 every 1:1:n:0:n:0 w p pt 7 ps 2

  pause 0.01

  n=n+1

  reread

ファイル名はそれぞれ，時間t：x座標：y座標が記録されたファイル名に書き換えること．原点から遠ざかっていく速度と．x軸とがなす角度を引力と斥力の場合でそれぞれ求める．角度は，datan(vy/vx)で求めることができる（これで出る角度は度でなくラジアン)．

プログラムと，引力の場合と斥力の場合の角度を送付すること．

日程表へ戻る <<

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