module Math.Root.Finder.FalsePosition
( FalsePosition, falsePosition
) where
import Math.Root.Finder
falsePosition :: (Ord a, Fractional a) => (a -> a) -> a -> a -> a -> Either (FalsePosition a a) a
falsePosition f x1 x2 xacc = fmap estimateRoot (findRoot f x1 x2 xacc)
data FalsePosition a b = FalsePosition
{ fpRoot :: !a
, fpDX :: !a
, _fpXL :: !a
, _fpFL :: !a
, _fpXH :: !a
, _fpFH :: !a
} deriving (Eq, Show)
instance (Fractional a, Ord a) => RootFinder FalsePosition a a where
initRootFinder f x1 x2
| f1 <= 0 && f2 >= 0
|| f2 <= 0 && f1 >= 0 = stepRootFinder f $ FalsePosition 0 0 x2 f2 x1 f1
| otherwise = error "FalsePosition: given interval does not bracket a root"
where
f1 = f x1
f2 = f x2
stepRootFinder f (FalsePosition _ _ xl fl xh fh) = case compare fNew 0 of
LT -> FalsePosition xNew (xl xNew) xNew fNew xh fh
EQ -> FalsePosition xNew 0 xNew fNew xNew fNew
GT -> FalsePosition xNew (xh xNew) xl fl xNew fNew
where
dx = xh xl
xNew = xl + dx * fl / (fl fh)
fNew = f xNew
estimateRoot = fpRoot
estimateError = fpDX