module Numeric.Tools.Equation (
Root(..)
, fromRoot
, solveBisection
, solveRidders
, solveNewton
) where
import Numeric.ApproxEq
import Control.Applicative
import Control.Monad (MonadPlus(..), ap)
import Data.Typeable (Typeable)
data Root a = NotBracketed
| SearchFailed
| Root a
deriving (Eq, Read, Show, Typeable)
instance Functor Root where
fmap _ NotBracketed = NotBracketed
fmap _ SearchFailed = SearchFailed
fmap f (Root a) = Root (f a)
instance Monad Root where
NotBracketed >>= _ = NotBracketed
SearchFailed >>= _ = SearchFailed
Root a >>= m = m a
return = Root
instance MonadPlus Root where
mzero = SearchFailed
r@(Root _) `mplus` _ = r
_ `mplus` p = p
instance Applicative Root where
pure = Root
(<*>) = ap
instance Alternative Root where
empty = SearchFailed
r@(Root _) <|> _ = r
_ <|> p = p
fromRoot :: a
-> Root a
-> a
fromRoot _ (Root a) = a
fromRoot a _ = a
solveBisection :: Double
-> (Double,Double)
-> (Double -> Double)
-> Root Double
solveBisection eps (lo,hi) f
| flo * fhi > 0 = NotBracketed
| flo == 0 = Root lo
| fhi == 0 = Root hi
| flo < 0 = worker 0 lo hi
| otherwise = worker 0 hi lo
where
flo = f lo
fhi = f hi
worker i a b
| within 1 a b = Root a
| abs (b a) <= eps = Root c
| fc == 0 = Root c
| i >= (100 :: Int) = SearchFailed
| fc < 0 = worker (i+1) c b
| otherwise = worker (i+1) a c
where
c = 0.5 * (a + b)
fc = f c
solveRidders :: Double
-> (Double,Double)
-> (Double -> Double)
-> Root Double
solveRidders eps (lo,hi) f
| flo == 0 = Root lo
| fhi == 0 = Root hi
| flo*fhi > 0 = NotBracketed
| otherwise = go lo flo hi fhi 0
where
go !a !fa !b !fb !i
| within 1 a b = Root a
| fm == 0 = Root m
| fn == 0 = Root n
| d < eps = Root n
| i >= (100 :: Int) = SearchFailed
| n == a || n == b = case () of
_| fm*fa < 0 -> go a fa m fm (i+1)
| otherwise -> go m fm b fb (i+1)
| fn*fm < 0 = go n fn m fm (i+1)
| fn*fa < 0 = go a fa n fn (i+1)
| otherwise = go n fn b fb (i+1)
where
d = abs (b a)
dm = (b a) * 0.5
!m = a + dm
!fm = f m
!dn = signum (fb fa) * dm * fm / sqrt(fm*fm fa*fb)
!n = m signum dn * min (abs dn) (abs dm 0.5 * eps)
!fn = f n
!flo = f lo
!fhi = f hi
solveNewton :: Double
-> (Double,Double)
-> (Double -> Double)
-> (Double -> Double)
-> Root Double
solveNewton eps (lo,hi) f f'
| flo == 0 = Root lo
| fhi == 0 = Root hi
| flo * fhi > 0 = NotBracketed
| otherwise = let d = 0.5*(hilo)
mid = 0.5*(lo+hi)
fun = worker 0 d mid (f mid) (f' mid)
in if flo < 0 then fun lo hi
else fun hi lo
where
flo = f lo
fhi = f hi
worker i dxOld x fx fx' a b
| fx == 0 = Root x
| within 1 a b = Root x
| abs (a b) < eps = Root x
| abs dx < eps = Root x'
| within 0 x x' = Root x'
| i > (100::Int) = SearchFailed
| not ((a x')*(b x') < 0) || abs (dxOld / dx) < 2 =
let dx' = 0.5 * (b a)
mid = 0.5 * (a + b)
in if fx < 0 then step dx' mid
else step dx' mid
| otherwise =
if fx < 0 then step dx x'
else step dx x'
where
dx = fx / fx'
x' = x dx
step dy y
| fy < 0 = worker (i+1) dy y fy fy' y b
| otherwise = worker (i+1) dy y fy fy' a y
where fy = f y
fy' = f' y