物理シミュレーションで体験する東大物理

きっかけ

物理シミュレーションゲームを作るには、物理の知識が必要です。
しかし、私の場合は逆に物理シミュレーションゲームを作ることで物理の理解が深まったという経験があります。
どうやら、「知識を深めること」と「物理シミュレーションをすること」は、車輪の両輪であり二人三脚で進歩していくようです。
今回は、物理シミュレーションを通じて物理の理解を深めることに挑戦してみようと思います。

問題

今回は高校物理の力学の問題を物理シミュレーションで再現してみます。
問題としては、東京大学の2003年前期入試の問題を取り扱ってみようと思います。
「いきなり東大！？」と思う方もいるかもしれませんが、大丈夫です。
今回は数式を立てて問題を解くのではなく、問題文の条件を物理シミュレーションに落とし込んで、実際に動かして体験してみることが目的なのです。

さて、肝心の問題文は以下の通りです。

問題文

図1

図 $1$ のように、質量 $2M$ の物体 $A$ と質量 $M$ の物体 $B$ が、ばね定数 $k$ の質量の無視できるばねによってつながれて、なめらかで水平な床の上に静止していた。また、物体 $A$ はかたい壁に接していた。床の上を左向きに進んできた物体 $C$ が、物体 $B$ に完全弾性衝突して、跳ね返された。右向きを正の向きと定めると、衝突直後の物体 $C$ の速度は $+u_1(u_1>0)$ ，物体 $B$ の速度は $-v_1(v_1>0)$ であった。その後、物体 $B$ と物体 $C$ が再び衝突することはなかった。

Ⅰ　まず、衝突前から物体 $A$ が壁から離れるまでの運動を考える。

衝突前の物体 $C$ の速度 $u_0(u_0 < 0)$ を $u_1$ と $v_1$ を用いて表せ。
ばねが最も縮んだときの自然長からの縮み $x(x > 0)$ を求めよ。
衝突してからばねの長さが自然長に戻るまでの時間 $T$ を求めよ。

Ⅱ　ばねの長さが自然長に戻ると、その直後に物体 $A$ が壁から離れた。

やがて、ばねの長さは最大値に達し、そのときの物体 $A$ と物体 $B$ の速度は等しくなった。その速度を求めよ。
ばねの長さが最大値に達したときの自然長からの伸び $y(y > 0)$ を求めよ。
その後ばねが縮んで、長さが再び自然長に戻ったとき、物体 $A$ の速度は最大値 $V$ に達した。 $V$ を求めよ。

Ⅲ　物体 $A$ が壁から離れた後、物体 $B$ と物体 $C$ の間隔は、ばねが伸び縮みを繰り返すたびに広がっていった。このことからわかる $u_1$ と $v_1$ の関係を、不等式で表せ。

シミュレーション

問題の概観

問題文を読むと、物体 $C$ が物体 $B$ に完全弾性衝突してバネが縮むところまでは想像できるのですが、バネが自然長に戻った直後に物体 $A$ が壁から離れた後の挙動があまり想像できません。
しかし、合格ラインを突破するにはⅡの問題までは解いておきたいです。これは困りましたね…
そこで、シミュレーションによって物体の挙動を観察してみようと思います。

シミュレーションの方針

実際に問題を解く際には、力学的エネルギー保存の法則を利用して各変数の関係式を立式することが有効です。
しかし、物理シミュレーションをする場合は、運動方程式 $ma=F$ に着目した方がやりやすいです。

私たちが見ている映像は紙芝居のようなものであり、大体は1秒間に60枚ほどの画像が使われており、それがものすごい速さで切り替わるおかげであたかも滑らかな動画に見えています。
「1秒間に何枚の画像が使われたか」を、fps(frame per second)という単位で表し、例えば1秒間に60枚ほどの画像が使われる場合は60fpsとなります。

物理シミュレーションでは、物体の時点 $t$ における加速度 $a(t)$ を用いて、 $\Delta t$ 秒経過後の速度 $v(t+\Delta t)$ を、

$v(t+\Delta t) \fallingdotseq v(t) + a(t)\Delta t$

