[Haskell-cafe] Debugging Newton's method for square roots
Tamas K Papp
tpapp at Princeton.EDU
Sun Oct 15 11:21:16 EDT 2006
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
More information about the Haskell-Cafe
mailing list