練習プログラムその2(Haskellでルンゲクッタ)

plotrk.png

rk.hs

xt :: Double -> Double -> Double -> Double
xt t x y = y

yt :: Double -> Double -> Double -> Double
yt t x y = -2*x + sin (2*t)

rk :: Int -> Double -> Double -> Double -> IO()
rk i t x y =
  let
    h = 0.01
    filename = "rk.dat"
    count = 10000

    k1x = xt t x y
    k1y = yt t x y
    k2x = xt (t+h/2.0) (x+k1x*h/2.0) (y+k1y*h/2.0)
    k2y = yt (t+h/2.0) (x+k1x*h/2.0) (y+k1y*h/2.0)
    k3x = xt (t+h/2.0) (x+k2x*h/2.0) (y+k2y*h/2.0)
    k3y = yt (t+h/2.0) (x+k2x*h/2.0) (y+k2y*h/2.0)
    k4x = xt (t+h) (x+k3x*h) (y+k3y*h)
    k4y = yt (t+h) (x+k3x*h) (y+k3y*h)
    totx = x + (k1x + 2*k2x + 2*k3x + k4x)*h/6.0
    toty = y + (k1y + 2*k2y + 2*k3y + k4y)*h/6.0
    tott = t + h
  in do
    let str = show tott ++"  "++ show totx ++"  "++ show toty
    appendFile filename (str++"\n");
    if i == count then return()
    else rk (i+1) tott totx toty

main :: IO()
main =
  let
    inix = 0.00
    iniy = 1.00
    init' = 0.00
    filename = "rk.dat"
  in do
    writeFile filename (show init' ++"  "++ show inix ++"  "++ show iniy ++"\n");
    rk 0 init' inix iniy