と近似します。
加速度 $a(t)$ とは、「 $t$ 秒時点の速度から1秒あたり $a(t)$ だけ速くなる」という意味ですので、 $\Delta t$ 秒経過後には速度が $v(t)$ から $a(t)\Delta t$ だけ加速している、ということを上記の数式は意味しています。
また、60fpsの場合は、 $\Delta t = 1 / 60$ となります。
これは一次近似なのでもちろん誤差がありますが、fpsが十分大きく、かつ、物体の挙動が単純な場合はほとんど問題になりません。

1フレームごとに、

衝突の判定
ばねの計算
各物体の速度の計算
各物体の位置の計算
グラフィックの計算

を繰り返すことでシミュレーションを実行します。

プログラム

今回はPhaserというJavaScriptのゲームエンジンを用いますが、非常に基礎的な機能だけを使うのでp5.jsでもUnityでも何でも大丈夫です。

まず、パラメータを設定します。

  //定数と変数を定義
  const k = 0.0001 //ばね定数
  const M = 1; //質量
  const l_0 = 100; //ばねの自然長
  let l = l_0;  //ばねの長さ
  let x = 0;  //自然長からの伸び（負なら縮み）
  
  //物体を定義
  const matterA = {
    m:2 * M,  //質量
    pos:0,  //x軸方向の座標
    a:0,  //加速度
    v:0,  //速度
  };
  const matterB = {
    m:M,  //質量
    pos:l_0,  //x軸方向の座標
    a:0,  //加速度
    v:0,  //速度
  };
  const matterC = {
    m:0.5 * M,  //質量 不明だが、「跳ね返された」という記述から、物体Bより小さい値のはず。
    pos:l_0 + 300,  //x軸方向の座標
    a:0,  //加速度 最初は等速直線運動をしているので0
    v:-0.5,  //速度 最初は等速直線運動をしているのでいったん負の適当な値でおく
  };

次に、上記で設定したパラメータを基に、物体などを描画します。

  //物体の描画
  const xLeft = 50; //グラフィック用 左側の余白
  const yPos = 500; //y座標 問題には関係なく、グラフィック用
  const r = 10; //球の半径 問題には関係なく、グラフィック用
  const spring = scene.add.sprite(xLeft + (matterA.pos + matterB.pos) / 2, yPos, "spring");
  spring.setScale((matterB.pos - matterA.pos) / spring.width);  //ばねのサイズ調整
  const colorA = 0xe15d60;
  const colorB = 0x00daa6;
  const colorC = 0xf2ad31;
  const colorCenter = 0xf94138;
  const circleA = scene.add.circle(xLeft + matterA.pos, yPos, r, colorA);
  const circleB = scene.add.circle(xLeft + matterB.pos, yPos, r, colorB);
  const circleC = scene.add.circle(xLeft + matterC.pos, yPos, r, colorC);
  //分かりやすさのために、AとBの重心を描画する。
  const centerX = (matterA.pos * matterA.m + matterB.pos * matterB.m) / (matterA.m + matterB.m);
  const circleCenter = scene.add.circle(xLeft + centerX, yPos, r * 0.5, colorCenter);

