Generalizing Heron’s method

Suppose we are given a positive rational number $y$ . We would like to compute an approximation of the $n$ -th root of $y$ , which we denote by $\sqrt[n]{y}$ . If $n = 2$ , we can use Heron’s method to approximate $\sqrt{y}$ , as we saw a few days ago. But, what if we have that $n > 2$ ? Is it possible to generalize Heron’s method?

Indeed, it is. Let us introduce a function $f_{n,y} : \mathbb{R} \to \mathbb{R}$ , defined by $f_{n,y} (x) = x^n - y$ , where $n \geq 2$ and $y > 0$ . By construction, the positive zero of function $f_{n,y}$ is $\sqrt[n]{y}$ . Do note that if $n$ is odd, then $-\sqrt[n]{y}$ is not a zero of $f_{n,y}$ . Recall that Newton’s method uses the recurrence relation $x_{k+1} = g_{n,y} ( x_k )$ , where function $g_{n,y} : \mathbb{R} \to \mathbb{R}$ is defined by

$g_{n,y} (x) = x - \displaystyle\frac{f_{n,y} (x)}{f_{n,y}' (x)}$

where $f_{n,y}'$ is the first derivative of $f_{n,y}$ . Hence, we obtain

$g_{n,y} (x) = x - \displaystyle\frac{f_{n,y} (x)}{f_{n,y}' (x)} = x - \displaystyle\left(\frac{x^n - y}{n x^{n-1}}\right) = \displaystyle\frac{1}{n} \left( (n-1) x + \frac{y}{x^{n-1}}\right)$

Thus, we have the generalized Heron’s method

$x_{k+1} = \displaystyle\frac{1}{n} \left( (n-1) x_k + \frac{y}{x_k^{n-1}}\right)$

If we make $n = 2$ , we obtain the standard Heron’s method, as we expected. The following Haskell script implements this generalized Heron’s method using arbitrary-precision rational numbers (included in the Data.Ratio library) to represent $y$ and the sequence of approximations $(x_k)_{k \in \mathbb{N}_0}$ :

import Data.Ratio

-- define generalized Heron map
g :: Integer -> Rational -> Rational -> Rational
g n y x = (1 % n) * ((fromIntegral (n-1) * x) + (y / (x^(n-1))))

-- initial approximation
x0 :: Rational
x0 = 1 % 1

-- list of approximations of the n-th root of y
roots :: Integer -> Rational -> [Rational]
roots n y | n  > 1 = iterate (g n y) x0
          | n <= 1 = error "Invalid value of n!!"

We load this script into GHCi. Let us now perform some numerical experiments.

__________

Experiment #1

We would like to compute a rational approximation of $\sqrt{2}$ . We have $n = 2$ and $y = 2$ . Here is a very brief GHCi session:

*Main Data.Ratio> take 6 (roots 2 2)
[1 % 1,3 % 2,17 % 12,577 % 408,665857 % 470832,
 886731088897 % 627013566048]

We obtain the following rational approximation after five iterations

$x_5 = \displaystyle\frac{886731088897}{627013566048} \approx \sqrt{2}$

which is the same approximation I obtained twice before.

__________

Experiment #2

We would like to compute a rational approximation of $\sqrt{5}$ . We have $n = 2$ and $y = 5$ . Here is a very brief GHCi session:

*Main Data.Ratio> take 7 (roots 2 5)
[1 % 1,3 % 1,7 % 3,47 % 21,2207 % 987,4870847 % 2178309,
 23725150497407 % 10610209857723]

We obtain the following rational approximation after six iterations

$x_6 = \displaystyle\frac{23725150497407}{10610209857723} \approx \sqrt{5}$

How good is this approximation? Let us check:

*Main Data.Ratio> (23725150497407 % 10610209857723)^2 - 5
4 % 112576553224922323902744729
*Main Data.Ratio> 4 / 112576553224922323902744729
3.5531377408652764e-26

Not bad! We could now obtain a rational approximation of the famous golden ratio $\varphi := \frac{1 + \sqrt{5}}{2}$ , as follows:

*Main Data.Ratio> let xs = roots 2 5
*Main Data.Ratio> let phis = map (/2) (map (+1) xs)
*Main Data.Ratio> take 7 phis
[1 % 1,2 % 1,5 % 3,34 % 21,1597 % 987,3524578 % 2178309,
 17167680177565 % 10610209857723]

We obtain the following rational approximation after six iterations

$\tilde{\varphi}_6 = \displaystyle\frac{17167680177565}{10610209857723} \approx \varphi$

How good is this approximation? Let us recall that the golden ratio is one of the solutions of the equation $x^2 - x - 1 = 0$ . Hence, we can check how close to zero $\tilde{\varphi}_6^2 - \tilde{\varphi}_6 - 1$ is:

*Main Data.Ratio> let phi6 = 17167680177565 % 10610209857723
*Main Data.Ratio> phi6^2 - phi6 - 1
1 % 112576553224922323902744729
*Main Data.Ratio> 1 / 112576553224922323902744729
8.882844352163191e-27

I would say that is fairly close.

__________

Experiment #3

We now would like to compute a rational approximation of $\sqrt[3]{10}$ . We have $n = 3$ and $y = 10$ . Here is a very brief GHCi session:

*Main Data.Ratio> take 6 (roots 3 10)
[1 % 1,4 % 1,23 % 8,4909 % 2116,55223315303 % 25495981298,
83759169926117983945469262167029 % 38876457805393768546966848104041]

We obtain the following rational approximation after five iterations

$x_5 = \displaystyle\frac{83759169926117983945469262167029}{38876457805393768546966848104041} \approx \sqrt[3]{10}$

Is this approximation any good? Let us check:

*Main Data.Ratio> let x5 = (roots 3 10) !! 5 
*Main Data.Ratio> x5^3
5876207108075025254603649658428809744645894398375582704
87986732898174481543266076676033790365389 % 58757060813
2677736272189042625165443146909261971794220978628909441
43919908404893015241356540921
*Main Data.Ratio> 5876207108075025254603649658428809744
6458943983755827048798673289817448154326607667603379036
5389 / 587570608132677736272189042625165443146909261971
79422097862890944143919908404893015241356540921
10.000852709004354

This is somewhat disappointing. The rational approximation $x_5$ is a ratio of two enormously long integers, but we do not have such a small error. Note that $x_5^3 - 10 \approx 10^{-3}$ . We can do some quick error analysis. Let $f_{3,10} (x) = x^3 - 10$ . The 1st order Taylor approximation of $f_{3,10}$ around $x_5$ evaluated at $x = \sqrt[3]{10}$ is

$f_{3,10} (\sqrt[3]{10}) \approx f_{3,10} (x_5) + f_{3,10}' (x_5) (\sqrt[3]{10} - x_5)$

and, since $f_{3,10} (\sqrt[3]{10}) = 0$ , we obtain the estimate of the error

$x_5 - \sqrt[3]{10} \approx \displaystyle\frac{f_{3,10} (x_5)}{f_{3,10}' (x_5)}$

In other words, the output deviation $f_{3,10} (x_5) - f_{3,10} (\sqrt[3]{10})$ should be divided by the slope $f_{3,10}' (x_5)$ to obtain an estimate of the magnitude of the approximation error $x_5 - \sqrt[3]{10}$ . Let us use GHCi yet once again:

*Main Data.Ratio> let error = (x5^3 - 10) / (3 * x5^2)
*Main Data.Ratio> error
1670089160826306272530773923851043922672595525468316978
5941152245094153072382174540074985393 % 272741620880842
8590518540423731512512814773109630109912907276686930689
12248623218095571481624481
*Main Data.Ratio> 1670089160826306272530773923851043922
6725955254683169785941152245094153072382174540074985393
 / 2727416208808428590518540423731512512814773109630109
91290727668693068912248623218095571481624481
6.123338108179484e-5

This is still quite a huge error, I would say.

__________

Bonus Experiment

Let us go crazy and attempt to compute a rational approximation of $\sqrt[100]{100}$ . We have $n = 100$ and $y = 100$ . We expect the rational numbers to quickly become ratios of astronomically long integers. Therefore, we will not be showing any rational approximations on GHCi. Here is something scary:

*Main Data.Ratio> let x3 = (roots 100 100) !! 3
*Main Data.Ratio> let error = (x3^100 - 100) / (100 * x3^99)
*Main Data.Ratio> let error_float = fromIntegral (numerator error) 
                                  / fromIntegral (denominator error)
*Main Data.Ratio> error_float

After some 20 minutes without an output, I decided to abort this experiment. If you want to show x3, you can, but you will see thousands of digits. Thousands! You have been warned.

About these ads

Tags: Arbitrary-Precision Arithmetic, Data.Ratio, Golden Ratio, Haskell, Heron’s Method, Newton’s Method, Numerical Methods, Rational Approximations, Rational Arithmetic

This entry was posted on November 25, 2012 at 10:53 and is filed under Haskell, Numerical Methods. You can follow any responses to this entry through the RSS 2.0 feed. You can leave a response, or trackback from your own site.

4 Responses to “Generalizing Heron’s method”

Simon Read Says:
November 29, 2012 at 13:54 | Reply
“rational approximation of \sqrt[3]{10}.
a very brief GHCi session:

*Main Data.Ratio> take 6 (roots 3 10)
1 % 1…

Note that x_5^3 – 10 \approx 10^{-3}.”

I see the problem: you started with 1/1 so used the first three iterations getting roughly near the root. If you start with (for example) 5/2, you’ll be much closer so allowing Newton-Raphson to do its best work. You’ll still have a lot of digits, though!

However, you can then use Haskell to get a rational approximation to your (super-multi-digit) approximation, so you might get something of similar accuracy but with fewer digits.
John Baez Says:
November 30, 2012 at 13:45 | Reply
Are 17167680177565 and 10610209857723 Fibonacci numbers? The best rational approximations to the golden ratio are ratios of Fibonacci numbers.
- Rod Carvalho Says:
  November 30, 2012 at 23:06 | Reply
  Prof. Baez,
  
  Thanks for the insightful comment. I had not thought of using the Fibonacci sequence to compute the golden ratio, though now that you mention it, it seems obvious.
  
  Indeed, those are Fibonacci numbers! Here’s the proof:
```
Prelude> let fibs = 1 : 1 : zipWith (+) fibs (tail fibs)
Prelude> take 20 fibs
[1,1,2,3,5,8,13,21,34,55,89,144,233,
 377,610,987,1597,2584,4181,6765]
Prelude> 17167680177565 `elem` fibs
True
Prelude> 10610209857723 `elem` fibs
True
```
Simon Read Says:
December 15, 2012 at 04:35 | Reply
Surely they are fibonacci numbers. If you use rational arithmetic with Newton-Raphson to solve the polynomial $x^2-x-1=0$ , then you get a selected tour of fractions using fibonacci numbers, *PROVIDED* you start with one pair in the first place, like $8/5$ or $5/3$ .

Rod Carvalho

Generalizing Heron’s method

Like this:

4 Responses to “Generalizing Heron’s method”

Leave a Reply Cancel reply

Follow “Rod Carvalho”