1Many numerical approximation methods compute infinite sequences of results; each, 2hopefully, more accurate than the previous one. 3 4<https://en.wikipedia.org/wiki/Newton's_method Newton's method> 5to find zeroes of a function is one such algorithm. 6 7> {-# LANGUAGE FlexibleInstances, MultiParamTypeClasses, UndecidableInstances #-} 8 9> import Control.Comonad.Trans.Coiter 10> import Control.Comonad.Env 11> import Control.Applicative 12> import Data.Foldable (toList, find) 13 14> data Function = Function { 15> -- Function to find zeroes of 16> function :: Double -> Double, 17> -- Derivative of the function 18> derivative :: Double -> Double 19> } 20> 21> data Result = Result { 22> -- Estimated zero of the function 23> value :: Double, 24> -- Estimated distance to the actual zero 25> xerror :: Double, 26> -- How far is value from being an actual zero; that is, 27> -- the difference between @0@ and @f value@ 28> ferror :: Double 29> } deriving (Show) 30> 31> data Outlook = Outlook { result :: Result, 32> -- Whether the result improves in future steps 33> progress :: Bool } deriving (Show) 34 35To make our lives easier, we will store the problem at hand using the Env 36environment comonad. 37 38> type Solution a = CoiterT (Env Function) a 39 40Problems consist of a function and its derivative as the environment, and 41an initial value. 42 43> type Problem = Env Function Double 44 45We can express an iterative algorithm using unfold over an initial environment. 46 47> newton :: Problem -> Solution Double 48> newton = unfold (\wd -> 49> let f = asks function wd in 50> let df = asks derivative wd in 51> let x = extract wd in 52> x - f x / df x) 53> 54> 55 56To estimate the error, we look forward one position in the stream. The next value 57will be much more precise than the current one, so we can consider it as the 58actual result. 59 60We know that the exact value of a function at one of it's zeroes is 0. So, 61@ferror@ can be computed exactly as @abs (f a - f 0) == abs (f a)@ 62 63> estimateError :: Solution Double -> Result 64> estimateError s = 65> let a:a':_ = toList s in 66> let f = asks function s in 67> Result { value = a, 68> xerror = abs $ a - a', 69> ferror = abs $ f a 70> } 71 72To get a sense of when the algorithm is making any progress, we can sample the 73future and check if the result improves at all. 74 75> estimateOutlook :: Int -> Solution Result -> Outlook 76> estimateOutlook sampleSize solution = 77> let sample = map ferror $ take sampleSize $ tail $ toList solution in 78> let result = extract solution in 79> Outlook { result = result, 80> progress = ferror result > minimum sample } 81 82To compute the square root of @c@, we solve the equation @x*x - c = 0@. We will 83stop whenever the accuracy of the result doesn't improve in the next 5 steps. 84 85The starting value for our algorithm is @c@ itself. One could compute a better 86estimate, but the algorithm converges fast enough that it's not really worth it. 87 88> squareRoot :: Double -> Maybe Result 89> squareRoot c = let problem = flip env c (Function { function = (\x -> x*x - c), 90> derivative = (\x -> 2*x) }) 91> in 92> fmap result $ find (not . progress) $ 93> newton problem =>> estimateError =>> estimateOutlook 5 94 95This program will output the result together with the error. 96 97> main :: IO () 98> main = putStrLn $ show $ squareRoot 3 99 100