# 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.

Cf. “Why Functional Programming Matters” by John Hughes

```
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¶

From “Why Functional Programming Matters” by John Hughes:

The remainder of a square root finder is a functionwithin, 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
```

- Discard a.
- Use
`x`

combinator to generate next term from`G`

. - Run
`within`

with`i`

(it is a`primrec`

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] primrec')
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
```