コイントスの結果を知る

最初の記事だからちょっと面白いことを書こうとがんばった結果がこれだよ....
さて突然だが, コインの裏表やらサイコロの出目やらの話になると, 必ずといっていいほど 「関係するパラメタが多すぎて実際上予測できない」といった文言に出くわす.
ぼくもそうなんだろうなとは思うが, 実際に確認したことはなかったのでやってみた.
以下の計算はつまり, コインの裏表を予測するのは非常に理想化した状況でもこんなに大変なんだなあ, ということを確認するためにある.

問題設定

半径 d, 厚さ 0, 質量 m のコインをトスする. 空気抵抗などは全て無視する.
コインの回転軸は常に紙面に垂直だとする().
コインの進行方向を x 軸正, 鉛直上向きを y 軸正, 角度は xy 平面での反時計回りを正にとる(図1).
コインは何度か床で跳ねた後, 表裏どちらかの面を上にして静止する.
床は完全に水平で表面に凹凸は一切なく, 床とコインの反発係数は e{0\lt e \lt 1} )とする.
コインの重心速度を {\boldsymbol{v} = (v_x, v_y) }, 角速度を \omega とする.

f:id:unv-twinrocket:20170330181301p:plain:w280 f:id:unv-twinrocket:20170330183000p:plain:w185
図1 床に衝突するコイン 図2 コイントスの瞬間

角度の追跡

