docs/Newton-Raphson.rst

   1 `Newton's method <https://en.wikipedia.org/wiki/Newton%27s_method>`__
   2 =====================================================================
   3
   4 Let's use the Newton-Raphson method for finding the root of an equation
   5 to write a function that can compute the square root of a number.
   6
   7 Cf. `"Why Functional Programming Matters" by John
   8 Hughes <https://www.cs.kent.ac.uk/people/staff/dat/miranda/whyfp90.pdf>`__
   9
  10 .. code:: ipython3
  11
  12     from notebook_preamble import J, V, define
  13
  14 A Generator for Approximations
  15 ------------------------------
  16
  17 To make a generator that generates successive approximations let’s start
  18 by assuming an initial approximation and then derive the function that
  19 computes the next approximation:
  20
  21 ::
  22
  23        a F
  24     ---------
  25         a'
  26
  27 A Function to Compute the Next Approximation
  28 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  29
  30 This is the equation for computing the next approximate value of the
  31 square root:
  32
  33 :math:`a_{i+1} = \frac{(a_i+\frac{n}{a_i})}{2}`
  34
  35 ::
  36
  37     a n over / + 2 /
  38     a n a    / + 2 /
  39     a n/a      + 2 /
  40     a+n/a        2 /
  41     (a+n/a)/2
  42
  43 The function we want has the argument ``n`` in it:
  44
  45 ::
  46
  47     F == n over / + 2 /
  48
  49 Make it into a Generator
  50 ~~~~~~~~~~~~~~~~~~~~~~~~
  51
  52 Our generator would be created by:
  53
  54 ::
  55
  56     a [dup F] make_generator
  57
  58 With n as part of the function F, but n is the input to the sqrt
  59 function we’re writing. If we let 1 be the initial approximation:
  60
  61 ::
  62
  63     1 n 1 / + 2 /
  64     1 n/1   + 2 /
  65     1 n     + 2 /
  66     n+1       2 /
  67     (n+1)/2
  68
  69 The generator can be written as:
  70
  71 ::
  72
  73     23 1 swap  [over / + 2 /] cons [dup] swoncat make_generator
  74     1 23       [over / + 2 /] cons [dup] swoncat make_generator
  75     1       [23 over / + 2 /]      [dup] swoncat make_generator
  76     1   [dup 23 over / + 2 /]                    make_generator
  77
  78 .. code:: ipython3
  79
  80     define('gsra 1 swap [over / + 2 /] cons [dup] swoncat make_generator')
  81
  82 .. code:: ipython3
  83
  84     J('23 gsra')
  85
  86
  87 .. parsed-literal::
  88
  89     [1 [dup 23 over / + 2 /] codireco]
  90
  91
  92 Let's drive the generator a few time (with the ``x`` combinator) and
  93 square the approximation to see how well it works...
  94
  95 .. code:: ipython3
  96
  97     J('23 gsra 6 [x popd] times first sqr')
  98
  99
 100 .. parsed-literal::
 101
 102     23.0000000001585
 103
 104
 105 Finding Consecutive Approximations within a Tolerance
 106 -----------------------------------------------------
 107
 108 From `"Why Functional Programming Matters" by John
 109 Hughes <https://www.cs.kent.ac.uk/people/staff/dat/miranda/whyfp90.pdf>`__:
 110
 111     The remainder of a square root finder is a function *within*, which
 112     takes a tolerance and a list of approximations and looks down the
 113     list for two successive approximations that differ by no more than
 114     the given tolerance.
 115
 116 (And note that by “list” he means a lazily-evaluated list.)
 117
 118 Using the *output* ``[a G]`` of the above generator for square root
 119 approximations, and further assuming that the first term a has been
 120 generated already and epsilon ε is handy on the stack...
 121
 122 ::
 123
 124        a [b G] ε within
 125     ---------------------- a b - abs ε <=
 126           b
 127
 128
 129        a [b G] ε within
 130     ---------------------- a b - abs ε >
 131        b [c G] ε within
 132
 133 Predicate
 134 ~~~~~~~~~
 135
 136 ::
 137
 138     a [b G]             ε [first - abs] dip <=
 139     a [b G] first - abs ε                   <=
 140     a b           - abs ε                   <=
 141     a-b             abs ε                   <=
 142     abs(a-b)            ε                   <=
 143     (abs(a-b)<=ε)
 144
 145 .. code:: ipython3
 146
 147     define('_within_P [first - abs] dip <=')
 148
 149 Base-Case
 150 ~~~~~~~~~
 151
 152 ::
 153
 154     a [b G] ε roll< popop first
 155       [b G] ε a     popop first
 156       [b G]               first
 157        b
 158
 159 .. code:: ipython3
 160
 161     define('_within_B roll< popop first')
 162
 163 Recur
 164 ~~~~~
 165
 166 ::
 167
 168     a [b G] ε R0 [within] R1
 169
 170 1. Discard a.
 171 2. Use ``x`` combinator to generate next term from ``G``.
 172 3. Run ``within`` with ``i`` (it is a "tail-recursive" function.)
 173
 174 Pretty straightforward:
 175
 176 ::
 177
 178     a [b G]        ε R0           [within] R1
 179     a [b G]        ε [popd x] dip [within] i
 180     a [b G] popd x ε              [within] i
 181       [b G]      x ε              [within] i
 182     b [c G]        ε              [within] i
 183     b [c G]        ε               within
 184
 185     b [c G] ε within
 186
 187 .. code:: ipython3
 188
 189     define('_within_R [popd x] dip')
 190
 191 Setting up
 192 ~~~~~~~~~~
 193
 194 The recursive function we have defined so far needs a slight preamble:
 195 ``x`` to prime the generator and the epsilon value to use:
 196
 197 ::
 198
 199     [a G] x ε ...
 200     a [b G] ε ...
 201
 202 .. code:: ipython3
 203
 204     define('within x 0.000000001 [_within_P] [_within_B] [_within_R] tailrec')
 205     define('sqrt gsra within')
 206
 207 Try it out...
 208
 209 .. code:: ipython3
 210
 211     J('36 sqrt')
 212
 213
 214 .. parsed-literal::
 215
 216     6.0
 217
 218
 219 .. code:: ipython3
 220
 221     J('23 sqrt')
 222
 223
 224 .. parsed-literal::
 225
 226     4.795831523312719
 227
 228
 229 Check it.
 230
 231 .. code:: ipython3
 232
 233     4.795831523312719**2
 234
 235
 236
 237
 238 .. parsed-literal::
 239
 240     22.999999999999996
 241
 242
 243
 244 .. code:: ipython3
 245
 246     from math import sqrt
 247
 248     sqrt(23)
 249
 250
 251
 252
 253 .. parsed-literal::
 254
 255     4.795831523312719
 256
 257