diff --git a/Numeric/ApproxEq.hs b/Numeric/ApproxEq.hs new file mode 100644 --- /dev/null +++ b/Numeric/ApproxEq.hs @@ -0,0 +1,81 @@ +-- | +-- Module : Numeric.ApproxEq +-- Copyright : (c) 2011 Aleksey Khudyakov, Bryan O'Sullivan +-- License : BSD3 +-- +-- Maintainer : Aleksey Khudyakov +-- Stability : experimental +-- Portability : portable +-- +-- Different implementations of approximate equality for floating +-- point values. There are multiple ways to implement approximate +-- equality. They have different semantics and it's up to programmer +-- to choose right one. +module Numeric.ApproxEq ( + eqRelative + , eqRelCompl + , eqAbsolute + , within + ) where + +import Control.Monad.ST (runST) + +import Data.Complex (Complex(..)) +import Data.Primitive.ByteArray (newByteArray, readByteArray, writeByteArray) +import Data.Int (Int64) + + + +-- | Relative difference between two numbers are less than predefined +-- value. For example 1 is approximately equal to 1.0001 with 1e-4 +-- precision. Same is true for 10000 and 10001. +-- +-- This method of camparison doesn't work for numbers which are +-- approximately 0. 'eqAbsolute' should be used instead. +eqRelative :: (Fractional a, Ord a) + => a -- ^ Relative precision + -> a -> a -> Bool +eqRelative eps a b = abs (a - b) <= abs (eps * a) +{-# INLINE eqRelative #-} + +-- | Relative equality for comlex numbers. +eqRelCompl :: (RealFloat a, Ord a) + => a -- ^ Relative precision + -> Complex a -> Complex a -> Bool +eqRelCompl eps a b = d <= eps * r + where + (d :+ _) = abs (a - b) + (r :+ _) = abs a +{-# INLINE eqRelCompl #-} + +-- | Difference between values is less than specified precision. +eqAbsolute :: (Num a, Ord a) + => a -- ^ Absolute precision + -> a -> a -> Bool +eqAbsolute eps a b = abs (a - b) <= eps +{-# INLINE eqAbsolute #-} + + +-- | Compare two 'Double' values for approximate equality, using +-- Dawson's method. +-- +-- The required accuracy is specified in ULPs (units of least +-- precision). If the two numbers differ by the given number of ULPs +-- or less, this function returns @True@. +-- +-- Algorithm is based on Bruce Dawson's \"Comparing floating point +-- numbers\": +-- +within :: Int -- ^ Number of ULPs of accuracy desired. + -> Double -> Double -> Bool +within ulps a b = runST $ do + buf <- newByteArray 8 + ai0 <- writeByteArray buf 0 a >> readByteArray buf 0 + bi0 <- writeByteArray buf 0 b >> readByteArray buf 0 + let big = 0x8000000000000000 :: Int64 + ai | ai0 < 0 = big - ai0 + | otherwise = ai0 + bi | bi0 < 0 = big - bi0 + | otherwise = bi0 + return $! abs (ai - bi) <= fromIntegral ulps + diff --git a/Numeric/Tools/Equation.hs b/Numeric/Tools/Equation.hs --- a/Numeric/Tools/Equation.hs +++ b/Numeric/Tools/Equation.hs @@ -1,3 +1,5 @@ +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE DeriveDataTypeable #-} -- | -- Module : Numeric.Tools.Equation -- Copyright : (c) 2011 Aleksey Khudyakov @@ -9,36 +11,204 @@ -- -- Numerical solution of ordinary equations. module Numeric.Tools.Equation ( - solveBisection + -- * Data type + Root(..) + , fromRoot + -- * Equations solversv + , solveBisection + , solveRidders + , solveNewton + -- * References + -- $references ) where -import Numeric.IEEE (epsilon) +import Numeric.ApproxEq +import Control.Applicative +import Control.Monad (MonadPlus(..), ap) +import Data.Typeable (Typeable) --- | Solve equation @f(x) = 0@ using bisection method. Function is --- must be continous. If function has different signs at the ends of --- initial interval answer is always returned. 'Nothing' is returned --- if function fails to find an answer. + +-- | 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 - -> Maybe Double -solveBisection eps (a,b) f - | a >= b = Nothing - | fa * fb > 0 = Nothing - | otherwise = Just $ bisectionWorker (abs eps) f a b fa fb + -> 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 - fa = f a - fb = f b + 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 -bisectionWorker :: Double -> (Double -> Double) -> Double -> Double -> Double -> Double -> Double -bisectionWorker eps f a b fa fb - | (b - a) <= eps = c - | (b - a) / b <= epsilon = c - | fa * fc < 0 = bisectionWorker eps f a c fa fc - | otherwise = bisectionWorker eps f c b fc fb + +-- | 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 - c = 0.5 * (a + b) - fc = f c + 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. diff --git a/Numeric/Tools/Integration.hs b/Numeric/Tools/Integration.hs --- a/Numeric/Tools/Integration.hs +++ b/Numeric/Tools/Integration.hs @@ -42,6 +42,7 @@ import qualified Data.Vector.Unboxed as U import qualified Data.Vector.Unboxed.Mutable as M +import qualified Numeric.IEEE as IEEE ---------------------------------------------------------------- @@ -75,6 +76,7 @@ -- | Result of numeric integration. data QuadRes = QuadRes { quadRes :: Maybe Double -- ^ Integraion result + , quadBestEst :: Double -- ^ Best estimate of integral , quadPrecEst :: Double -- ^ Rough estimate of attained precision , quadNIter :: Int -- ^ Number of iterations } @@ -99,14 +101,15 @@ eps = quadPrecision param -- Requred precision maxN = maxIter param -- Maximum allowed number of iterations worker n nPoints q - | n > 5 && d < eps = ret (Just q') - | n >= maxN = ret Nothing - | otherwise = worker (n+1) (nPoints*2) q' + | n > 5 && converged eps q q' = done True + | n >= maxN = done False + | otherwise = worker (n+1) (nPoints*2) q' where q' = nextTrapezoid a b nPoints f q -- New approximation - d = abs (q' - q) / abs q -- Precision estimate - ret = \x -> QuadRes x d n + done True = QuadRes (Just q') q' (estimatePrec q q') n + done False = QuadRes Nothing q' (estimatePrec q q') n + -- | Integration using Simpson rule. It should be more efficient than -- 'quadTrapezoid' if function being integrated have finite fourth -- derivative. @@ -119,14 +122,14 @@ eps = quadPrecision param -- Requred precision maxN = maxIter param -- Maximum allowed number of points for evaluation worker n nPoints s st - | n > 5 && d < eps = ret (Just s') - | n >= maxN = ret Nothing - | otherwise = worker (n+1) (nPoints*2) s' st' + | n > 5 && converged eps s s' = done True + | n >= maxN = done False + | otherwise = worker (n+1) (nPoints*2) s' st' where st' = nextTrapezoid a b nPoints f st s' = (4*st' - st) / 3 - d = abs (s' - s) / abs s - ret = \x -> QuadRes x d n + done True = QuadRes (Just s') s' (estimatePrec s s') n + done False = QuadRes Nothing s' (estimatePrec s s') n -- | Integration using Romberg rule. For sufficiently smooth functions -- (e.g. analytic) it's a fastest of three. @@ -151,11 +154,12 @@ let worker n nPoints st s = do let st' = nextTrapezoid a b nPoints f st s' <- M.write arr 0 st >> nextAppr n st' - let d = abs (s' - s) / abs s + let done True = return $ QuadRes (Just s') s' (estimatePrec s s') n + done False = return $ QuadRes Nothing s' (estimatePrec s s') n case () of - _ | n > 5 && d < eps -> return $ QuadRes (Just s') d n - | n >= maxN -> return $ QuadRes Nothing d n - | otherwise -> worker (n+1) (nPoints*2) st' s' + _ | n > 5 && converged eps s s' -> done True + | n >= maxN -> done False + | otherwise -> worker (n+1) (nPoints*2) st' s' -- Calculate integral worker 1 1 st0 st0 where st0 = trapGuess a b f @@ -182,3 +186,26 @@ sep = (b - a) / fromIntegral n -- Separation between points x0 = a + 0.5 * sep -- Starting point s = U.sum $ U.map f $ U.iterateN n (+sep) x0 -- Sum of all points + + +-- Check for convergence. Convergence to zero is tricky case since we +-- use relative error and checking for convergence to zero requires +-- absolute one. Instead iteration have converged if two successive +-- operations produced zero +converged :: Double -> Double -> Double -> Bool +converged eps q q' + -- Iterations yielded same numbers. + | q == q' = True + -- Both numbers have identical IEEE representation It's not the same + -- as previous clause because of excess precision. And previous + -- clause is required because of negative zero. + | IEEE.identicalIEEE q q' = True + -- Usual check for convergence. + | otherwise = abs (q' - q) < eps * abs q +{-# INLINE converged #-} + +-- Estimate precision +estimatePrec :: Double -> Double -> Double +estimatePrec q q' + | q == 0 && q' == 0 = 0 + | otherwise = abs (q - q') / abs q diff --git a/numeric-tools.cabal b/numeric-tools.cabal --- a/numeric-tools.cabal +++ b/numeric-tools.cabal @@ -1,6 +1,6 @@ Name: numeric-tools -Version: 0.1.0.1 -Cabal-Version: >= 1.6 +Version: 0.2.0.0 +Cabal-Version: >= 1.8 License: BSD3 License-File: LICENSE Author: Aleksey Khudyakov @@ -15,21 +15,43 @@ Package provides function to perform numeric integration and differentiation, function interpolation. . - Changes in 0.1.0.1 + Changes in 0.2.0.0 . - * Fixed bug in quadRomberg which caused crash if algorithm fails to - converge + * Equation solvers now use custom return type. + . + * Function to solve equations using Ridder and Newton methods. + . + * New function to test approximate equality for doubles. + . + * QuadRes contains best approximation achieved even if required + accuracy is not obtained + . + * Improve convergence test when integral converges to + zero. Convergence is still poor + source-repository head type: hg location: https://bitbucket.org/Shimuuar/numeric-tools +test-suite tests + type: exitcode-stdio-1.0 + hs-source-dirs: test + main-is: test.hs + ghc-options: -Wall -threaded -rtsopts + build-depends: + base, + numeric-tools, + HUnit + Library Build-Depends: base >=3 && <5, ieee754 >= 0.7.3, + primitive, vector >= 0.7.0.1 Exposed-modules: Control.Monad.Numeric + Numeric.ApproxEq Numeric.Classes.Indexing Numeric.Tools.Equation Numeric.Tools.Differentiation diff --git a/test/test.hs b/test/test.hs new file mode 100644 --- /dev/null +++ b/test/test.hs @@ -0,0 +1,123 @@ + +import Test.HUnit +import Text.Printf + +import Numeric.Tools.Integration +import Numeric.Tools.Differentiation +import Numeric.Tools.Equation +import Numeric.ApproxEq + + +---------------------------------------------------------------- +-- Integration +---------------------------------------------------------------- + +-- Functions together with indefinite integral and list of ranges +type FunctionInt = ( Double -> Double -- Function + , Double -> Double -- Integral + , [(Double,Double)] -- List of ranges + , String -- Name + ) + +-- Test integrator +integratorTest :: String + -> (QuadParam -> (Double,Double) -> (Double -> Double) -> QuadRes) + -> QuadParam + -> FunctionInt + -> [Test] +integratorTest name quad param (f,f',ranges,fname) = + [ case quadRes $ quad param (a,b) f of + Nothing -> TestCase $ assertFailure (printf "%s: convergence for %s failed (%f,%f)" name fname a b) + Just appr -> + let exact = f' b - f' a + in TestCase $ assertBool (printf "%s: poor convergence for %s (%f,%f) %g instead of %g" + name fname a b appr exact) + (eqRelative (quadPrecision param) exact appr) + | (a,b) <- ranges + ] + + +testIntegration :: [Test] +testIntegration = concat + [ integratorTest "Trapeze" quadTrapezoid defQuad =<< [funBlamg,funExp,funLog] + , integratorTest "Simpson" quadSimpson defQuad =<< [funBlamg,funExp,funLog] + , integratorTest "Romberg" quadRomberg defQuad =<< [funBlamg,funExp,funLog] + ] + where + funBlamg = ( \x -> x^4 * log(x + sqrt (x*x + 1)) + , \x -> 1/5*x^5 * log(x + sqrt (x*x + 1)) - 1/75*sqrt(x*x+1)*(3*x^4 - 4*x^2 + 8) + , [(0,2), (1,3), (-2,3)] + , "x^4·log(x + sqrt(x^2 + 1))" + ) + funExp = ( exp, exp + , [(0,2), (1,3), (-2,3)] + , "exp" + ) + funLog = ( log + , \x -> x * (log x - 1) + , [(1,2), (0.3,3)] + , "log" + ) + + + +---------------------------------------------------------------- +-- +---------------------------------------------------------------- + +type FunctionDiff = ( Double -> Double -- Function + , Double -> Double -- Derivative + , [(Double,Double)] -- Points and delta to evaluate + , String -- Name + ) + +differentiationTest :: String + -> ((Double -> Double) -> Double -> Double -> DiffRes) + -> FunctionDiff + -> [Test] +differentiationTest name diff (f,f',xs,fname) = + [ let DiffRes appr err = diff f h x + exact = f' x + in TestCase $ + assertBool + (printf "%s: poor precision for %s, got %g instead of %g" name fname appr exact) + (eqRelative 1e-13 appr exact) + | (x,h) <- xs + ] + +testDifferentiation :: [Test] +testDifferentiation = concat + [ differentiationTest "richardson" diffRichardson =<< [funSqr,funExp,funLog] + ] + where + funSqr = ( \x -> x*x + , \x -> 2*x + , zip ([-10 .. -1]++[1..10]) (repeat 1) + , "square" + ) + funExp = ( exp, exp + , zip [-10..10] (repeat 1) + , "exp" + ) + funLog = ( log, recip + , map (\x -> (x,x/3)) [0.1,0.2 .. 2] + , "log" + ) + +testEquation :: [Test] +testEquation = + [ TestCase $ assertBool "Bisection" $ ok (pi/2) (solveBisection 0 (1,2) cos) + , TestCase $ assertBool "Ridders" $ ok (pi/2) $ solveRidders 0 (1,2) cos + , TestCase $ assertBool "Newton" $ ok (pi/2) $ solveNewton 0 (1,2) cos sin + ] + where + ok exact (Root x) = within 1 exact x + ok _ _ = False + +main :: IO () +main = do + res <- runTestTT $ TestList $ concat [ testDifferentiation + , testIntegration + , testEquation + ] + print res