まず, コインが床と衝突する際のコインの角度 \theta を求める.
あらかじめコインの両端に点 \mathrm{A} , \mathrm{B} を定めておく(コイン上のこの2点のみが床と衝突する).
n 回目( {n \geq 1} )の衝突が時刻 {t=T_n} に起こるとし, その直後のコインの角度, y 方向の速度, 角速度をそれぞれ \theta_n ,  v_{y,n+1} , \omega_{n+1} で表すことにする.
床に衝突したしたほうのコインの端をいったん \mathrm{A} とし, \mathrm{A} , \mathrm{B} の高さを  h_\mathrm{A}(t) ,  h_\mathrm{B}(t) で表すと, \begin{align} h_\mathrm{A}(T_n)=0, \qquad h_\mathrm{B}(T_n)=2d \sin\theta'_n. \qquad \left(\theta':=\theta-\pi\left\lfloor\frac{\theta}{\pi}\right\rfloor\right) \end{align} となる.
ただし {\lfloor x\rfloor} で実数 x の整数部分を表し, 衝突の際の角度を {0 \lt\theta'\lt\pi} におさめるために \theta' を導入した( {0 \lt\theta'\lt\pi} の条件がなければ2通りの衝突を1つの式で表せない).
重力に引かれて v'_y は次第に減少するが, 角運動量保存が成り立つので空中での \omega は一定である.
よって {t\geq T_n} での h は \begin{align} h_\mathrm{A}(t) &= d\sin\theta'_n + v_{y,n+1}(t-T_n)-\frac{1}{2}g(t-T_n)^2 - d\sin(\theta'_n+\omega_{n+1}(t-T_n)), \tag{1}\\ h_\mathrm{B}(t) &= d\sin\theta'_n + v_{y,n+1}(t-T_n)-\frac{1}{2}g(t-T_n)^2 + d\sin(\theta'_n+\omega_{n+1}(t-T_n)).\tag{2} \end{align} {h(t)=0}{t \gt T_n} で解けば, その解 t の最小値が次の衝突時間 {t=T_{n+1}} となる.
すると {n+1} 回目の衝突の際のコインの角度 \theta_{n+1} は \begin{align*} \theta_{n+1} = \theta_n+\omega_{n+1}(T_{n+1}-T_n) \tag{3} \end{align*} によって求まる.

速度と角速度の追跡

次に, 衝突の際の v_y\omega を求める.
コインが床に衝突した瞬間, 床からコインには力積  \boldsymbol{J}_n が働く.
床からコインに働く力は正の y 成分のみを持つので, 力積  \boldsymbol{J}_ny 成分 {J_n \gt 0} のみを持つ.
よって運動量変化と角運動量変化は次の式で表される: \begin{align} m\boldsymbol{v}_{n+1} &= m\boldsymbol{v}'_n+\boldsymbol{J}_n, \\ I\omega_{n+1} &= I\omega_n - J_nd\cos\theta'_n. \end{align} ただし, \boldsymbol{v}'_n は衝突直前のコインの重心速度で, x 方向の運動量保存と, y 方向のエネルギー保存から \begin{align} \boldsymbol{v}'_n = \begin{pmatrix} v_{x,n} \\ -v_{y,n} \end{pmatrix} \end{align} である.
また, I はコインの慣性モーメント (Wikipedia: 慣性モーメント) で, コインが回転軸の周りにどの程度よく回るかを表す定数である(定数であることがわかればそれ以上のことが分からなくても問題ない).
更に, 反発係数の定義より \begin{align*} e = -\frac{(\boldsymbol{v}_{n+1}-d\boldsymbol{\omega}_{n+1})_y}{(\boldsymbol{v}'_n-d\boldsymbol{\omega}_n)_y} = -\frac{v_{y,n+1}-d\omega_{n+1}\cos\theta'_n}{-v_{y,n}-d\omega_n\cos\theta'_n} \end{align*} も成立する(ここで注意すべきは, 反発係数は床と衝突点の相対速度の鉛直成分について考える).
以上をまとめて \begin{align} mv_{x,n+1} &= mv_{x,n}, \\ mv_{y,n+1} &= -mv_{y,n}+J_n, \\ I\omega_{n+1} &= I\omega_n - J_nd\cos\theta'_n, \\ -e &= \frac{v_{y,n+1}-d\omega_{n+1}\cos\theta'_n}{-v_{y,n}-d\omega_n\cos\theta'_n} \end{align} となる. これを解いて \begin{align} v_{x,n+1} &= v_{x,n}, \tag{4}\\ v_{y,n+1} &= -v_{y,n} + (1+e)I\,\frac{v_{y,n}+d\omega_n\cos\theta'_n}{I+d^2m\cos^2\theta'_n}, \tag{5}\\ \omega_{n+1} &= \omega_n - (1+e)dm\,\frac{v_{y,n}+d\omega_n\cos\theta'_n}{I+d^2m\cos^2\theta'_n}\cos\theta'_n. \tag{6} \end{align} これで v_y\omega は漸化式のように順次求まる.

コインの裏表の判定

最後に, コイントスした瞬間と, コインの裏表が決定する瞬間について述べる.
コイントスの瞬間には, 水平に静止したコインに図2のように力積 \boldsymbol{J} が加わる.
この力積が加わる点を \mathrm{A} とするのがわかりやすい.
運動量および角運動量は最初の議論と同様である.
{\boldsymbol{J}=(J_x, J_y)} として \begin{align} mv_{x,0}=J_x, \qquad mv_{y,0}=J_y, \qquad I\omega(=I\omega_1) = -J_yd. \end{align} もちろん \boldsymbol{J} は人間の目には見えないが, コインの高さの最大値 h_M およびそのときのコインと手の(水平)距離 l から逆算できる: \begin{align} v_{y,0}^2 =2g(h_M-h_0), \qquad l = v_{x,0}\,v_{y,0}/g. \end{align} ( h_0 は手の高さ)よってここから {\omega_1=-dm\sqrt{2g(h_M-h_0)}\,/I} とわかり, 1回目の衝突までの  h_\mathrm{A}(t) ,  h_\mathrm{B}(t) が \begin{align} h_\mathrm{A}(t) &= h_0 + v_{y,0} t-\frac{1}{2}gt^2 - d\sin\omega_1 t, \\ h_\mathrm{B}(t) &= h_0 + v_{y,0} t-\frac{1}{2}gt^2 + d\sin\omega_1 t \end{align} とわかる.
なお,  v_{y,1} はエネルギー保存から \begin{align} v_{y,1} = \sqrt{v_{y,0}^2+2gh_0} = \sqrt{2gh_M} \end{align} と求まる.
コインの裏表が決定するのは, コインがひっくり返れなくなったときである.
エネルギーで見ると \begin{align} \frac{1}{2} mv_{y,N}^2 + mgd\sin\theta_N \lt mgd \end{align} となる最小の {N \gt 0} がコインの裏表を決める.
このときの \theta_N を見て, \begin{align} \Theta := \theta_N - 2\pi\left\lfloor\frac{\theta_N}{2\pi}\right\rfloor \end{align} が {\pi/2\lt\Theta\lt3\pi/2} ならばコイントスのときに下だった面が上を向き, そうでないならばトスのときに上だった面が上を向いている.

実際に計算してみた

これを瞬時に予測できる能力があるなら別のことに生かしたほうがよさそうではあるが, 参考までに h(t) を適当なパラメタで計算したものを添付しておく.
トスの瞬間に力積が加わった点(\mathrm{A})が青く描画されている.

f:id:unv-twinrocket:20170331180721p:plain:w195 f:id:unv-twinrocket:20170331180738p:plain:w195 f:id:unv-twinrocket:20170331180755p:plain:w195
図3 コイン両端の高さ 図4 トス直後の拡大 図5 衝突の瞬間の拡大

設定したパラメタ: 重力加速度 {g=9.80} \mathrm{m/s^2}, 反発係数 {e=0.20}.
コインは {d=1.13\times10^{-2}} \mathrm{m}, {m=4.80\times10^{-3}} \mathrm{kg}百円硬貨).
この場合コインの慣性モーメントは {I =(d^2m/4=) 1.53\times10^{-7}} \mathrm{kg\cdot m^2}.
また, 手の高さ {h_0=1.50} \mathrm{m}, 最大高さ {h_M=1.70} \mathrm{m}, 打点の影響などで最初の角速度 \omega_1{|\omega_1|=200} \mathrm{/s}(1秒間に約30回転)と小さく見ておく.
すると {v_{y,0} =1.98} \mathrm{m/s}, {\omega_1=-200} \mathrm{/s}, {v_{y,1}=5.77} \mathrm{m/s}, \begin{align*} &T_1 = 0.790 ~\mathrm{s}, \quad \theta_1 = -157.92, \quad v_{y,2} = -2.63~\mathrm{m/s}, \quad \omega_2 = 541.16~\mathrm{/s}. \\ &T_2 = 0.792 ~\mathrm{s}, \quad \theta_2 = -156.81, \quad v_{y,3} = 3.46~\mathrm{m/s}, \quad \omega_3 = 258.14~\mathrm{/s}, \\ & \vdots \\ &T_8 = 1.962 ~\mathrm{s}, \quad \theta_8 = -11.72, \quad v_{y,9} = 0.57~\mathrm{m/s}, \quad \omega_9 = 73.02~\mathrm{/s}, \\ &T_9 = 2.085 ~\mathrm{s}, \quad \theta_9 = -2.72, \quad v_{y,10} = 0.21~\mathrm{m/s}, \quad \omega_{10} = -45.46~\mathrm{/s}, \end{align*} となり, 9回目の接触によって裏が出る(トスのとき下を向いていた面が上になる)とわかる.

注について

実際問題として, コイントスが今回の設定のように上手くいくことはありえないであろう.
コインの回転軸と飛ぶ方向は, 打点(コインに力積が加わる点)と力積の向きによって決定される.
これを考慮に入れた場合も少し考えてみたのだが, これは難しくしかもかなり長くなりそうである.

参考文献

[1] http://www.th.phys.titech.ac.jp/~muto/lectures/Gmech08/chap12.pdf, 第12章 剛体の運動(2), 武藤研究室.
また, 図2は http://icooon-mono.com/ から, 図3,4,5は http://www.osaka-kyoiku.ac.jp/\~tomodak/grapes/ を用いて作成した.