docs/Newton-Raphson.ipynb

   1 {
   2  "cells": [
   3   {
   4    "cell_type": "markdown",
   5    "metadata": {},
   6    "source": [
   7     "# [Newton's method](https://en.wikipedia.org/wiki/Newton%27s_method)\n",
   8     "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.\n",
   9     "\n",
  10     "Cf. [\"Why Functional Programming Matters\" by John Hughes](https://www.cs.kent.ac.uk/people/staff/dat/miranda/whyfp90.pdf)"
  11    ]
  12   },
  13   {
  14    "cell_type": "code",
  15    "execution_count": 1,
  16    "metadata": {},
  17    "outputs": [],
  18    "source": [
  19     "from notebook_preamble import J, V, define"
  20    ]
  21   },
  22   {
  23    "cell_type": "markdown",
  24    "metadata": {},
  25    "source": [
  26     "## A Generator for Approximations\n",
  27     "\n",
  28     "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:\n",
  29     "\n",
  30     "       a F\n",
  31     "    ---------\n",
  32     "        a'"
  33    ]
  34   },
  35   {
  36    "cell_type": "markdown",
  37    "metadata": {},
  38    "source": [
  39     "### A Function to Compute the Next Approximation\n",
  40     "\n",
  41     "This is the equation for computing the next approximate value of the square root:\n",
  42     "\n",
  43     "$a_{i+1} = \\frac{(a_i+\\frac{n}{a_i})}{2}$"
  44    ]
  45   },
  46   {
  47    "cell_type": "markdown",
  48    "metadata": {},
  49    "source": [
  50     "    a n over / + 2 /\n",
  51     "    a n a    / + 2 /\n",
  52     "    a n/a      + 2 /\n",
  53     "    a+n/a        2 /\n",
  54     "    (a+n/a)/2\n",
  55     "\n",
  56     "The function we want has the argument `n` in it:\n",
  57     "\n",
  58     "    F == n over / + 2 /"
  59    ]
  60   },
  61   {
  62    "cell_type": "markdown",
  63    "metadata": {},
  64    "source": [
  65     "### Make it into a Generator\n",
  66     "\n",
  67     "Our generator would be created by:\n",
  68     "\n",
  69     "    a [dup F] make_generator\n",
  70     "\n",
  71     "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:\n",
  72     "\n",
  73     "    1 n 1 / + 2 /\n",
  74     "    1 n/1   + 2 /\n",
  75     "    1 n     + 2 /\n",
  76     "    n+1       2 /\n",
  77     "    (n+1)/2\n",
  78     "\n",
  79     "The generator can be written as:\n",
  80     "\n",
  81     "    23 1 swap  [over / + 2 /] cons [dup] swoncat make_generator\n",
  82     "    1 23       [over / + 2 /] cons [dup] swoncat make_generator\n",
  83     "    1       [23 over / + 2 /]      [dup] swoncat make_generator\n",
  84     "    1   [dup 23 over / + 2 /]                    make_generator"
  85    ]
  86   },
  87   {
  88    "cell_type": "code",
  89    "execution_count": 2,
  90    "metadata": {
  91     "scrolled": true
  92    },
  93    "outputs": [],
  94    "source": [
  95     "define('gsra 1 swap [over / + 2 /] cons [dup] swoncat make_generator')"
  96    ]
  97   },
  98   {
  99    "cell_type": "code",
 100    "execution_count": 3,
 101    "metadata": {},
 102    "outputs": [
 103     {
 104      "name": "stdout",
 105      "output_type": "stream",
 106      "text": [
 107       "[1 [dup 23 over / + 2 /] codireco]\n"
 108      ]
 109     }
 110    ],
 111    "source": [
 112     "J('23 gsra')"
 113    ]
 114   },
 115   {
 116    "cell_type": "markdown",
 117    "metadata": {},
 118    "source": [
 119     "Let's drive the generator a few time (with the `x` combinator) and square the approximation to see how well it works..."
 120    ]
 121   },
 122   {
 123    "cell_type": "code",
 124    "execution_count": 4,
 125    "metadata": {},
 126    "outputs": [
 127     {
 128      "name": "stdout",
 129      "output_type": "stream",
 130      "text": [
 131       "23.0000000001585\n"
 132      ]
 133     }
 134    ],
 135    "source": [
 136     "J('23 gsra 6 [x popd] times first sqr')"
 137    ]
 138   },
 139   {
 140    "cell_type": "markdown",
 141    "metadata": {},
 142    "source": [
 143     "## Finding Consecutive Approximations within a Tolerance\n",
 144     "\n",
 145     "From [\"Why Functional Programming Matters\" by John Hughes](https://www.cs.kent.ac.uk/people/staff/dat/miranda/whyfp90.pdf):\n",
 146     "\n",
 147     "\n",
 148     "> 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.\n",
 149     "\n",
 150     "(And note that by “list” he means a lazily-evaluated list.)\n",
 151     "\n",
 152     "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...\n",
 153     "\n",
 154     "       a [b G] ε within\n",
 155     "    ---------------------- a b - abs ε <=\n",
 156     "          b\n",
 157     "\n",
 158     "\n",
 159     "       a [b G] ε within\n",
 160     "    ---------------------- a b - abs ε >\n",
 161     "       b [c G] ε within\n",
 162     "\n"
 163    ]
 164   },
 165   {
 166    "cell_type": "markdown",
 167    "metadata": {},
 168    "source": [
 169     "### Predicate\n",
 170     "\n",
 171     "    a [b G]             ε [first - abs] dip <=\n",
 172     "    a [b G] first - abs ε                   <=\n",
 173     "    a b           - abs ε                   <=\n",
 174     "    a-b             abs ε                   <=\n",
 175     "    abs(a-b)            ε                   <=\n",
 176     "    (abs(a-b)<=ε)"
 177    ]
 178   },
 179   {
 180    "cell_type": "code",
 181    "execution_count": 5,
 182    "metadata": {},
 183    "outputs": [],
 184    "source": [
 185     "define('_within_P [first - abs] dip <=')"
 186    ]
 187   },
 188   {
 189    "cell_type": "markdown",
 190    "metadata": {},
 191    "source": [
 192     "### Base-Case\n",
 193     "\n",
 194     "    a [b G] ε roll< popop first\n",
 195     "      [b G] ε a     popop first\n",
 196     "      [b G]               first\n",
 197     "       b"
 198    ]
 199   },
 200   {
 201    "cell_type": "code",
 202    "execution_count": 6,
 203    "metadata": {},
 204    "outputs": [],
 205    "source": [
 206     "define('_within_B roll< popop first')"
 207    ]
 208   },
 209   {
 210    "cell_type": "markdown",
 211    "metadata": {},
 212    "source": [
 213     "### Recur\n",
 214     "\n",
 215     "    a [b G] ε R0 [within] R1\n",
 216     "\n",
 217     "1. Discard a.\n",
 218     "2. Use `x` combinator to generate next term from `G`.\n",
 219     "3. Run `within` with `i` (it is a \"tail-recursive\" function.)\n",
 220     "\n",
 221     "Pretty straightforward:\n",
 222     "\n",
 223     "    a [b G]        ε R0           [within] R1\n",
 224     "    a [b G]        ε [popd x] dip [within] i\n",
 225     "    a [b G] popd x ε              [within] i\n",
 226     "      [b G]      x ε              [within] i\n",
 227     "    b [c G]        ε              [within] i\n",
 228     "    b [c G]        ε               within\n",
 229     "\n",
 230     "    b [c G] ε within"
 231    ]
 232   },
 233   {
 234    "cell_type": "code",
 235    "execution_count": 7,
 236    "metadata": {},
 237    "outputs": [],
 238    "source": [
 239     "define('_within_R [popd x] dip')"
 240    ]
 241   },
 242   {
 243    "cell_type": "markdown",
 244    "metadata": {},
 245    "source": [
 246     "### Setting up\n",
 247     "\n",
 248     "The recursive function we have defined so far needs a slight preamble: `x` to prime the generator and the epsilon value to use:\n",
 249     "\n",
 250     "    [a G] x ε ...\n",
 251     "    a [b G] ε ..."
 252    ]
 253   },
 254   {
 255    "cell_type": "code",
 256    "execution_count": 8,
 257    "metadata": {},
 258    "outputs": [],
 259    "source": [
 260     "define('within x 0.000000001 [_within_P] [_within_B] [_within_R] tailrec')\n",
 261     "define('sqrt gsra within')"
 262    ]
 263   },
 264   {
 265    "cell_type": "markdown",
 266    "metadata": {},
 267    "source": [
 268     "Try it out..."
 269    ]
 270   },
 271   {
 272    "cell_type": "code",
 273    "execution_count": 9,
 274    "metadata": {
 275     "scrolled": true
 276    },
 277    "outputs": [
 278     {
 279      "name": "stdout",
 280      "output_type": "stream",
 281      "text": [
 282       "6.0\n"
 283      ]
 284     }
 285    ],
 286    "source": [
 287     "J('36 sqrt')"
 288    ]
 289   },
 290   {
 291    "cell_type": "code",
 292    "execution_count": 10,
 293    "metadata": {
 294     "scrolled": true
 295    },
 296    "outputs": [
 297     {
 298      "name": "stdout",
 299      "output_type": "stream",
 300      "text": [
 301       "4.795831523312719\n"
 302      ]
 303     }
 304    ],
 305    "source": [
 306     "J('23 sqrt')"
 307    ]
 308   },
 309   {
 310    "cell_type": "markdown",
 311    "metadata": {},
 312    "source": [
 313     "Check it."
 314    ]
 315   },
 316   {
 317    "cell_type": "code",
 318    "execution_count": 11,
 319    "metadata": {
 320     "scrolled": true
 321    },
 322    "outputs": [
 323     {
 324      "data": {
 325       "text/plain": [
 326        "22.999999999999996"
 327       ]
 328      },
 329      "execution_count": 11,
 330      "metadata": {},
 331      "output_type": "execute_result"
 332     }
 333    ],
 334    "source": [
 335     "4.795831523312719**2"
 336    ]
 337   },
 338   {
 339    "cell_type": "code",
 340    "execution_count": 12,
 341    "metadata": {},
 342    "outputs": [
 343     {
 344      "data": {
 345       "text/plain": [
 346        "4.795831523312719"
 347       ]
 348      },
 349      "execution_count": 12,
 350      "metadata": {},
 351      "output_type": "execute_result"
 352     }
 353    ],
 354    "source": [
 355     "from math import sqrt\n",
 356     "\n",
 357     "sqrt(23)"
 358    ]
 359   }
 360  ],
 361  "metadata": {
 362   "kernelspec": {
 363    "display_name": "Python 2",
 364    "language": "python",
 365    "name": "python2"
 366   },
 367   "language_info": {
 368    "codemirror_mode": {
 369     "name": "ipython",
 370     "version": 3
 371    },
 372    "file_extension": ".py",
 373    "mimetype": "text/x-python",
 374    "name": "python",
 375    "nbconvert_exporter": "python",
 376    "pygments_lexer": "ipython3",
 377    "version": "3.8.3"
 378   }
 379  },
 380  "nbformat": 4,
 381  "nbformat_minor": 2
 382 }