Haskellで常微分方程式(オイラー法)
初速 250 km/h で, 45°の角度で打ったボールの軌跡をオイラー法で計算する
(空気抵抗係数を 0.01 で計算)
重力による鉛直方向の減速分は, 重力加速度を g, 時間を t とすると,
空気抵抗による水平方向の減速分は,速度を v, 速度の水平方向成分を vx, 空気抵抗係数を k とすると,
同様に, 鉛直方向の減速分は, 速度の鉛直方向成分を vy とすると,
初期値 x0 から次の式によって, 順次 x1, x2, … を求める.
import Text.Printf -- 重力加速度 g = -9.8 :: Double -- 空気抵抗係数 k = -0.01 :: Double -- 時間間隔(秒) h = 0.01 :: Double -- 空気抵抗による水平方向の減速分 fx::Double->Double->Double fx vx vy = let v = sqrt(vx * vx + vy * vy) in k * v * vx -- 空気抵抗による鉛直方向の減速分 fy::Double->Double->Double fy vx vy = let v = sqrt(vx * vx + vy * vy) in g + (k * v * vy) -- Euler法 euler::Integer->Double->Double->Double->Double->IO () euler i vx vy x y = let t = (fromIntegral i) * h wx = x + h * vx wy = y + h * vy wvx = vx + h * (fx vx vy) wvy = vy + h * (fy vx vy) in if y >= 0.0 then do printf "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f\n" t vx vy wx wy euler (i+1) wvx wvy wx wy else return () main = do -- 角度 let degree = 45.0 :: Double let radian = degree * pi / 180.0 -- 初速 250 km/h -> 秒速に変換 let v = (fromIntegral (250 * 1000 `div` 3600)) -- 水平方向の速度 let vx = v * cos(radian) -- 鉛直方向の速度 let vy = v * sin(radian) -- 位置 let x = 0.0 let y = 0.0 -- Euler法 euler 1 vx vy x y
Z:\>runghc Hs0801.hs 0.01 48.79037 48.79037 0.48790 0.48790 0.02 48.45371 48.35571 0.97244 0.97146 0.03 48.12203 47.92670 1.45366 1.45073 0.04 47.79520 47.50319 1.93161 1.92576 0.05 47.47312 47.08509 2.40634 2.39661 0.06 47.15570 46.67226 2.87790 2.86333 0.07 46.84284 46.26460 3.34633 3.32598 0.08 46.53443 45.86201 3.81167 3.78460 0.09 46.23039 45.46436 4.27398 4.23924 0.10 45.93064 45.07157 4.73328 4.68996 省略 6.20 9.24481 -23.75680 125.68321 2.37768 6.21 9.22125 -23.79424 125.77543 2.13974 6.22 9.19771 -23.83152 125.86740 1.90143 6.23 9.17422 -23.86864 125.95915 1.66274 6.24 9.15076 -23.90561 126.05065 1.42368 6.25 9.12734 -23.94242 126.14193 1.18426 6.26 9.10395 -23.97907 126.23297 0.94447 6.27 9.08060 -24.01557 126.32377 0.70431 6.28 9.05728 -24.05191 126.41435 0.46379 6.29 9.03401 -24.08809 126.50469 0.22291 6.30 9.01076 -24.12412 126.59479 -0.01833
参考文献