ONLY DO WHAT ONLY YOU CAN DO

こけたら立ちなはれ 立ったら歩きなはれ

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
参考文献