[Haskell-cafe] Debugging Newton's method for square roots

[Tamas K Papp]

On Sun, Oct 15, 2006 at 12:17:20PM +0000, Vraj Mohan wrote:
> I am new to Haskell and need help in debugging my code.
> 
> I wrote the following function for calculating square roots using Newton's 
> method:
> 
> my_sqrt :: Float -> Float
> my_sqrt x = improve 1 x
>          where improve y x = if abs (y * y - x) < epsilon 
>                                 then y 
>                                 else improve ((y + (x/y))/ 2) x
>                epsilon = 0.00001
> 
> 
> 
> This works for several examples that I tried out but goes into an infinite loop
> for my_sqrt 96. How do I go about debugging this code in GHC or Hugs?
> 
> (The equivalent code is well-behaved on MIT Scheme)

Hi Vraj,

1. First, generally about debugging: I asked this question a month
ago, on this mailing list.  See this thread:

http://www.haskell.org/pipermail/haskell-cafe/2006-September/017858.html

2. Newton's method is not guaranteed to converge.  For examples, and a
nice introduction to uni- and multivariate rootfinding, see Numerical
Methods in Economics, Kenneth L Judd, Chapter 5.

I suggest that you use bisection, it is much more robust.  I have
written a bisection routine which I have been planning to post after a
bit of testing, maybe you will find it useful:

bisection tol reltol f a b | tol < 0 || reltol < 0 =
                               error "bisection needs positive tolerances"
                           | abs fa <= tol = Just a -- a is a root
                           | abs fb <= tol = Just b -- b is a root
                           | fa*fb > 0 = Nothing -- IVT not guaranteed
                           | a > b = bis b a fb fa -- exchange endpoints
                           | otherwise = bis a b fa fb -- standard case
    where
      fa = f a                  -- function evaluated only once at each
      fb = f b                  -- endpoint, store them here
      bis a b fa fb = let m = (a+b)/2
                          fm = f m in
                      if (b-a) <= reltol*(1+abs(a)+abs(b)) || abs fm <= tol
                      then Just m
                      -- implicit assumption: fa * fb > 0
                      else if fm*fa < 0 -- new interval: [a,m]
                           then bis a m fa fm
                           else bis m b fm fb -- new interval: [m,b]

-- uniroot is a specialisation of bisection, with commonly used
-- tolerance values
uniroot = bisection 1e-20 1e-15


Untested code:

mysqrt x = uniroot f 0 x
       where f y = y**2-x

It will give you a Maybe a type value, ie Nothing for negative
numbers.  Comments about my code are of course welcome, I am a newbie
myself.

Best,

Tamas