# numeric-tools / Numeric / Tools / Integration.hs

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211``` ```{-# LANGUAGE BangPatterns #-} {-# LANGUAGE DeriveDataTypeable #-} -- | -- Module : Numeric.Tools.Integration -- Copyright : (c) 2011 Aleksey Khudyakov -- License : BSD3 -- -- Maintainer : Aleksey Khudyakov -- Stability : experimental -- Portability : portable -- -- Funtions for numerical integration. 'quadRomberg' or 'quadSimpson' -- are reasonable choices in most cases. For non-smooth function they -- converge poorly and 'quadTrapezoid' should be used then. -- -- For example this code intergrates exponent from 0 to 1: -- -- >>> let res = quadRomberg defQuad (0, 1) exp -- -- >>> quadRes res -- Integration result -- Just 1.718281828459045 -- -- >>> quadPrecEst res -- Estimate of precision -- 2.5844957590474064e-16 -- -- >>> quadNIter res -- Number of iterations performed -- 6 module Numeric.Tools.Integration ( -- * Integration parameters and results QuadParam(..) , defQuad , QuadRes(..) -- * Integration functions , quadTrapezoid , quadSimpson , quadRomberg ) where import Control.Monad.ST import Data.Data (Data,Typeable) import qualified Data.Vector.Unboxed as U import qualified Data.Vector.Unboxed.Mutable as M import qualified Numeric.IEEE as IEEE ---------------------------------------------------------------- -- Data types ---------------------------------------------------------------- -- | Integration parameters for numerical routines. Note that each -- additional iteration doubles number of function evaluation required -- to compute integral. -- -- Number of iterations is capped at 30. data QuadParam = QuadParam { quadPrecision :: Double -- ^ Relative precision of answer , quadMaxIter :: Int -- ^ Maximum number of iterations } deriving (Show,Eq,Data,Typeable) -- Number of iterations limited to 30 maxIter :: QuadParam -> Int maxIter = min 30 . quadMaxIter -- | Default parameters for integration functions -- -- * Maximum number of iterations = 20 -- -- * Precision is 10⁻⁹ defQuad :: QuadParam defQuad = QuadParam { quadPrecision = 1e-9 , quadMaxIter = 20 } -- | 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 } deriving (Show,Eq,Data,Typeable) ---------------------------------------------------------------- -- Different integration methods ---------------------------------------------------------------- -- | Integration of using trapezoids. This is robust algorithm and -- place and useful for not very smooth. But it is very slow. It -- hundreds times slower then 'quadRomberg' if function is -- sufficiently smooth. quadTrapezoid :: QuadParam -- ^ Parameters -> (Double, Double) -- ^ Integration limits -> (Double -> Double) -- ^ Function to integrate -> QuadRes quadTrapezoid param (a,b) f = worker 1 1 (trapGuess a b f) where eps = quadPrecision param -- Requred precision maxN = maxIter param -- Maximum allowed number of iterations worker n nPoints 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 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. quadSimpson :: QuadParam -- ^ Parameters -> (Double, Double) -- ^ Integration limits -> (Double -> Double) -- ^ Function to integrate -> QuadRes quadSimpson param (a,b) f = worker 1 1 0 (trapGuess a b f) where eps = quadPrecision param -- Requred precision maxN = maxIter param -- Maximum allowed number of points for evaluation worker n nPoints 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 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. quadRomberg :: QuadParam -- ^ Parameters -> (Double, Double) -- ^ Integration limits -> (Double -> Double) -- ^ Function to integrate -> QuadRes quadRomberg param (a,b) f = runST \$ do let eps = quadPrecision param maxN = maxIter param arr <- M.new (maxN + 1) -- Calculate new approximation let nextAppr n = runNextAppr 0 4 where runNextAppr i fac s = do x <- M.read arr i M.write arr i s if i >= n then return s else runNextAppr (i+1) (fac*4) \$ s + (s - x) / (fac - 1) -- Maine loop 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 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 && 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 ---------------------------------------------------------------- -- Helpers ---------------------------------------------------------------- -- Initial guess for trapezoid rule trapGuess :: Double -> Double -> (Double -> Double) -> Double trapGuess !a !b f = 0.5 * (b - a) * (f b + f a) -- Refinement of guess using trapeziod algorithms nextTrapezoid :: Double -- Lower integration limit -> Double -- Upper integration limit -> Int -- Number of additional points -> (Double -> Double) -- Function to integrate -> Double -- Approximation -> Double nextTrapezoid !a !b !n f !q = 0.5 * (q + sep * s) where 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 ```