惑星の運動のシミュレーション

惑星の運動のシミュレーション

先ほどの丸が移動してゆくプログラムをもとに、太陽系を動く惑星のような、重力によって速度が変わる系のシミュレーションをしてみましょう。

さきほどのプログラムを一部変更し、○の位置を表す変数(x, y)に加え、○の速度を表す変数(vx, vy)を定義します。

x = 0.0
y = 0.0
vx = 1.0
vy = 1.0

def setup():
    size(400, 400)
    smooth()  #円の縁を滑らかにする
    stroke(255)
    background(0)
def draw():
    global x, y, vx, vy
    ellipse(x, y, 5, 5)
    x = x + vx
    y = y + vy

このままだと、全く前と同じく、○が右下に向かって移動していきます。(等速直線運動)
これから、画面の中心(200,200)に太陽のような重力中心があり、○が重力によってその周りを回るような運動のシミュレーションを考えます。

vxやvyは○の速度を表していますが、その単位は m/s などではなく、
pixel/frame となっていることに注意してください。
pixelは画面上の画素を単位とした長さ、
frameはアニメーションの一コマ分(draw1回分)の時間です。

さて、○を宇宙空間を漂う物体とみなして、これに力が加わる時の運動について考えてみましょう。
物体の運動と力は、次の運動方程式で記述されます。

 F = ma

a は加速度で、単位時間当たりの速度の変化を表しています。
また、2物質間に働く重力は

 \displaystyle F = -G\frac{M m}{r^2}

以上2式から、

 \displaystyle a = -G\frac{M}{r^2}

であることがわかります。

r2 は べき乗を計算する命令 pow(a,b)を使い、
(pixel単位)で次のように求められます。
計算結果は rs という変数に入れておきましょう。
下の命令は、関数drawの中に入れます。

rs = pow(x-200,2)+pow(y-200,2);

これで物体への加速度 a を求めることができます。
G と M を決めておかなくてはなりませんが、これは GM = 100 としておいてください。

a = -100/rs;

a はベクトルなので、x 方向と y 方向に分けておかなくてはなりません。
例えば a の x成分(ax) は a に r と x の比(x/r)をかけることで求めることができます。
(ここでの x はx軸方向の○と中心の間の距離)
r は sqrt(rs),
x は (x-200)
なので、ax は
ax = a * x / r;
ですね。符号が中心に向かう方向になるように 注意しなくてはなりません。

あとは次のように、draw中に 加速度を加える処理をつければ惑星の運動がシミュレートできます。
x, y の初期値を (200, 100)、
vx, vyの初期値を (1, 0) として円運動をすれば正解です。

def draw():
  ellipse(x, y, 5, 5)

  #ここでax, ayを求める

  vx = vx + ax
  vy = vy + ay

  x = x + vx
  y = y + vy

演習

  • 上のプログラムの #ここでax, ay を求める のところに正しい命令を入れて、○を円運動させてください。
  • 初期速度を少し変えて、惑星が楕円運動をすることなどを確かめてください。

サンプル