numeric-tools 0.1.0.1 → 0.2.0.0

raw patch · 5 files changed

+463/−40 lines, 5 filesdep +HUnitdep +numeric-toolsdep +primitivedep ~basePVP ok

version bump matches the API change (PVP)

Dependencies added: HUnit, numeric-tools, primitive

Dependency ranges changed: base

API changes (from Hackage documentation)

Files

• Numeric/ApproxEq.hs +81/−0
• Numeric/Tools/Equation.hs +191/−21
• Numeric/Tools/Integration.hs +41/−14
• numeric-tools.cabal +27/−5
• test/test.hs +123/−0

@@ -0,0 +1,81 @@+-- |+-- Module    : Numeric.ApproxEq+-- Copyright : (c) 2011 Aleksey Khudyakov, Bryan O'Sullivan+-- License   : BSD3+--+-- Maintainer  : Aleksey Khudyakov <alexey.skladnoy@gmail.com>+-- 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\":+-- <http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm>+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+

@@ -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&#8211;980. 

@@ -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

@@ -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 <alexey.skladnoy@gmail.com>@@ -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

@@ -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