はじめに

Processingを使って万有引力をシミュレーションします。

万有引力

質量を持つ物体の間に働く力です。
物体1の位置ベクトルを$\vec{r_1}$、質量を$m_1$とします。物体2の位置ベクトルを$\vec{r_2}$、質量を$m_2$とします。
$G$は万有引力定数です。物体1が受ける力の方向は$\vec{r_2}-\vec{r_1}$です。物体2が受ける力の方向は$\vec{r_1}-\vec{r_2}$です。
$r = |\vec{r}_1 - \vec{r}_2|$です。2つの物体に働く万有引力は次の式で与えられます。

物体1が物体2から受ける万有引力

\vec{F} = G \frac{m_1 m_2}{r^2} \vec{e}_{12} \\
\vec{e}_{12} = \frac{\vec{r}_2 - \vec{r}_1}{r}

物体2が物体1から受ける万有引力

\vec{F} = G \frac{m_1 m_2}{r^2} \vec{e}_{21} \\
\vec{e}_{21} = \frac{\vec{r}_1 - \vec{r}_2}{r}

運動方程式

万有引力のみで外力のない運動方程式は以下の通りです。
質量$m_1$の物体1が$m_2$の物体2から万有引力を受けているとします。
物体1のみ記載します。物体2もほとんど同じため省略します。

\frac{d^2 \vec{r_1}}{dt^2} = G \frac{m_2}{|\vec{r_1}-\vec{r_2}|^2} \vec{e}_{12} \\

差分化

差分化
$r_{1x}$は物体1のx座標です。その他の物体、座標も同様の記号とします。

r_{1x}(t+2Δt)-2r_{1x}(t+Δt)+r_{1x}(t) = Gm_2 r^{-2} e_{x12} \\
r_{1y}(t+2Δt)-2r_{1y}(t+Δt)+r_{1y}(t) = Gm_2 r^{-2} e_{y12} \\
r_{2x}(t+2Δt)-2r_{2x}(t+Δt)+r_{2x}(t) = Gm_1 r^{-2} e_{x21} \\
r_{2y}(t+2Δt)-2r_{2y}(t+Δt)+r_{2y}(t) = Gm_1 r^{-2} e_{y21} \\
r = \sqrt{(r_{1x}(t+2Δt)-r_{2x}(t+2Δt))^2+(r_{1y}(t+2Δt)-r_{2y}(t+2Δt))^2}\\
e_{x12} = \frac{r_{2x}-r_{1x}}{r} \\
e_{y12} = \frac{r_{2y}-r_{1y}}{r} \\
e_{x21} = \frac{r_{1x}-r_{2x}}{r} \\
e_{y21} = \frac{r_{1y}-r_{2y}}{r} \\

Processingで可視化

float G = 100.0;
float m1 = 100.0;
float m2 = 100.0;
PVector[] r1; // r1[time]
PVector[] r2;

void setup() {
  size(800, 500);
  init();
  frameRate(30);  
  fill(255, 255, 255, 255);
  noStroke();
  background(0, 0, 0);
}

void draw() {
  fill(0, 20);
  rect(0, 0, width, height);
  fill(255);
  ellipse(r1[2].x, r1[2].y, 10, 10);
  ellipse(r2[2].x, r2[2].y, 10, 10);
  for (int i=0; i<4; i++) // speed up
    updateTime();
}

void init()
{
  r1 = new PVector[3];
  r2 = new PVector[3];
  float left_mgn = 100;
  float right_mgn = 100;
  float centor = height/2;

  r1[2] = new PVector(left_mgn, centor);
  r2[2] = new PVector(width-right_mgn, centor);  
  r1[1] = new PVector(r1[2].x, r1[2].y+1);
  r2[1] = new PVector(r2[2].x, r2[2].y-1);  
  r1[0] = new PVector(r1[1].x, r1[1].y+2);
  r2[0] = new PVector(r2[1].x, r2[1].y-2);
}

void updateTime()
{
  float r = r1[1].dist(r2[1]);
  float rr = pow(r,2);
  PVector e12 = PVector.mult(PVector.sub(r2[1], r1[1]), 1/r);
  PVector e21 = PVector.mult(PVector.sub(r1[1], r2[1]), 1/r);
  r1[2].x = 2*r1[1].x - r1[0].x + G*m2/rr * e12.x;
  r1[2].y = 2*r1[1].y - r1[0].y + G*m2/rr * e12.y;
  r2[2].x = 2*r2[1].x - r2[0].x + G*m1/rr * e21.x;
  r2[2].y = 2*r2[1].y - r2[0].y + G*m1/rr * e21.y;

  r1[0] = r1[1].copy();
  r2[0] = r2[1].copy();
  r1[1] = r1[2].copy();
  r2[1] = r2[2].copy();
}