積分をあれこれ

calc.hs

import Control.Monad

linearScale :: (Fractional a, Integral b) => (a, a) -> b -> ([a], a)
linearScale (begin, end) n = (xs, dx)
  where
    dx = (end - begin) / (fromIntegral n)
    xs = [(fromIntegral n) * dx + begin | n <- [0 .. ]]

forwardEuler :: (Fractional a, Integral b) => (a, a) -> b -> (a -> a) -> [a]
forwardEuler range n f = ys
  where
    (xs, dx) = linearScale range n
    ys = scanl (\l r -> l + (f r) * dx) 0.0 xs

backwardEuler :: (Fractional a, Integral b) => (a, a) -> b -> (a -> a) -> [a]
backwardEuler range n f = ys
  where
    (xs, dx) = linearScale range n
    ys = scanl (\l r -> l + (f r) * dx) 0.0 $ tail xs

trapezoid :: (Fractional a, Integral b) => (a, a) -> b -> (a -> a) -> [a]
trapezoid range n f = ys
  where
    (xs, dx) = linearScale range n
    g = \l (r0, r1) -> l + ((f r0) + (f r1)) * 0.5 * dx
    xx = zip xs (tail xs)
    ys = scanl g 0 xx

example f range n = [exs, fes, bes]
  where
    exs   = [x * x | x <- fst $ linearScale range n]
    fes   = forwardEuler  range n f
    bes   = backwardEuler range n f
    
printResults :: (Show a) => Int -> [[a]] -> IO ()
printResults n (xs:ys:zs:[]) = forM_ [0 .. n] $ \i -> do
    let x = xs !! i
        y = ys !! i
        z = zs !! i
    putStrLn $ (show x) ++ " " ++ (show y) ++ " " ++ (show z)

main = do
    let f     = \x -> 2.0 * x
        range = (0, 1)
        n     = 10
        results = example f range n
        mids = trapezoid range n f
    forM_ [0 .. n] $ \i -> do
        putStrLn $ show i ++ " " ++ show (mids !! i)