各種オイラー法

test0.hs

import Graphics.Gnuplot.Simple
import Data.IORef
import Control.Monad

forwardEulerMethod :: Fractional a => (a -> a) -> Integer -> (a, a) -> [(a, a)]
forwardEulerMethod f divNum range = zip (tail xs) accums
  where
    xs = linearScale divNum range
    ys = map f xs
    dx = xs !! 1 - xs !! 0
    ss = map (\y -> y * dx) ys
    accums = scanl1 (+) ss

backwardEulerMethod :: Fractional a => (a -> a) -> Integer -> (a, a) -> [(a, a)]
backwardEulerMethod f divNum range = zip xs accums
  where
    xs = linearScale divNum range
    ys = map f xs
    dx = xs !! 1 - xs !! 0
    ss = map (\y -> y * dx) ys
    accums = scanl1 (+) ss

middlePointMethod :: Fractional a => (a -> a) -> Integer -> (a, a) -> [(a, a)]
middlePointMethod f divNum range = zip (tail xs) accums
  where
    xs = linearScale divNum range
    ys = map f xs
    ms = map (\(y0, y1) -> (y1 + y0) / 2.0) $ zip (tail ys) ys
    dx = xs !! 1 - xs !! 0
    ss = map (\y -> y * dx) ms
    accums = scanl1 (+) ss

trapezoidalApproximation :: Fractional a => (a -> a) -> Integer -> (a, a) -> [(a, a)]
trapezoidalApproximation f divNum range = zip (tail xs) accums
  where
    xs = linearScale divNum range
    dx = xs !! 1 - xs !! 0
    ys = map f xs
    ss = map (\(u,l) -> (u + l) * dx / 2.0) $ zip (tail ys) ys
    accums = scanl1 (+) ss

runTest :: IO ()
runTest = do
    let forwardResults   = forwardEulerMethod        f 10 (0,1)
        backwardResults  = backwardEulerMethod       f 10 (0,1)
        middPointResults = middlePointMethod         f 10 (0,1)
        trapResults      = trapezoidalApproximation  f 10 (0,1)
    putStrLn "forward euler method:"
    mapM_ printResult forwardResults
    putStrLn "backward euler method:"
    mapM_ printResult backwardResults
    putStrLn "middle point method:"
    mapM_ printResult middPointResults
    putStrLn "trapezoidal Approximation:"
    mapM_ printResult trapResults
  where
    printResult (x, y) = putStrLn $ (show x) ++ " " ++ (show y)

f :: Fractional a => a -> a
f t = 2.0 * t

doForwardEulerMethod f n = do
    let dt = 1.0 / n
        x0 = 0.0
    x <- newIORef x0
    putStrLn $ "0.0 " ++ (show x0)
    forM_ [0 .. n - 1] $ \i -> do
        let t1 = i * dt
            t2 = (i + 1) * dt
        modifyIORef x (\x -> (f t1) * dt + x)
        xNew <- readIORef x
        putStrLn $ (show t2) ++ " " ++ (show xNew)

main = runTest