ボールの軌跡

 投げたボールの軌跡を計算しグラフィックスで表示をするのが目的です。空気抵抗がない場合、ボールの軌跡は放物線(時間に対する2次式)になります(これが「放物線」の名前の意味です)。ここでは、空気抵抗を考慮したボールに軌跡の表示を試みます。

  1. 運動方程式

     vx を x 方向の速度、dx を時間 dt の間に進む x 方向の距離とすると、速度は時間当たり移動する距離ですから次の関係があります。
         dx=vx * dt ;速度の定義
     物体に力を加えると、速度が変化します。時間当たりの速度の変化率を加速度といいます。加速度と力の間には次の関係があります。
         質量 * 加速度 = 力  :ニュートンの法則
    これが、力と質量、加速度の間の関係式になります。質量 m の物体に力Fを加えると、加速度aは次のようになります。
      m * a = F :力と加速度
    投げたボールに働く力は重力のみです。重力はボールに y(垂直)方向の加速度 g = 9.8 を与えます。

  2. 空中のボールの動き

    ボールの動きは、水平方向と垂直方向に分けて別々に計算することができます。y方向の速度をvy、x(水平)方向の速度をvx、とすると、dt時間で移動する x、y方向の距離 dx、dyは
           dy = vy * dt, dx = vx * dt 
    となります。


    手から離れた直後のボールの x,y 方向の速度を vx0 と vy0 とします。これを初速度といいます。
     手から離れた直後のボールの速度をv、投げた角度を ang とすると、
            vx0 = v*cos(ang)、vy0 = v*sin(ang)
    となります。cos(ang)、sin(ang)は三角関数で、vのx成分とy成分を計算します。
     空中のボールには引力に基づく重力のみが作用します。したがって、水平方向の速度は変化しません。一方、垂直方向には下方向に重力 F = m * g が作用します。g は実験で求められている数値(9.8m/ss)を利用します。投げ上げたボールは引力により下向きの加速度gを受け、秒当たりgだけ速度が減っていきます。したがって、t 秒後の速度変化は
            vy(t)=vy0 - g * dt
    となります。vy0/g 秒後にはy方向の速度は0になり、以後落下をはじめます。

       y方向の速度の変化

     ボールが最高点に達する時刻tは、ボールのy方向の速度が0になる時刻ですから、
     vy0 - g * t = 0  したがって、
           t=vy0/g=v・sin(ang)/g 
    となります。

  3. 軌跡の計算

     時刻tに対する式を求めるのでなく、計算機で数値計算を行う方法をシミュレーションといいます。ここでは、数値計算の基礎的な手法を紹介します。
     veを初速度、vx0、vy0をX,Y方向の初速度とすると。
         vx0=ve*(float)cos((3.14/180.0)*ang);
         vy0=ve*(float)sin((3.14/180.0)*ang);
    となります。ここで、(3.14/180.0)*ang で角度をラジアン単位に直しています(180度がπラジアンです)。
     次に、重力の加速度をg、計算する時間間隔dtを設定します。すると、以下の関係式が成立します。
       x = x + vx * dt;//位置変化
       y = y + vy * dt;
       vy = vy - g * dt;//速度変化
     ここで、dtは十分短くする必要があります。ここでの、計算は「dtの間はvx,vyは変化しない」、ことを前提としています。dtを大きくすると、この条件が成立しなくなるからです、
     あとは、繰り返しを用いて、y>0 の間、dt時間後の位置x,yを求め、y方向の速度を更新します(この数値計算方法はオイラー法とも呼ばれます)。
    {
    // 軌跡の計算
     float ve,ang,vx0,vy0;
     float g,dt,t1,t,x,y,vx,vy; 
     //速度設定、角度読み込み、初速度計算
     ve=30.0;
     vx0=ve*(float)cos((3.14/180.0)*ang);
     vy0=ve*(float)sin((3.14/180.0)*ang);
    
     //重力加速度、刻み時間、落下時刻
     g=9.8f;
     dt=0.01f;
     //初期設定
     x=0;y=0;
     vx=vx0;vy=vy0;
    
     for( t=0.0;y>=0;t=t+dt){
      x = x + vx * dt;//位置変化
      y = y + vy * dt;
      vy = vy - g * dt;//速度変化
     }
    }

  4. 空気抵抗がある場合の計算法

     一般には投げ上げた物体は空気の抵抗を受けて減速します。空気の抵抗は、形状が一定で回転を考慮しないと、速度の2乗に比例します。v を速さ、k を空気抵抗の定数とすると空気抵抗による加速度はいかのようになります。
        -k*v*v
     空気抵抗があると、vx,vyとも変化します。v は 速度成分 vx,vy から、三平方の定理で合成します。合成した加速度 v から 各加速度の成分を求めるには k*v*v*cos A を計算する必要がありますが、v*cos A = vx ですから、加速度の x 成分は k*v*vx で計算できます。
     したがって、空気抵抗がある場合の dt 時間後の速度は
           v=sqrt(vx*vx+vy*vy);
           vx=vx - k*v*vx*dt
           vy=vy - g *dt - k*v*vy*dt
    となります。

  5. 軌跡の計算プログラム

     setup()で定数を設定、draw()で軌跡を描画します。開始点は(10,10)とします(描画するy座標は、height-10 になります)。時間は0.1秒間隔、速度の初期値は60m/s、空気抵抗は0とします。
     フレーム毎に、経過した位置を円で描画し、y の値が10以下になったら、noloop() で描画を停止します。キーイベントで、初期速度、角度、空気抵抗定数、が変更可能です。空気抵抗の値はゴルフボールの値を参考にしましたが、実際より少し大きい値です(実際には 0.00002程度)。
     ディスプレイウインドウの最下行に各パラメータの値を描画しています。
    // 軌跡の計算
    int ve,ang;
    float vx0,vy0;
    float g,dt,t1,t,x,y,vx,vy; 
    float k;
    //速度設定、角度読み込み、初速度計算
    
    
    void setup(){
      ve=60;
      ang=45;
      init(ve,ang);
      g=9.8f;
      dt=0.1f;
      k=0.0;
      size(500,200);
    }
    
    void draw(){
      float v;
      v=sqrt(vx*vx+vy*vy);
      x = x + vx * dt;//位置変化
      y = y + vy * dt;
      vy = vy - g * dt-k * v * vy * dt;//速度変化
      vx=vx - k * v * vx * dt;
      ellipse(x,height-y,5,5);
      if(y<10) noLoop();
      fill(150,150,150);
      rect(0,height-10,width,height);
      fill(250,250,250);
      text("v="+ve+" a="+ang+" k="+k,10,height);
    }
    
    void init(int ve,int ang){
      vx0=ve*(float)cos((3.14/180.0)*ang);
      vy0=ve*(float)sin((3.14/180.0)*ang);
      //初期設定
      x=10;
      y=10;
      vx=vx0;
      vy=vy0;
    }
    
    void keyPressed(){
      if(key=='v') ve=ve+3;
      if(key=='V') if(ve>40) ve=ve-3;
      if(key=='a') ang=ang+3;
      if(key=='A') if(ang>0) ang=ang-3;
      if(key=='k') k=k+0.0001;
      if(key=='K') if(k>0.0001) k=k-0.0001;
      if(key=='c') background(200);
      init(ve,ang);
      loop(); 
    }
     実行画面(v=60,a=39,kを変化))


  6. まとめ

     空気抵抗を含む軌跡の理論式はかなり高度な数学を必要としますが、数値計算の場合は計算モデルさえわかれば、比較的簡単に計算できます。オイラー法による時間進行に伴う計算方法は、計算時間が長くなると誤差が累積するため、不正確になります。