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