Aleksey Khudyakov avatar Aleksey Khudyakov committed b46d0ff

Implement root finding using Newton-Raphson method

Comments (0)

Files changed (1)

Numeric/Tools/Equation.hs

     -- * Equations solversv
   , solveBisection
   , solveRidders  
+  , solveNewton
   -- * References
   -- $references
   ) where
     !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
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.