roots (empty) → 0.1
raw patch · 12 files changed
+636/−0 lines, 12 filesdep +basedep +taggedsetup-changed
Dependencies added: base, tagged
Files
- Setup.lhs +5/−0
- roots.cabal +40/−0
- src/Math/Root/Bracket.hs +46/−0
- src/Math/Root/Finder.hs +88/−0
- src/Math/Root/Finder/Bisection.hs +38/−0
- src/Math/Root/Finder/Brent.hs +143/−0
- src/Math/Root/Finder/Dekker.hs +53/−0
- src/Math/Root/Finder/FalsePosition.hs +46/−0
- src/Math/Root/Finder/InverseQuadratic.hs +40/−0
- src/Math/Root/Finder/Newton.hs +33/−0
- src/Math/Root/Finder/Ridders.hs +58/−0
- src/Math/Root/Finder/Secant.hs +46/−0
+ Setup.lhs view
@@ -0,0 +1,5 @@+#!/usr/bin/env runhaskell++> import Distribution.Simple+> main = defaultMain+
+ roots.cabal view
@@ -0,0 +1,40 @@+name: roots+version: 0.1+stability: experimental++cabal-version: >= 1.6+build-type: Simple++author: James Cook <mokus@deepbondi.net>+maintainer: James Cook <mokus@deepbondi.net>+license: PublicDomain+homepage: /dev/null++category: Math, Numerical+synopsis: Root-finding algorithms (1-dimensional)+description: Framework for and a few implementations of+ (1-dimensional) numerical root-finding algorithms.++tested-with: GHC == 6.8.3,+ GHC == 6.10.4,+ GHC == 6.12.1, GHC == 6.12.3++source-repository head+ type: darcs+ location: http://code.haskell.org/~mokus/roots++Library+ ghc-options: -Wall+ hs-source-dirs: src+ exposed-modules: Math.Root.Bracket+ Math.Root.Finder+ Math.Root.Finder.Bisection+ Math.Root.Finder.Brent+ Math.Root.Finder.Dekker+ Math.Root.Finder.FalsePosition+ Math.Root.Finder.InverseQuadratic+ Math.Root.Finder.Newton+ Math.Root.Finder.Ridders+ Math.Root.Finder.Secant+ + build-depends: base >= 3 && <5, tagged
+ src/Math/Root/Bracket.hs view
@@ -0,0 +1,46 @@+{-# OPTIONS_GHC -fno-warn-name-shadowing #-}+module Math.Root.Bracket where++-- |Predicate that returns 'True' whenever the given pair of points brackets+-- a root of the given function.+brackets :: (Eq a, Num b) => (a -> b) -> (a,a) -> Bool+brackets f (x1,x2)+ | x1 == x2 = f x1 == 0+ | otherwise = signum (f x1) /= signum (f x2)++-- |@bracket f x1 x2@: Given a function and an initial guessed range x1 to x2,+-- this function expands the range geometrically until a root is bracketed by +-- the returned values, returning a list of the successively expanded ranges. +-- The list will be finite if and only if the sequence yields a bracketing pair.+bracket :: (Fractional a, Num b, Ord b) =>+ (a -> b) -> a -> a -> [(a, a)]+bracket f x1 x2+ | x1 == x2 = error "bracket: empty range"+ | otherwise = go x1 (f x1) x2 (f x2)+ where+ factor = 1.618 -- golden ratio (close enough to it, anyway)+ go x1 f1 x2 f2+ | signum f1 /= signum f2 = [(x1, x2)]+ | abs f1 < abs f2 = (x1, x2) : go x1' (f x1') x2 f2+ | otherwise = (x1, x2) : go x1 f1 x2' (f x2')+ where + x1' = x1 - factor * w+ x2' = x2 + factor * w+ w = x2 - x1++-- |@subdivideAndBracket f x1 x2 n@: Given a function defined on the interval+-- [x1,x2], subdivide the interval into n equally spaced segments and search +-- for zero crossings of the function. The returned list will contain all +-- bracketing pairs found.+subdivideAndBracket :: (Num b, Fractional a, Integral c) =>+ (a -> b) -> a -> a -> c -> [(a, a)]+subdivideAndBracket f x1 x2 n = + [ (x1, x2)+ | ((x1, y1), (x2, y2)) <- zip xys (tail xys)+ , signum y1 /= signum y2+ ]+ where+ dx = (x2 - x1) / fromIntegral n+ xs = x1 : [x1 + dx * fromIntegral i | i <- [1..n]]+ xys = map (\x -> (x, f x)) xs+
+ src/Math/Root/Finder.hs view
@@ -0,0 +1,88 @@+{-# LANGUAGE MultiParamTypeClasses, ScopedTypeVariables, FlexibleContexts #-}+module Math.Root.Finder where++import Control.Monad.Instances ()+import Data.Tagged++-- |General interface for numerical root finders.+class RootFinder r a b where+ -- |@initRootFinder f x0 x1@: Initialize a root finder for the given+ -- function with the initial bracketing interval (x0,x1).+ initRootFinder :: (a -> b) -> a -> a -> r a b+ + -- |Step a root finder for the given function (which should generally + -- be the same one passed to @initRootFinder@), refining the finder's+ -- estimate of the location of a root.+ stepRootFinder :: (a -> b) -> r a b -> r a b+ + -- |Extract the finder's current estimate of the position of a root.+ estimateRoot :: r a b -> a+ + -- |Extract the finder's current estimate of the upper bound of the + -- distance from @estimateRoot@ to an actual root in the function.+ -- + -- Generally, @estimateRoot r@ +- @estimateError r@ should bracket+ -- a root of the function.+ estimateError :: r a b -> a+ + -- |Test whether a root finding algorithm has converged to a given + -- relative accuracy.+ converged :: (Num a, Ord a) => a -> r a b -> Bool+ converged xacc r = abs (estimateError r) <= abs xacc+ + -- |Default number of steps after which root finding will be deemed + -- to have failed. Purely a convenience used to control the behavior+ -- of built-in functions such as 'findRoot' and 'traceRoot'. The+ -- default value is 250.+ defaultNSteps :: Tagged (r a b) Int+ defaultNSteps = Tagged 250++-- |@traceRoot f x0 x1 mbEps@ initializes a root finder and repeatedly+-- steps it, returning each step of the process in a list. When the algorithm+-- terminates or the 'defaultNSteps' limit is exceeded, the list ends.+-- Termination criteria depends on @mbEps@; if it is of the form @Just eps@ +-- then convergence to @eps@ is used (using the @converged@ method of the+-- root finder). Otherwise, the trace is not terminated until subsequent+-- states are equal (according to '=='). This is a stricter condition than+-- convergence to 0; subsequent states may have converged to zero but as long+-- as any internal state changes the trace will continue.+traceRoot :: (Eq (r a b), RootFinder r a b, Num a, Ord a) =>+ (a -> b) -> a -> a -> Maybe a -> [r a b]+traceRoot f a b xacc = go nSteps start (stepRootFinder f start)+ where+ Tagged nSteps = (const :: Tagged a b -> a -> Tagged a b) defaultNSteps start+ start = initRootFinder f a b+ + -- lookahead 1; if tracing with no convergence test, apply a+ -- naive test to bail out if the root stops changing. This is+ -- provided because that's not always the same as convergence to 0,+ -- and the main purpose of this function is to watch what actually+ -- happens inside the root finder.+ go n x next+ | maybe (x==next) (flip converged x) xacc = [x]+ | n <= 0 = []+ | otherwise = x : go (n-1) next (stepRootFinder f next)++-- |@findRoot f x0 x1 eps@ initializes a root finder and repeatedly+-- steps it. When the algorithm converges to @eps@ or the 'defaultNSteps'+-- limit is exceeded, the current best guess is returned, with the @Right@ +-- constructor indicating successful convergence or the @Left@ constructor +-- indicating failure to converge.+findRoot :: (RootFinder r a b, Num a, Ord a) =>+ (a -> b) -> a -> a -> a -> Either (r a b) (r a b)+findRoot f a b xacc = go nSteps start+ where+ Tagged nSteps = (const :: Tagged a b -> a -> Tagged a b) defaultNSteps start+ start = initRootFinder f a b+ + go n x+ | converged xacc x = Right x+ | n <= 0 = Left x+ | otherwise = go (n-1) (stepRootFinder f x)++-- |A useful constant: 'eps' is (for most 'RealFloat' types) the smallest+-- positive number such that @1 + eps /= 1@.+eps :: RealFloat a => a+eps = eps'+ where+ eps' = encodeFloat 1 (1 - floatDigits eps')
+ src/Math/Root/Finder/Bisection.hs view
@@ -0,0 +1,38 @@+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}+module Math.Root.Finder.Bisection+ ( Bisect, bisection+ ) where++import Math.Root.Finder++-- |Bisect an interval in search of a root. At all times, @f (estimateRoot _)@+-- is less than or equal to 0 and @f (estimateRoot _ + estimateError _)@ is +-- greater than or equal to 0.+data Bisect a b = Bisect { _bisX :: !a, _bisF :: !b , _bisDX :: !a }+ deriving (Eq, Ord, Show)++instance (Fractional a, Ord b, Num b) => RootFinder Bisect a b where+ initRootFinder f x1 x2+ | f1 < 0 = Bisect x1 f1 (x2-x1)+ | otherwise = Bisect x2 f2 (x1-x2)+ where f1 = f x1; f2 = f x2+ + stepRootFinder f orig@(Bisect x fx dx) = case fMid `compare` 0 of+ LT -> Bisect xMid fMid dx2+ EQ -> orig+ GT -> Bisect x fx dx2 + where+ dx2 = dx * 0.5+ xMid = x + dx2+ fMid = f xMid+ + estimateRoot (Bisect x _ _) = x+ + estimateError (Bisect _ 0 _) = 0+ estimateError (Bisect _ _ dx) = dx++-- |Using bisection, return a root of a function known to lie between x1 and x2.+-- The root will be refined till its accuracy is +-xacc. If convergence fails,+-- returns the final state of the search.+bisection :: (Ord a, Fractional a, Ord b, Num b) => (a -> b) -> a -> a -> a -> Either (Bisect a b) a+bisection f x1 x2 xacc = fmap estimateRoot (findRoot f x1 x2 xacc)
+ src/Math/Root/Finder/Brent.hs view
@@ -0,0 +1,143 @@+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, FlexibleContexts #-}+{-# OPTIONS_GHC -fno-warn-name-shadowing #-}+module Math.Root.Finder.Brent+ ( Brent+ , brent+ ) where++import Math.Root.Bracket+import Math.Root.Finder+import Data.Maybe+import Text.Printf++-- Invariants:+-- 1) B and C bracket the root+-- 2) |f(B)| <= |f(C)|+-- 3) min(f(B),f(C)) <= f(A) <= max(f(B),f(C))+-- 4) e >= 0+-- |Working state for Brent's root-finding method.+data Brent a b = Brent+ { brA :: !a+ , brFA :: !b+ , brB :: !a+ , brFB :: !b+ , brC :: !a+ , brFC :: !b+ , brE :: a+ } deriving (Eq, Show)++-- TODO: clean up this mess!+instance (RealFloat a, Real b, Fractional b) => RootFinder Brent a b where+ initRootFinder f x1 x2 = fixMagnitudes (Brent x1 f1 x2 f2 x1 f1 dx)+ where f1 = f x1; f2 = f x2; dx = x2 - x1+ + stepRootFinder f r@(Brent a fa b fb c fc e)+ | e >= 2 * min tol1 abs_s -- require that the method be making progress, overall+ && 1.5 * m * signum s >= tol1 + abs_s -- require that the proposed step is getting closer to 'b' - specifically, s should be between 0 and 0.75*(c - b)+ = advance s abs_s+ | otherwise = advance m (abs (b - a))+ where+ -- Minimum step size to continue with inverse-quadratic interpolation+ tol1 = eps * (abs b + 0.5)+ abs_s = abs s+ + -- midpoint for bisection step+ m = 0.5 * (c - b)+ + -- subdivision point for inverse quadratic interpolation step+ s | fa /= fc && fa /= fb+ = let a' = realToFrac (fa / (fc - fb))+ b' = realToFrac (fb / (fc - fa))+ c' = realToFrac (fc / (fb - fa))+ in (a' * b' * c) - ((a' * c' + 1) * b) + (a * b' * c')+ | otherwise+ -- Fall back to linear interpolation when quadratic+ -- interpolation will yield nonsensical results.+ = (c - b) * realToFrac (fb / (fb - fc))+ + -- |Moves the current estimate by 'd' (or by tol1, whichever+ -- is greater) and sets 'brE' to 'e', maintaining all invariants.+ -- Ensuring that at least some tiny jump is made allows quick + -- discovery and termination in the case where the current best+ -- estimate is already nearly on top of the root. Without such+ -- a check, the method would repeatedly tighten the 'c' bound+ -- by bisection every other step, which is really rather stupid+ -- if 'b' is already sitting on a root.+ advance d newE = update b' (f b') newE r+ where+ b' = if abs d > tol1 then b + d else b + tol1 * signum m+++ estimateRoot = brB+ estimateError = brE+ converged _ Brent{brFB = 0} = True+ converged tol Brent{brB = b, brE = e} = + abs e <= 4 * eps * abs b + tol++-- |@brent f x1 x2 xacc@: attempt to find a root of a function known to +-- lie between x1 and x2, using Brent's method. The root will be refined+-- till its accuracy is +-xacc. If convergence fails, returns the final+-- state of the search.+brent :: RealFloat a => (a -> a) -> a -> a -> a -> Either (Brent a a) a+brent f x1 x2 xacc = fmap estimateRoot (findRoot f x1 x2 xacc)++-- |Updates the state by incorporating a new estimate and setting 'brE',+-- maintaining all invariants.+update :: (Num a, Num b, Ord b) => a -> b -> a -> Brent a b -> Brent a b+update b fb e r@Brent{brB = a, brFB = fa} + = fixMagnitudes (fixSigns r{brA = a, brFA = fa, brB = b, brFB = fb, brE = e})++-- Establish invariant (1) that b and c bracket the root,+-- based on precondition that (a,c) already does.+-- +-- (a,c) brackets implies that either (b,c) or (a,b) brackets. In the +-- former case, nothing needs to be done as (by construction) either fb is already+-- between fa and fc or b is already between a and c (depending which kind of +-- step was taken). In the latter case, discard C and use A in its place, because+-- c and fc are both (by the existing invariants - (a,c) bracket, |f(c)| >= |f(a)|) +-- outside the new region of interest.+fixSigns :: (Num a, Num b, Ord b) => Brent a b -> Brent a b+fixSigns br@Brent{ brA = a+ , brFA = fa, brFB = fb, brFC = fc }+ | (fb > 0 && fc > 0) || (fb < 0 && fc < 0)+ = br { brC = a, brFC = fa }+ | otherwise + = br++-- Establish invariant (2) that |f(c)| >= |f(b)| and invariant (3) that+-- 'fa' falls between fb and fc.+fixMagnitudes :: (Num b, Ord b) => Brent a b -> Brent a b+fixMagnitudes br@Brent{ brC = c, brB = b+ , brFC = fc, brFB = fb }+ | abs fc < abs fb+ = br { brA = b, brFA = fb+ , brB = c, brFB = fc+ , brC = b, brFC = fb+ }+ | otherwise + = br++-- |debugging function to show a nice trace of the progress of the algorithm+_traceBrent :: (PrintfArg a, RealFloat a,+ PrintfArg b, Ord b, Num b,+ RootFinder Brent a b) =>+ (a -> b) -> Maybe (a, a) -> IO ()+_traceBrent f mbRange = do+ xs <- sequence+ [ put br >> return br+ | br <- traceRoot f x0 x1 (Just eps)+ ]++ putStrLn "(converged)"+ go (last xs)+ where + (x0,x1) = fromMaybe (last (bracket f 0 1)) mbRange+ put Brent{brA=a, brB=b, brC=c, brFA=fa, brFB=fb, brFC=fc} = + putStrLn . map fixPlus $+ printf (concat (replicate 6 "%-+25g")) a b c fa fb fc+ fixPlus '+' = ' '+ fixPlus c = c+ go x + | x == x' = return ()+ | otherwise = put x >> go x'+ where x' = stepRootFinder f x
+ src/Math/Root/Finder/Dekker.hs view
@@ -0,0 +1,53 @@+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}+module Math.Root.Finder.Dekker (Dekker, dekker) where++import Math.Root.Finder++-- fields: a fa b fb oldB oldFb+-- invariants: +-- 1) signum fA /= signum fB+-- 2) abs fB <= abs fA+-- 3) (oldB-a)*(oldB-b) >= 0+data Dekker a b = Dekker !a !b !a !b !a !b deriving (Eq, Show)++-- |@dekker f x1 x2 xacc@: attempt to find a root of a function known to +-- lie between x1 and x2, using Dekker's method. The root will be refined+-- till its accuracy is +-xacc. If convergence fails, returns the final+-- state of the search.+dekker :: RealFloat a => (a -> a) -> a -> a -> a -> Either (Dekker a a) a+dekker f x1 x2 xacc = fmap estimateRoot (findRoot f x1 x2 xacc)++instance (Fractional a, Ord a, Real b, Fractional b, Ord b) => RootFinder Dekker a b where+ initRootFinder f x0 x1 + | signum f0 == signum f1 = error "initRootFinder/Dekker: starting points do not (obviously) bracket a root"+ | abs f0 <= abs f1 = Dekker x1 f1 x0 f0 x1 f1+ | otherwise = Dekker x0 f0 x1 f1 x0 f0+ where f0 = f x0; f1 = f x1+ + stepRootFinder f orig@(Dekker a _ b fb oldB oldFb)+ | fb == 0 = orig+ | oldFb /= fb + && s `between` (a,b) = step s (f s) orig+ | otherwise = step m (f m) orig+ where+ s = b - (b * oldB) * realToFrac (fb / (fb - oldFb))+ m = 0.5 * (a + b)+ + estimateRoot (Dekker _ _ b _ _ _) = b+ estimateError (Dekker a _ b _ _ _) = a - b++between :: Ord a => a -> (a,a) -> Bool+a `between` (x,y) = (a > min x y) && (a < max x y)++-- |Incorporates a new point, maintaining invariant 1, assuming invariant 3,+-- and using 'accept' to restore invariant 2.+step :: (Num b, Ord b) => a -> b -> Dekker a b -> Dekker a b+step x fx orig@(Dekker a fa b fb _ _)+ | signum fx /= signum fa = accept a fa x fx orig+ | otherwise = accept x fx b fb orig++-- |Re-establishes invariant 2 (abs fb <= abs fa) without affecting invariants 1 and 3.+accept :: (Num b, Ord b) => a -> b -> a -> b -> Dekker a b -> Dekker a b+accept a fa b fb (Dekker _ _ oldB oldFb _ _)+ | abs fb <= abs fa = Dekker a fa b fb oldB oldFb+ | otherwise = Dekker b fb a fa oldB oldFb
+ src/Math/Root/Finder/FalsePosition.hs view
@@ -0,0 +1,46 @@+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}+module Math.Root.Finder.FalsePosition+ ( FalsePosition, falsePosition+ ) where++import Math.Root.Finder++-- | @falsePosition f x1 x2 xacc@: Using the false-position method, return a+-- root of a function known to lie between x1 and x2. The root is refined +-- until its accuracy is += xacc.+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)++-- |Iteratively refine a bracketing interval [x1, x2] of a root of f+-- until total convergence (which may or may not ever be achieved) using +-- the false-position method.+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+ -- step once to compute first estimate+ | 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
+ src/Math/Root/Finder/InverseQuadratic.hs view
@@ -0,0 +1,40 @@+{-# LANGUAGE+ MultiParamTypeClasses,+ FlexibleInstances+ #-}+module Math.Root.Finder.InverseQuadratic (InverseQuadratic, inverseQuadratic) where++import Math.Root.Finder++data InverseQuadratic a b = InverseQuadratic !a !b !a !b !a !b+ deriving (Eq, Show)++-- |@inverseQuadratic f x1 x2 xacc@: attempt to find a root of a function +-- known to lie between x1 and x2, using the inverse quadratic interpolation +-- method. The root will be refined till its accuracy is +-xacc. If+-- convergence fails, returns the final state of the search.+inverseQuadratic :: RealFloat a => (a -> a) -> a -> a -> a -> Either (InverseQuadratic a a) a+inverseQuadratic f x1 x2 xacc = fmap estimateRoot (findRoot f x1 x2 xacc)++instance (Fractional a, Ord a, Real b, Fractional b) => RootFinder InverseQuadratic a b where+ initRootFinder f x1 x2 = InverseQuadratic x0 (f x0) x1 (f x1) x2 (f x2)+ where x0 = 0.5 * (x1 + x2)+ stepRootFinder f orig@(InverseQuadratic x0 f0 x1 f1 x2 f2)+ | f1 /= f2 = InverseQuadratic newX newF x0 f0 x1 f1+ | otherwise = orig+ where+ newX + | f0 /= f1 && f0 /= f2 + = let a = realToFrac (f0 / (f2 - f1))+ b = realToFrac (f1 / (f2 - f0))+ c = realToFrac (f2 / (f1 - f0))+ in (a * b * x2) - (a * c * x1) + (b * c * x0)+ | otherwise+ -- Fall back to secant method (linear interpolation)+ -- when quadratic interpolation will yield nonsensical results.+ = x1 - realToFrac f1 * (x1 - x2) / realToFrac (f1 - f2)+ newF = f newX+ + estimateRoot (InverseQuadratic x0 _ _ _ _ _) = x0+ estimateError (InverseQuadratic x0 _ x1 _ x2 _) = + maximum [x0, x1, x2] - minimum [x0, x1, x2]
+ src/Math/Root/Finder/Newton.hs view
@@ -0,0 +1,33 @@+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}+module Math.Root.Finder.Newton+ ( Newton, newton+ ) where++import Math.Root.Finder++data Newton a b = Newton+ { newtRTN :: !a+ , newtDX :: a+ } deriving (Eq, Show)++instance Fractional a => RootFinder Newton a (a,a) where+ initRootFinder f x1 x2 = stepRootFinder f (Newton rtn undefined)+ where+ rtn = 0.5 * (x1 + x2)+ + stepRootFinder f Newton{newtRTN = rtn} = Newton (rtn - dx) dx+ where+ (y,dy) = f rtn+ dx = y / dy+ + estimateRoot Newton{newtRTN = rtn} = rtn+ estimateError Newton{newtDX = dx} = dx ++-- | @newton f x1 x2 xacc@: using Newton's method, return a root of a+-- function known to lie between x1 and x2. The root is refined until its+-- accuracy is += xacc.+-- +-- The function passed should return a pair containing the value of the+-- function and its derivative, respectively.+newton :: (Ord a, Fractional a) => (a -> (a, a)) -> a -> a -> a -> Either (Newton a (a,a)) a+newton f x1 x2 xacc = fmap estimateRoot (findRoot f x1 x2 xacc)
+ src/Math/Root/Finder/Ridders.hs view
@@ -0,0 +1,58 @@+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}+{-# OPTIONS_GHC -fno-warn-name-shadowing #-}+module Math.Root.Finder.Ridders+ ( RiddersMethod, ridders+ ) where++import Math.Root.Finder++-- |@ridders f x1 x2 xacc@: attempt to find a root of a function known to +-- lie between x1 and x2, using Ridders' method. The root will be refined+-- till its accuracy is +-xacc. If convergence fails, returns the final+-- state of the search.+ridders :: (Ord a, Floating a) => (a -> a) -> a -> a -> a -> Either (RiddersMethod a a) a+ridders f x1 x2 xacc = fmap estimateRoot (findRoot f x1 x2 xacc)++data RiddersMethod a b+ = ConvergedRidders !a+ | RiddersMethod+ { ridXL :: !a+ , _ridFL :: !b+ , ridXH :: !a+ , _ridFH :: !b+ } deriving (Eq, Show)++instance (Floating a, Ord a) => RootFinder RiddersMethod a a where+ initRootFinder f x1 x2+ | f1 < 0 && f2 < 0+ || f2 > 0 && f1 > 0 = error "riddersMethod: interval does not bracket a root"+ | otherwise = RiddersMethod x1 f1 x2 f2+ where+ f1 = f x1+ f2 = f x2+ stepRootFinder _ orig@ConvergedRidders{} = orig+ stepRootFinder f (RiddersMethod xl fl xh fh)+ | signNEQ fm fNew = finish xNew fNew xm fm+ | signNEQ fl fNew = finish xNew fNew xl fl+ | signNEQ fh fNew = finish xNew fNew xh fh+ | otherwise = error "RiddersMethod: encountered singularity"+ where+ xm = 0.5 * (xl + xh)+ fm = f xm+ s = sqrt (fm*fm - fl*fh)+ xNew = xm + (xm-xl)*((if fl >= fh then id else negate) fm / s)+ fNew = f xNew+ + signNEQ a b = a /= 0 && signum b /= signum a+ + finish xl fl xh fh+ | xl == xh = ConvergedRidders xl+ | fl == 0 = ConvergedRidders xl+ | fh == 0 = ConvergedRidders xh+ | otherwise = RiddersMethod xl fl xh fh+ + estimateRoot (ConvergedRidders x) = x+ estimateRoot RiddersMethod{ridXL = x} = x+ + estimateError ConvergedRidders{} = 0+ estimateError RiddersMethod{ridXL = xl, ridXH = xh} = xl - xh
+ src/Math/Root/Finder/Secant.hs view
@@ -0,0 +1,46 @@+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}+module Math.Root.Finder.Secant+ ( SecantMethod, secant+ ) where++import Math.Root.Finder++-- | @secant f x1 x2 xacc@: Using the secant method, return the root of a+-- function thought to lie between x1 and x2. The root is refined until its+-- accuracy is +-xacc.+secant :: (Ord a, Fractional a) => (a -> a) -> a -> a -> a -> Either (SecantMethod a a) a+secant f x1 x2 xacc = fmap estimateRoot (findRoot f x1 x2 xacc)++-- |Iteratively refine 2 estimates x1, x2 of a root of f until total +-- convergence (which may or may not ever be achieved) using the+-- secant method.+data SecantMethod a b+ = ConvergedSecantMethod !a+ | SecantMethod+ { secDX :: !a+ , secXL :: !a+ , _secFL :: !b+ , _secRTS :: !a+ , _secFRTS :: !b+ } deriving (Eq, Show)++instance (Fractional a, Ord a) => RootFinder SecantMethod a a where+ initRootFinder f x1 x2+ | abs f1 < abs f2 = stepRootFinder f $ SecantMethod 0 x2 f2 x1 f1+ | otherwise = stepRootFinder f $ SecantMethod 0 x1 f1 x2 f2+ where f1 = f x1; f2 = f x2+ + stepRootFinder _ orig@ConvergedSecantMethod{} = orig+ stepRootFinder f (SecantMethod _ xl fl rts fRts)+ | fNew == 0 = ConvergedSecantMethod xNew+ | otherwise = SecantMethod dx rts fRts xNew fNew+ where+ dx = (xl - rts) * fRts / (fRts - fl)+ xNew = rts + dx+ fNew = f xNew+ + estimateRoot (ConvergedSecantMethod x) = x+ estimateRoot SecantMethod{secXL = x} = x+ + estimateError ConvergedSecantMethod{} = 0+ estimateError SecantMethod{secDX = dx} = dx