最後に、毎フレームの更新処理を記述します。
ここでは、updateという関数が毎フレームで呼び出されるものとし、先ほど述べた1から5までの処理をコードに落とし込みます。

  update = (time:number, delta:number) => {
      const deltaT = delta * 0.5;
      //1.衝突の判定
      if(matterB.pos >= matterC.pos){
        //物体BとCが衝突しているとき、完全弾性衝突の式に従って速度を更新
        const newB_v = ((matterB.m - matterC.m) * matterB.v + 2 * matterC.m * matterC.v)
                       / (matterB.m + matterC.m);
        const newC_v = ((matterC.m - matterB.m) * matterC.v + 2 * matterB.m * matterB.v)
                      / (matterB.m + matterC.m);
        matterB.v = newB_v;
        matterC.v = newC_v;
      }

      //2.ばねの計算
      l = matterB.pos - matterA.pos;
      x = l - l_0;
      //運動方程式ma=Fに従い、加速度を更新する。
      //運動方程式ma=Fとフックの法則F=kxより、a=kx/mとなる。
      matterA.a =   k * x / matterA.m;
      matterB.a = - k * x / matterB.m;
      //物体Aが壁に接していて、かつ、壁方向の加速度があるるとき、加速度を0にする
      if(matterA.pos === 0 && matterA.a < 0){
        matterA.a = 0;
      }

      //3.速度の更新
      matterA.v += matterA.a * deltaT;
      matterB.v += matterB.a * deltaT;
      matterC.v += matterC.a * deltaT;

      //4.位置の更新
      matterA.pos += matterA.v * deltaT;
      matterB.pos += matterB.v * deltaT;
      matterC.pos += matterC.v * deltaT; 

      //5.グラフィックの更新
      circleA.x = xLeft + matterA.pos;
      circleB.x = xLeft + matterB.pos;
      circleC.x = xLeft + matterC.pos;
      circleCenter.x = xLeft + (matterA.pos * matterA.m + matterB.pos * matterB.m) 
        / (matterA.m + matterB.m);
      spring.x = xLeft + (matterA.pos + matterB.pos) / 2;
      spring.setScale((matterB.pos - matterA.pos) / spring.width, spring.scaleY);
  }

結果

シミュレーションの結果は以下の通りです。

運動する物体の上にあるグラフは、横軸に時間を、縦軸に速度をとって、物体の時間と速度の関係を表現したものです。
正の速度のとき、すなわち右方向に進むときは、グラフの横軸の上側に点をプロットしています。
また、考察のためにばねに繋がれた物体 $A-B$ を一つの物体と見なしたときの重心を赤色の円で表示しています。

考察

単振動

上記のコードで分かる通り、今回のシミュレーションに用いているのは運動方程式 $ma=F$ 、フックの法則 $F=kx$ 、そして完全弾性衝突の式だけです。
そこにはもちろん三角関数が全く出ていないのですが、物体 $A$ と物体 $B$ の速度のグラフは綺麗なsinカーブを描いております。（これはとても神秘的なことだと思います。）
運動方程式とフックの法則から単振動の式を導出するには大学で習う微分方程式の知識が必要なのですが、微分方程式を解かずともこのように視覚的に体感できるのは物理シミュレーションの利点だと思います。

重心の運動

物体 $A$ が壁を離れてからは、物体 $A-B$ を一つの物体と見なしたときの重心は等速直線運動をします。
そして、物体 $C$ も等速直線運動をしているので、Ⅲの問題は重心の速度と物体 $C$ の速度の大小関係を求める問題に帰着しそうだな、と想像することができます。
上記の場合、物体 $C$ の速度の方が重心の速度よりも大きいことがグラフから分かり、したがって時間が経つにつれて間隔が広がっていきます。

速度の分解

物体 $A$ が壁を離れた後の物体 $A$ および物体 $B$ の速度は、「重心の速度」と「重心からの相対速度」に分解することができます。
このことを数式で表現すると、

$v(t) = v_{center} + A\mathrm{sin}(\omega t + \alpha)$

となります。
ここで、

$v_{center}$ : 重心の速度
$A$ : 振幅
$\omega$ : 振動数
$\alpha$ : 初期位相

としています。
重心は等速直線運動をするので、その速度は時間 $t$ に依存しない定数となります。
$v(t)$ のグラフは、sinカーブを $v_{center}$ だけ上にシフトさせたグラフになり、それは上記のシミュレーションからも視覚的に分かります。

最後に

実際に問題を解く際には力学的エネルギー保存の法則に着目して変数の関係を整理して素早く答えを導き出すのでしょうが、その前段階として「そもそもどういう運動をするのだろうか？」ということを想像するためには物理的直観を養う必要があります。
物理シミュレーションによって物理的直観を養うことができれば、応用的な問題が出たとしても柔軟に対処できるでしょうし、何よりとても面白いと思います。
個人的には、フックの法則を使用した物理シミュレーションは初めてだったのですが、実装の簡単さに対して非常に奥深い動きをして楽しかったので、またやってみたいです。

前の記事
【MBTI】IN型が会社で他タイプとうまくやっていくには

次の記事
藤木雅彦展 ―幾何学的アートの試み― で学んだこと

Feedback

あなたの一言が大きなはげみとなります！