소개

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|$입니다. 두 물체에 작용하는 만유 인력은 다음 식으로 주어진다.

물체 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();
}