野球における打球の軌跡

課題

コード

import numpy as np
import math
import matplotlib.pylab as plt

# 加速度の関数を定義
def a_x(v_x, v_y):
    return - k * math.sqrt(v_x * v_x + v_y * v_y) * v_x
def a_y(v_x, v_y):
    return - g - k * math.sqrt(v_x * v_x + v_y * v_y) * v_y

def sim(degree, v0, dt, n_steps, g, k):
    
    # ラジアンに
    radian = degree * math.pi / 180.0

    # 各方向の初速度と距離を定義
    v_x0 = v0 * math.cos(radian)
    # print(f"v_x0 {v_x0:.10f}")
    v_y0 = v0 * math.sin(radian)
    # print(f"v_y0 {v_y0:.10f}")
    d_x0 = 0
    d_y0 = 0

    # 刻み幅とステップ数を定義
    dt = 0.01
    n_steps = 1000

    # 速度と距離の配列を定義
    v_x = np.zeros(n_steps)
    v_y = np.zeros(n_steps)
    d_x = np.zeros(n_steps)
    d_y = np.zeros(n_steps)

    # 初期値を設定
    v_x[0] = v_x0
    v_y[0] = v_y0
    d_x[0] = d_x0
    d_y[0] = d_y0

    # 2次のRunge-Kutta法を使って微分方程式を解く
    for i in range(1, n_steps):
    
        if d_y[i - 1] >= 0.0:
        
           # dv_x/dt = a_x(v_x, v_y)
           # dv_y/dt = a_y(v_x, v_y)
           k1_vx = a_x(v_x[i-1], v_y[i-1]) * dt
           k1_vy = a_y(v_x[i-1], v_y[i-1]) * dt
           k2_vx = a_x(v_x[i-1] + k1_vx/2, v_y[i-1] + k1_vy/2) * dt
           k2_vy = a_y(v_x[i-1] + k1_vx/2, v_y[i-1] + k1_vy/2) * dt
           v_x[i] = v_x[i-1] + k2_vx
           v_y[i] = v_y[i-1] + k2_vy
        
           d_x[i] = d_x[i-1] + v_x[i-1] * dt
           d_y[i] = d_y[i-1] + v_y[i-1] * dt
        
           # print(f"d_x[i] d_y[i] {d_x[i]:.10f} {d_y[i]:.10f}")

        # 地面に埋もれてるので直前の値を設定
        else:
           v_x = np.delete(v_x, np.s_[i - 1:])
           v_y = np.delete(v_y, np.s_[i - 1:])
           d_x = np.delete(d_x, np.s_[i - 1:])
           d_y = np.delete(d_y, np.s_[i - 1:])
        
           break
        
    # print(d_y)

    # 飛距離と速さを計算
    v = np.sqrt(v_x[-1]**2 + v_y[-1]**2)
    d = d_x[-1]

    plt.plot(d_x, d_y)  # x座標とy座標の組
    plt.xlim(0, 140) 
    plt.ylim(0, 80) 
    plt.show()

    # 結果を表示    
    return v, d, radian

degreeList = []
distanceList = []
if __name__ == '__main__':
    
   # 初速度(m/s)
   v0 = 55
    
   # 重力加速度
   g = 9.8
    
   #  空気抵抗係数
   k = 0.008
    
   # 刻み幅とステップ数を定義
   dt = 0.01
   n_steps = 1000

   for degree in range(10, 91, 5):
       
       v, d, radian = sim(degree, v0, dt, n_steps, g, k)
    
       print(f"角度 {degree:.0f} ラジアン {radian:.5f}")
       print(f"最終速さ {v:.2f} m/s")
       print(f"飛距離 {d:.2f} m")
       print()
       
       degreeList.append(degree)
       distanceList.append(d)
    
   plt.plot(degreeList, distanceList, marker='.',markersize=10)
   plt.xlabel("角度(度)", fontname="MS Gothic")
   plt.ylabel("飛距離(m)", fontname="MS Gothic")
   plt.grid()
   plt.show()

   # 結果を表示
   print("最大飛距離{:.2f}mは角度{:.2f}度で達成されます".format(distanceList[np.argmax(distanceList)], degreeList[np.argmax(distanceList)]))

実行結果の一部

参考

• オイラー法/ホイン法/ルンゲクッタ法をつかった常微分方程式の数値解析超入門
  • https://qiita.com/Dr_ASA/items/d44560cf1367290aee6d