Ridders method

Ridders.hs

module Ridders where

import qualified Data.Map as M
import Control.Monad.State

type R = Double

ridders :: R -> R -> R -> (R -> R) -> Maybe R
ridders ε a b fn = evalState (go a b) M.empty
  where
    calc x = do
      mem <- get
      case M.lookup x mem of
        Nothing -> do
                   let y = fn x
                   put $ M.insert x y mem
                   return y
        Just y ->  return y

    go :: R -> R -> State (M.Map R R) (Maybe R)
    go one two = do
      let x0 = min one two
          x2 = max one two
          x1 = (x0 + x2)/2
          d  = x1 - x0
      f0 <- calc x0
      f1 <- calc x1
      f2 <- calc x2
      let sgn = if f0 > f2 then 1 else -1
          x3 = x1 + d*sgn*f1/sqrt(f1*f1 - f0*f2)
      f3 <- calc x3
      case () of
        _ | abs f0 < ε -> return $ Just x0
        _ | abs f1 < ε -> return $ Just x1
        _ | abs f2 < ε -> return $ Just x2
        _ | abs f3 < ε -> return $ Just x3
        _ | f3 * f0 < 0 -> go x3 x0
        _ | f3 * f1 < 0 -> go x3 x1
        _ | f3 * f2 < 0 -> go x3 x2
        _ -> return Nothing