Thun/docs/Newton-Raphson.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [Newton's method](https://en.wikipedia.org/wiki/Newton%27s_method)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from notebook_preamble import J, V, define"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cf. [\"Why Functional Programming Matters\" by John Hughes](https://www.cs.kent.ac.uk/people/staff/dat/miranda/whyfp90.pdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$a_{i+1} = \\frac{(a_i+\\frac{n}{a_i})}{2}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's define a function that computes the above equation:\n",
    "\n",
    "         n a Q\n",
    "    ---------------\n",
    "       (a+n/a)/2\n",
    "\n",
    "    n a tuck / + 2 /\n",
    "    a n a    / + 2 /\n",
    "    a n/a      + 2 /\n",
    "    a+n/a        2 /\n",
    "    (a+n/a)/2\n",
    "\n",
    "We want it to leave n but replace a, so we execute it with `unary`:\n",
    "\n",
    "    Q == [tuck / + 2 /] unary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "define('Q == [tuck / + 2 /] unary')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And a function to compute the error:\n",
    "\n",
    "    n a sqr - abs\n",
    "    |n-a**2|\n",
    "\n",
    "This should be `nullary` so as to leave both n and a on the stack below the error.\n",
    "\n",
    "    err == [sqr - abs] nullary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "define('err == [sqr - abs] nullary')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can define a recursive program that expects a number `n`, an initial estimate `a`, and an epsilon value `ε`, and that leaves on the stack the square root of `n` to within the precision of the epsilon value.  (Later on we'll refine it to generate the initial estimate and hard-code an epsilon value.)\n",
    "\n",
    "    n a ε square-root\n",
    "    -----------------\n",
    "          √n\n",
    "\n",
    "\n",
    "If we apply the two functions `Q` and `err` defined above we get the next approximation and the error on the stack below the epsilon.\n",
    "\n",
    "    n a ε [Q err] dip\n",
    "    n a Q err ε \n",
    "    n a'  err ε \n",
    "    n a' e    ε\n",
    "\n",
    "Let's define the recursive function from here.  Start with `ifte`; the predicate and the base case behavior are obvious:\n",
    "\n",
    "    n a' e ε [<] [popop popd] [J] ifte\n",
    "\n",
    "Base-case\n",
    "\n",
    "    n a' e ε popop popd\n",
    "    n a'           popd\n",
    "      a'\n",
    "\n",
    "The recursive branch is pretty easy.  Discard the error and recur.\n",
    "\n",
    "    w/ K == [<] [popop popd] [J] ifte\n",
    "\n",
    "    n a' e ε J\n",
    "    n a' e ε popd [Q err] dip [K] i\n",
    "    n a'   ε      [Q err] dip [K] i\n",
    "    n a' Q err ε              [K] i\n",
    "    n a''  e   ε               K\n",
    "\n",
    "This fragment alone is pretty useful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "define('K == [<] [popop popd] [popd [Q err] dip] primrec')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5.000000232305737\n"
     ]
    }
   ],
   "source": [
    "J('25 10 0.001 dup K')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5.000000000000005\n"
     ]
    }
   ],
   "source": [
    "J('25 10 0.000001 dup K')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So now all we need is a way to generate an initial approximation and an epsilon value:\n",
    "\n",
    "    square-root == dup 3 / 0.000001 dup K"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "define('square-root == dup 3 / 0.000001 dup K')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.000000000000007\n"
     ]
    }
   ],
   "source": [
    "J('36 square-root')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2212475.6192184356\n"
     ]
    }
   ],
   "source": [
    "J('4895048365636 square-root')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4895048365636.0"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "2212475.6192184356 * 2212475.6192184356"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}