Processingでシミュレーション～放物運動

はじめに

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))$の3成分があります。
$\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階差分を計算する

ここでは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;
}