numeric-tools / Numeric / Tools / Equation.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``` ```{-# LANGUAGE BangPatterns #-} {-# LANGUAGE DeriveDataTypeable #-} -- | -- Module : Numeric.Tools.Equation -- Copyright : (c) 2011 Aleksey Khudyakov -- License : BSD3 -- -- Maintainer : Aleksey Khudyakov -- 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 -- | 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. 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 | otherwise = worker 0 lo flo hi fhi where flo = f lo fhi = f hi worker i a fa b fb | within 1 a b = Root a | (b - a) <= eps = Root c | fc == 0 = Root c | i >= (100 :: Int) = SearchFailed | fa * fc < 0 = worker (i+1) a fa c fc | otherwise = worker (i+1) c fc b fb 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 :: Int -- ^ Maximum number of iterations -> Double -- ^ Absolute error tolerance -> (Double,Double) -- ^ Lower and upper bounds for root -> (Double -> Double) -- ^ Function -> (Double -> Double) -- ^ Function's derivative -> Root Double solveNewton nMax 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 > nMax = 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 step dy y | fy < 0 = worker (i+1) dy y (f y) (f' y) y b | otherwise = worker (i+1) dy y (f y) (f' y) 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. ```