numeric-tools-0.2.0.0: Numeric/Tools/Equation.hs
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE DeriveDataTypeable #-}
-- |
-- Module : Numeric.Tools.Equation
-- Copyright : (c) 2011 Aleksey Khudyakov
-- License : BSD3
--
-- Maintainer : Aleksey Khudyakov <alexey.skladnoy@gmail.com>
-- Stability : experimental
-- Portability : portable
--
-- Numerical solution of ordinary equations.
module Numeric.Tools.Equation (
-- * Data type
Root(..)
, fromRoot
-- * Equations solversv
, solveBisection
, solveRidders
, solveNewton
-- * References
-- $references
) where
import Numeric.ApproxEq
import Control.Applicative
import Control.Monad (MonadPlus(..), ap)
import Data.Typeable (Typeable)
-- | The result of searching for a root of a mathematical function.
data Root a = NotBracketed
-- ^ The function does not have opposite signs when
-- evaluated at the lower and upper bounds of the search.
| SearchFailed
-- ^ The search failed to converge to within the given
-- error tolerance after the given number of iterations.
| Root a
-- ^ A root was successfully found.
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
-- | Returns either the result of a search for a root, or the default
-- value if the search failed.
fromRoot :: a -- ^ Default value.
-> Root a -- ^ Result of search for a root.
-> a
fromRoot _ (Root a) = a
fromRoot a _ = a
-- | Use bisection method to compute root of function.
--
-- The function must have opposite signs when evaluated at the lower
-- and upper bounds of the search (i.e. the root must be bracketed).
solveBisection :: Double -- ^ Required absolute precision
-> (Double,Double) -- ^ Range
-> (Double -> Double) -- ^ Equation
-> 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 function. Preconditions:
-- f a < 0
-- f b > 0
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
-- | Use the method of Ridders to compute a root of a function.
--
-- The function must have opposite signs when evaluated at the lower
-- and upper bounds of the search (i.e. the root must be bracketed).
solveRidders :: Double -- ^ Absolute error tolerance.
-> (Double,Double) -- ^ Lower and upper bounds for the search.
-> (Double -> Double) -- ^ Function to find the roots of.
-> Root Double
solveRidders eps (lo,hi) f
| flo == 0 = Root lo
| fhi == 0 = Root hi
| flo*fhi > 0 = NotBracketed -- root is not bracketed
| otherwise = go lo flo hi fhi 0
where
go !a !fa !b !fb !i
-- Root is bracketed within 1 ulp. No improvement could be made
| within 1 a b = Root a
-- Root is found. Check that f(m) == 0 is nessesary to ensure
-- that root is never passed to 'go'
| fm == 0 = Root m
| fn == 0 = Root n
| d < eps = Root n
-- Too many iterations performed. Fail
| i >= (100 :: Int) = SearchFailed
-- Ridder's approximation coincide with one of old
-- bounds. Revert to bisection
| n == a || n == b = case () of
_| fm*fa < 0 -> go a fa m fm (i+1)
| otherwise -> go m fm b fb (i+1)
-- Proceed as usual
| 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
-- | Solve equation using Newton-Raphson method. Root must be
-- bracketed. If Newton's step jumps outside of bracket or do not
-- converge sufficiently fast function reverts to bisection.
solveNewton :: Double -- ^ Absolute error tolerance
-> (Double,Double) -- ^ Lower and upper bounds for root
-> (Double -> Double) -- ^ Function
-> (Double -> Double) -- ^ Function's derivative
-> 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*(hi-lo)
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 function
worker i dxOld x fx fx' a b
-- Convergence achieved
| fx == 0 = Root x -- x is a root
| within 1 a b = Root x -- Bracket is too tight. No improvements could be made
| abs (a - b) < eps = Root x --
| abs dx < eps = Root x' -- Precision achieved
| within 0 x x' = Root x' -- Newton step doesn't improve solution
-- Too many iterations
| i > (100::Int) = SearchFailed
-- Newton step jumped out of range or step decreases too slow.
-- Revert to bisection.
-- NOTE this guard is selected if dx evaluated to NaN or ±∞.
| 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
-- Normal newton step
| otherwise =
if fx < 0 then step dx x'
else step dx x'
where
dx = fx / fx'
x' = x - dx
-- Perform one step
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
{-# INLINABLE solveNewton #-}
-- $references
--
-- * Ridders, C.F.J. (1979) A new algorithm for computing a single
-- root of a real continuous function.
-- /IEEE Transactions on Circuits and Systems/ 26:979–980.