# Newton’s method¶

Let’s use the Newton-Raphson method for finding the root of an equation to write a function that can compute the square root of a number.

from notebook_preamble import J, V, define


## A Generator for Approximations¶

To make a generator that generates successive approximations let’s start by assuming an initial approximation and then derive the function that computes the next approximation:

   a F
---------
a'


### A Function to Compute the Next Approximation¶

This is the equation for computing the next approximate value of the square root:

$$a_{i+1} = \frac{(a_i+\frac{n}{a_i})}{2}$$

a n over / + 2 /
a n a    / + 2 /
a n/a      + 2 /
a+n/a        2 /
(a+n/a)/2


The function we want has the argument n in it:

F == n over / + 2 /


### Make it into a Generator¶

Our generator would be created by:

a [dup F] make_generator


With n as part of the function F, but n is the input to the sqrt function we’re writing. If we let 1 be the initial approximation:

1 n 1 / + 2 /
1 n/1   + 2 /
1 n     + 2 /
n+1       2 /
(n+1)/2


The generator can be written as:

23 1 swap  [over / + 2 /] cons [dup] swoncat make_generator
1 23       [over / + 2 /] cons [dup] swoncat make_generator
1       [23 over / + 2 /]      [dup] swoncat make_generator
1   [dup 23 over / + 2 /]                    make_generator

define('gsra 1 swap [over / + 2 /] cons [dup] swoncat make_generator')

J('23 gsra')

[1 [dup 23 over / + 2 /] codireco]


Let’s drive the generator a few time (with the x combinator) and square the approximation to see how well it works…

J('23 gsra 6 [x popd] times first sqr')

23.0000000001585


## Finding Consecutive Approximations within a Tolerance¶

The remainder of a square root finder is a function within, which takes a tolerance and a list of approximations and looks down the list for two successive approximations that differ by no more than the given tolerance.

(And note that by “list” he means a lazily-evaluated list.)

Using the output [a G] of the above generator for square root approximations, and further assuming that the first term a has been generated already and epsilon ε is handy on the stack…

   a [b G] ε within
---------------------- a b - abs ε <=
b

a [b G] ε within
---------------------- a b - abs ε >
b [c G] ε within


### Predicate¶

a [b G]             ε [first - abs] dip <=
a [b G] first - abs ε                   <=
a b           - abs ε                   <=
a-b             abs ε                   <=
abs(a-b)            ε                   <=
(abs(a-b)<=ε)

define('_within_P [first - abs] dip <=')


### Base-Case¶

a [b G] ε roll< popop first
[b G] ε a     popop first
[b G]               first
b

define('_within_B roll< popop first')


### Recur¶

a [b G] ε R0 [within] R1


2. Use x combinator to generate next term from G.

3. Run within with i (it is a “tail-recursive” function.)

Pretty straightforward:

a [b G]        ε R0           [within] R1
a [b G]        ε [popd x] dip [within] i
a [b G] popd x ε              [within] i
[b G]      x ε              [within] i
b [c G]        ε              [within] i
b [c G]        ε               within

b [c G] ε within

define('_within_R [popd x] dip')


### Setting up¶

The recursive function we have defined so far needs a slight preamble: x to prime the generator and the epsilon value to use:

[a G] x ε ...
a [b G] ε ...

define('within x 0.000000001 [_within_P] [_within_B] [_within_R] tailrec')
define('sqrt gsra within')


Try it out…

J('36 sqrt')

6.0

J('23 sqrt')

4.795831523312719


Check it.

4.795831523312719**2

22.999999999999996

from math import sqrt

sqrt(23)

4.795831523312719