소개

Processing을 사용하여 포물 운동을 시뮬레이션합니다.

미분에서 차분으로

물리 법칙의 대부분은 미분 방정식으로 설명됩니다.
미분은 무한소를 다루기 때문에 컴퓨터에서는 다룰 수 없습니다.
무한소는 포기하고 허용할 수 있는 작은 변화를 사용하여 미분을 차이로 바꿉니다.

미분 정의

\frac{d}{dt}r(t) = \lim_{Δt→0}\frac{r(t)-r(t-Δt)}{Δt}

1층 차분
lim을 잃었을 뿐입니다. Δt가 충분히 작으면 이것으로 충분하다고 가정합니다.

\frac{d}{dt}r(t) \sim \frac{r(t)-r(t-Δt)}{Δt} = v(t)

2층 차분
$v(t)$에 대해 1층 차분을 더 계산합니다.

\frac{d^2}{dt^2}r(t) \sim \frac{d}{dt}v(t) \\
\sim \frac{v(t)-v(t-Δt)}{Δt} \\
 = \frac{r(t)-r(t-Δt)-(r(t-Δt)-r(t-2Δt))}{Δt^2} \\
 = \frac{r(t)-2r(t-Δt)+r(t-2Δt)}{Δt^2} 

이산화

$r(t)$와 $t$는 연속량이지만 컴퓨터에서는 취급할 수 없습니다.
연속량을 이산량으로 대체하여 생각합니다.

t는 0,1,2,3, ...로 가정합니다. Δt는 1이 됩니다.
r(t)는 r(0),r(1),r(2), ...입니다.

운동 방정식

물체의 위치 $\vec{r}$ 는 시간에 따라 변화하므로 시간의 함수가 됩니다.
$\vec{r} = (r_x(t), r_y(t), r_z(t))$의 세 가지 구성 요소가 있습니다.
$\vec F$는 힘, $\vec a$는 가속도, $\vec v$는 속도, $\vec r$는 위치, $m$는 질량입니다.

벡터 표기

\vec{F} = m\vec{a} = m\frac{d}{dt}\vec{v} = m\frac{d^2}{dt^2}\vec{r}

성분 표기

F_x = ma_x = m\frac{d}{dt}v_x = m\frac{d^2}{dt^2}r_x \\
F_y = ma_y = m\frac{d}{dt}v_y = m\frac{d^2}{dt^2}r_y \\
F_z = ma_z = m\frac{d}{dt}v_z = m\frac{d^2}{dt^2}r_z \\

중력

수직 하향으로 일정한 중력이 걸리면 다음과 같이 됩니다.
수직 하향을 y축, 수평 방향을 x축으로 합니다. z축은 생략합니다. m, g는 정수로 합니다.

\vec{F} = (F_x, F_y) = (0, mg)

운동방정식은

\frac{d^2}{dt^2}r_x(t) = 0 \\
\frac{d^2}{dt^2}r_y(t) = g \\

됩니다.

이 식에서 $r_x,r_y$를 구합니다. 다음과 같은 방법이 가능합니다.

두 번 적분

2 층 차분 계산

여기에서는 2층 차분을 사용하여 근사값을 구합니다.

운동 방정식의 차분화·이산화

차분화

\frac{1}{Δt^2}(r_x(t)-2r_x(t-Δt)+r_x(t-2Δt)) = 0 \\
\frac{1}{Δt^2}(r_y(t)-2r_y(t-Δt)+r_y(t-2Δt)) = g \\

이산화
Δt=1, t=0,1,2,3...로 합니다.
t=0일 때

r_x(0)-2r_x(-1)+r_x(-2) = 0 \\
r_y(0)-2r_y(-1)+r_y(-2) = g \\

입니다. $r_x(-1),r_x(-2),r_y(-1),r_y(-2)$가 나옵니다.
이것은 초기 조건입니다. 첫 번째 위치, 속도에 따라 결정해야합니다.
우선 0으로 합니다.

r_x(0) = 0 \\
r_y(0) = g \\

t=1일 때

r_x(1) = 0 \\
r_y(1) = 2r_y(0)-r_y(-1)+g = 3g \\

t=2일 때

r_x(2) = 0 \\
r_y(2) = 2r_x(1)-r_x(0)+g  = 6g \\

t=3일 때

r_x(3) = 0 \\
r_y(3) = 2r_x(2)-r_x(1)+g  = 10g \\

이렇게 잇달아 시간 발전을 쫓아갈 수 있습니다.
g 앞의 숫자의 줄은 1,3,6,10,15,.

Processing으로 시각화

현재의 위치를 ​​rx2, ry2로 합니다.
rx1,ry1,rx0,ry0은 각각 1시간 전, 2시간 전의 xy 좌표입니다.
프레임 업데이트가 단위 시간 업데이트에 해당합니다.
프레임 속도를 줄이면 시간이 느려집니다.
초기위치(rx0, ry0), 초속도(vx0,vy0), 중력가속도(g)는 조정값이 됩니다.

float g = 1;
float vx0 = 10;
float vy0 = -30;
float rx0 = 0;
float rx1 = vx0 + rx0;
float ry0 = 480;
float ry1 = vy0 + ry0;
float rx2 = 2*rx1 - rx0;
float ry2 = 2*ry1 - ry0 + g;

void setup() {
  size(640,480);
  frameRate(30);
}

void draw() {
  ellipse(rx2, ry2, 10, 10);
  updateTime();
}

void updateTime()
{
  rx0 = rx1;
  ry0 = ry1;
  rx1 = rx2;
  ry1 = ry2;
  rx2 = 2*rx1 - rx0;
  ry2 = 2*ry1 - ry0 + g;
}