docs/Newton-Raphson.md

   1 # [Newton's method](https://en.wikipedia.org/wiki/Newton%27s_method)
   2 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.
   3
   4 Cf. ["Why Functional Programming Matters" by John Hughes](https://www.cs.kent.ac.uk/people/staff/dat/miranda/whyfp90.pdf)
   5
   6
   7 ```python
   8 from notebook_preamble import J, V, define
   9 ```
  10
  11 ## A Generator for Approximations
  12
  13 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:
  14
  15        a F
  16     ---------
  17         a'
  18
  19 ### A Function to Compute the Next Approximation
  20
  21 This is the equation for computing the next approximate value of the square root:
  22
  23 $a_{i+1} = \frac{(a_i+\frac{n}{a_i})}{2}$
  24
  25     a n over / + 2 /
  26     a n a    / + 2 /
  27     a n/a      + 2 /
  28     a+n/a        2 /
  29     (a+n/a)/2
  30
  31 The function we want has the argument `n` in it:
  32
  33     F == n over / + 2 /
  34
  35 ### Make it into a Generator
  36
  37 Our generator would be created by:
  38
  39     a [dup F] make_generator
  40
  41 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:
  42
  43     1 n 1 / + 2 /
  44     1 n/1   + 2 /
  45     1 n     + 2 /
  46     n+1       2 /
  47     (n+1)/2
  48
  49 The generator can be written as:
  50
  51     23 1 swap  [over / + 2 /] cons [dup] swoncat make_generator
  52     1 23       [over / + 2 /] cons [dup] swoncat make_generator
  53     1       [23 over / + 2 /]      [dup] swoncat make_generator
  54     1   [dup 23 over / + 2 /]                    make_generator
  55
  56
  57 ```python
  58 define('gsra 1 swap [over / + 2 /] cons [dup] swoncat make_generator')
  59 ```
  60
  61
  62 ```python
  63 J('23 gsra')
  64 ```
  65
  66     [1 [dup 23 over / + 2 /] codireco]
  67
  68
  69 Let's drive the generator a few time (with the `x` combinator) and square the approximation to see how well it works...
  70
  71
  72 ```python
  73 J('23 gsra 6 [x popd] times first sqr')
  74 ```
  75
  76     23.0000000001585
  77
  78
  79 ## Finding Consecutive Approximations within a Tolerance
  80
  81 From ["Why Functional Programming Matters" by John Hughes](https://www.cs.kent.ac.uk/people/staff/dat/miranda/whyfp90.pdf):
  82
  83
  84 > 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.
  85
  86 (And note that by “list” he means a lazily-evaluated list.)
  87
  88 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...
  89
  90        a [b G] ε within
  91     ---------------------- a b - abs ε <=
  92           b
  93
  94
  95        a [b G] ε within
  96     ---------------------- a b - abs ε >
  97        b [c G] ε within
  98
  99
 100
 101 ### Predicate
 102
 103     a [b G]             ε [first - abs] dip <=
 104     a [b G] first - abs ε                   <=
 105     a b           - abs ε                   <=
 106     a-b             abs ε                   <=
 107     abs(a-b)            ε                   <=
 108     (abs(a-b)<=ε)
 109
 110
 111 ```python
 112 define('_within_P [first - abs] dip <=')
 113 ```
 114
 115 ### Base-Case
 116
 117     a [b G] ε roll< popop first
 118       [b G] ε a     popop first
 119       [b G]               first
 120        b
 121
 122
 123 ```python
 124 define('_within_B roll< popop first')
 125 ```
 126
 127 ### Recur
 128
 129     a [b G] ε R0 [within] R1
 130
 131 1. Discard a.
 132 2. Use `x` combinator to generate next term from `G`.
 133 3. Run `within` with `i` (it is a "tail-recursive" function.)
 134
 135 Pretty straightforward:
 136
 137     a [b G]        ε R0           [within] R1
 138     a [b G]        ε [popd x] dip [within] i
 139     a [b G] popd x ε              [within] i
 140       [b G]      x ε              [within] i
 141     b [c G]        ε              [within] i
 142     b [c G]        ε               within
 143
 144     b [c G] ε within
 145
 146
 147 ```python
 148 define('_within_R [popd x] dip')
 149 ```
 150
 151 ### Setting up
 152
 153 The recursive function we have defined so far needs a slight preamble: `x` to prime the generator and the epsilon value to use:
 154
 155     [a G] x ε ...
 156     a [b G] ε ...
 157
 158
 159 ```python
 160 define('within x 0.000000001 [_within_P] [_within_B] [_within_R] tailrec')
 161 define('sqrt gsra within')
 162 ```
 163
 164 Try it out...
 165
 166
 167 ```python
 168 J('36 sqrt')
 169 ```
 170
 171     6.0
 172
 173
 174
 175 ```python
 176 J('23 sqrt')
 177 ```
 178
 179     4.795831523312719
 180
 181
 182 Check it.
 183
 184
 185 ```python
 186 4.795831523312719**2
 187 ```
 188
 189
 190
 191
 192     22.999999999999996
 193
 194
 195
 196
 197 ```python
 198 from math import sqrt
 199
 200 sqrt(23)
 201 ```
 202
 203
 204
 205
 206     4.795831523312719
 207
 208