lectures/09-Secant-method.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Secant method\n",
    "\n",
    "The secant method fall between bisection and Newton methods as a way to solve 1d nonlinear equations $f(x)=0$. The idea is best explained with a picture.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"600\" height=\"400\" viewBox=\"0 0 600 400\">\n",
       "<defs>\n",
       "  <clipPath id=\"clip00\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"600\" height=\"400\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip00)\" points=\"\n",
       "0,400 600,400 600,0 0,0 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip01\">\n",
       "    <rect x=\"120\" y=\"0\" width=\"421\" height=\"400\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip00)\" points=\"\n",
       "45.0072,360.065 580.315,360.065 580.315,11.811 45.0072,11.811 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip02\">\n",
       "    <rect x=\"45\" y=\"11\" width=\"536\" height=\"349\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  45.0072,354.842 45.0072,17.0348 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  178.834,354.842 178.834,17.0348 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  312.661,354.842 312.661,17.0348 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  446.488,354.842 446.488,17.0348 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  580.315,354.842 580.315,17.0348 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  53.0368,360.065 572.285,360.065 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  53.0368,243.981 572.285,243.981 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  53.0368,127.896 572.285,127.896 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  53.0368,11.811 572.285,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  45.0072,360.065 580.315,360.065 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  45.0072,360.065 45.0072,354.842 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  178.834,360.065 178.834,354.842 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  312.661,360.065 312.661,354.842 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  446.488,360.065 446.488,354.842 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  580.315,360.065 580.315,354.842 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  45.0072,360.065 45.0072,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  45.0072,360.065 53.0368,360.065 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  45.0072,243.981 53.0368,243.981 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  45.0072,127.896 53.0368,127.896 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  45.0072,11.811 53.0368,11.811 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 45.0072, 373.865)\" x=\"45.0072\" y=\"373.865\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 178.834, 373.865)\" x=\"178.834\" y=\"373.865\">1</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 312.661, 373.865)\" x=\"312.661\" y=\"373.865\">2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 446.488, 373.865)\" x=\"446.488\" y=\"373.865\">3</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 580.315, 373.865)\" x=\"580.315\" y=\"373.865\">4</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 39.0072, 364.565)\" x=\"39.0072\" y=\"364.565\">-5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 39.0072, 248.481)\" x=\"39.0072\" y=\"248.481\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 39.0072, 132.396)\" x=\"39.0072\" y=\"132.396\">5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 39.0072, 16.311)\" x=\"39.0072\" y=\"16.311\">10</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:16; text-anchor:middle;\" transform=\"rotate(0, 312.661, 397.6)\" x=\"312.661\" y=\"397.6\">x</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:16; text-anchor:middle;\" transform=\"rotate(-90, 14.4, 185.938)\" x=\"14.4\" y=\"185.938\">y</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#0000ff; stroke-width:2; stroke-opacity:1; fill:none\" points=\"\n",
       "  111.921,331.044 120.114,329.536 128.308,327.853 136.501,325.997 144.695,323.966 152.888,321.761 161.082,319.383 169.275,316.83 177.469,314.103 185.662,311.202 \n",
       "  193.856,308.127 202.049,304.878 210.242,301.455 218.436,297.858 226.629,294.087 234.823,290.141 243.016,286.022 251.21,281.729 259.403,277.261 267.597,272.62 \n",
       "  275.79,267.804 283.984,262.815 292.177,257.651 300.371,252.313 308.564,246.802 316.758,241.116 324.951,235.256 333.145,229.222 341.338,223.014 349.532,216.632 \n",
       "  357.725,210.076 365.919,203.346 374.112,196.442 382.306,189.364 390.499,182.111 398.693,174.685 406.886,167.085 415.08,159.31 423.273,151.362 431.467,143.239 \n",
       "  439.66,134.943 447.854,126.472 456.047,117.827 464.241,109.008 472.434,100.016 480.628,90.8488 488.821,81.5078 497.015,71.9928 505.208,62.3038 513.401,52.4407 \n",
       "  \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  111.921,243.981 120.114,243.981 128.308,243.981 136.501,243.981 144.695,243.981 152.888,243.981 161.082,243.981 169.275,243.981 177.469,243.981 185.662,243.981 \n",
       "  193.856,243.981 202.049,243.981 210.242,243.981 218.436,243.981 226.629,243.981 234.823,243.981 243.016,243.981 251.21,243.981 259.403,243.981 267.597,243.981 \n",
       "  275.79,243.981 283.984,243.981 292.177,243.981 300.371,243.981 308.564,243.981 316.758,243.981 324.951,243.981 333.145,243.981 341.338,243.981 349.532,243.981 \n",
       "  357.725,243.981 365.919,243.981 374.112,243.981 382.306,243.981 390.499,243.981 398.693,243.981 406.886,243.981 415.08,243.981 423.273,243.981 431.467,243.981 \n",
       "  439.66,243.981 447.854,243.981 456.047,243.981 464.241,243.981 472.434,243.981 480.628,243.981 488.821,243.981 497.015,243.981 505.208,243.981 513.401,243.981 \n",
       "  \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#0000ff; stroke-width:1; stroke-opacity:1; fill:none\" stroke-dasharray=\"8, 5\" points=\"\n",
       "  178.834,313.631 446.488,127.896 \n",
       "  \"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"178.834\" cy=\"313.631\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#0000ff; stroke:none; fill-opacity:1\" cx=\"178.834\" cy=\"313.631\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"446.488\" cy=\"127.896\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#0000ff; stroke:none; fill-opacity:1\" cx=\"446.488\" cy=\"127.896\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"178.834\" cy=\"243.981\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"178.834\" cy=\"243.981\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"446.488\" cy=\"243.981\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"446.488\" cy=\"243.981\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"279.204\" cy=\"243.981\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"279.204\" cy=\"243.981\" r=\"3\"/>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:21; text-anchor:middle;\" transform=\"rotate(0, 178.834, 228.264)\" x=\"178.834\" y=\"228.264\">a</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:21; text-anchor:middle;\" transform=\"rotate(0, 446.488, 228.264)\" x=\"446.488\" y=\"228.264\">b</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:21; text-anchor:middle;\" transform=\"rotate(0, 279.204, 228.264)\" x=\"279.204\" y=\"228.264\">c</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:21; text-anchor:middle;\" transform=\"rotate(0, 178.834, 344.348)\" x=\"178.834\" y=\"344.348\">f(a)</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:21; text-anchor:middle;\" transform=\"rotate(0, 446.488, 105.214)\" x=\"446.488\" y=\"105.214\">f(b)</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:21; text-anchor:middle;\" transform=\"rotate(0, 526.784, 42.528)\" x=\"526.784\" y=\"42.528\">f(x)</text>\n",
       "</g>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using Plots\n",
    "gr()\n",
    "\n",
    "x = linspace(0.5,3.5)\n",
    "f(x) = x.^2 - 4\n",
    "\n",
    "a = 1\n",
    "b = 3\n",
    "fa = f(a)\n",
    "fb = f(b)\n",
    "\n",
    "c = a - fa*(b - a)/(fb - fa)\n",
    "\n",
    "δ = 1 # a small displacement for ploting labels\n",
    "\n",
    "plot(x, f(x), color=\"blue\",  linestyle=:solid, label=\"\", linewidth=2)\n",
    "plot!(x, 0*x, color=\"black\", linestyle=:solid, label=\"\")\n",
    "plot!([a, b], [fa, fb], color=\"blue\", linestyle=:dash, label=\"\")\n",
    "\n",
    "scatter!([a, b],    [fa, fb],  marker=:circ, color=\"blue\",  label=\"\")\n",
    "scatter!([a, b, c], [0, 0, 0], marker=:circ, color=\"black\", label=\"\")\n",
    "\n",
    "annotate!(a, δ, text(\"a\"))\n",
    "annotate!(b, δ, text(\"b\"))\n",
    "annotate!(c, δ, text(\"c\"))\n",
    "annotate!(a, f(a)-δ, text(\"f(a)\"))\n",
    "annotate!(b, f(b)+1.3δ, text(\"f(b)\"))\n",
    "annotate!(3.6, 9, text(\"f(x)\"))\n",
    "plot!(xlabel=\"x\", ylabel=\"y\", ylims=(-5, 10), xlims=(0,4))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We picture the secant method beginning $a,b$ bracketing a root, like the bisection method. Only instead of taking the midpoint as the next iterate, we draw the secant to the curve (the line that connects the two points $(a, f(a))$ and $(b, f(b))$), and we choose the point $c$ where the secant crosses the $y$ axis. \n",
    "\n",
    "Let's derive a formula for the value of $c$. Think of the secant line $\\hat{f}(x)$ as a linear approximation to $f(x)$ that starts at point $(a, f(a))$ and moves upward from there with slope $(f(b) - f(a))/(b-a)$. The formula for such a line is\n",
    "\n",
    "\\begin{equation*}\n",
    "\\hat{f}(x) = f(a) + (x-a) \\, \\frac{f(b) - f(a)}{b-a}\n",
    "\\end{equation*}\n",
    "\n",
    "At the intersection point with $y=0$ we have $\\hat{f}(c) = 0$. Plug that in to the above\n",
    "\n",
    "\\begin{equation*}\n",
    "0 = f(a) + (c-a) \\, \\frac{f(b) - f(a)}{b-a}\n",
    "\\end{equation*}\n",
    "\n",
    "and solve for $c$ to get\n",
    "\n",
    "\\begin{equation*}\n",
    "c = a - f(a) \\, \\frac{b-a}{f(b) - f(a)}\n",
    "\\end{equation*}\n",
    "\n",
    "Unlike the bisection algorithm, $c$ is not guaranteed to lie between $a$ and $b$, so we need some other way of choosing which of $a$ and $b$ to replace for the iteration. Standard practice is to replace whichever of $|f(a)|$ and $|f(b)|$ is larger. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pseudocode for secant algorithm\n",
    "\n",
    "```\n",
    "secantsearch(f, a, b)\n",
    "    start with a,b preferably bracketing a root of f(x)\n",
    "    while approximate solution is still not good enough\n",
    "       compute c = a - f(a) (b-a)/(f(b) - f(a))    \n",
    "       if |f(c)| is not smaller than the lesser of |f(a)| and |f(b)|\n",
    "          then algorithm isn't converging, so exit\n",
    "       end    \n",
    "       if |f(a)| > |f(b)|\n",
    "          let a = c\n",
    "       else\n",
    "          let b = c\n",
    "       end      \n",
    "    end      \n",
    "    return c\n",
    "end\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Julia implementation of secant algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "secantsearch (generic function with 2 methods)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function secantsearch(f, a::AbstractFloat, b::AbstractFloat, diagnostics=false) \n",
    "# return solution x of f(x) == 0, with a,b serving as initial guesses\n",
    "    \n",
    "    c = (a+b)/2         # c will take on bigger float type of a,b\n",
    "    fa = f(a)           # cache values of f(a), f(b), f(c)      \n",
    "    fb = f(b)\n",
    "    fc = f(c)\n",
    "    ϵ = eps(typeof(c))  # machine epsilon\n",
    "    n = 1               # iteration count\n",
    "    N = 1024            # max iterations\n",
    "    \n",
    "    # stop when f(c) is small, |b-a| is small relative to a and b\n",
    "    # or when we reach maximium number of iterations\n",
    "    while abs(fc) > 10ϵ && abs(b-a)/(abs(a)+abs(b)) > 10ϵ && n < N\n",
    "    \n",
    "        # compute new values for c, fc by secant approximation\n",
    "        c = a - fa*(b-a)/(fb-fa)\n",
    "        fc = f(c)\n",
    "        \n",
    "        # print diagnostics if they're turned on\n",
    "        diagnostics && println(\"n = $n, c = $c, f(c) = $(f(c))\")\n",
    "        \n",
    "        # error and exit if f(c) is bigger than both f(a) and f(b)\n",
    "        if abs(fc) > abs(fa) && abs(fc) > abs(fb)\n",
    "            println(\"secantsearch error: f(c) is increasing\")\n",
    "            return c\n",
    "        end\n",
    "        \n",
    "        # replace whichever of a, b has bigger f(a), f(b)\n",
    "        if abs(fa) > abs(fb)\n",
    "            a,fa = c,fc\n",
    "        else\n",
    "            b,fb = c,fc\n",
    "        end      \n",
    "        \n",
    "        # increment iteration counter\n",
    "        n += 1\n",
    "    end\n",
    "    \n",
    "    # return c as best estimate of root\n",
    "    c\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Verify algorithm on a simple test case\n",
    "\n",
    "Let's try out `secantsearch` on a simple function whose root we know, say, $f(x) = x^2 - 4$, which has a root $x=2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "f (generic function with 1 method)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f(x) = x^2 - 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n = 1, c = 1.6, f(c) = -1.4399999999999995\n",
      "n = 2, c = 2.1538461538461537, f(c) = 0.6390532544378695\n",
      "n = 3, c = 1.9836065573770492, f(c) = -0.06530502553077122\n",
      "n = 4, c = 1.9993904297470282, f(c) = -0.00243790943599409\n",
      "n = 5, c = 2.000002508903471, f(c) = 1.003562017931614e-5\n",
      "n = 6, c = 1.9999999996176037, f(c) = -1.5295853472707677e-9\n",
      "n = 7, c = 1.9999999999999998, f(c) = -8.881784197001252e-16\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "1.9999999999999998"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = secantsearch(f,1.0,4.0,true)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Oh, sweet! It converges from order-1 error to nearly double-precision accuracy in only seven iterations. Should we be concerned that it returns 1.9999999999999998 instead of 2.0 exactly? No, in general, we can't expect to find exact answers with floating-point computations. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting convergence\n",
    "\n",
    "HW1 asks you to modify your `bisectsearch` and `newtonsearch` codes to plot the error versus iteration number. I'll illustrate here how to do that with `secantsearch`.\n",
    "\n",
    "Main ideas of the revision:\n",
    "  - Add a vector `x` to store the value of the best guess $x_n$ as you iterate $n=1,2,\\ldots$.\n",
    "  - Return both the final answer `r` (held in variable `c`) and the vector $x$ of iterates.\n",
    "  - After the function returns, compute the error $e_n = |x_n - r|$ and plot it logarithmically. \n",
    "  - Remove the last value of $x$ before plotting so you don't take the log of zero.\n",
    "  - Take some care to get a sensible value for the first iterate $x_1$, even if the search fails at $n=1$. \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "secantsearch (generic function with 2 methods)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function secantsearch(f, a::AbstractFloat, b::AbstractFloat, diagnostics=false) \n",
    "# return solution x of f(x) == 0, with a,b serving as initial guesses\n",
    "    \n",
    "    fa = f(a)           # cache values of f(a), f(b)     \n",
    "    fb = f(b)\n",
    "    \n",
    "    # set c to better of a,b (and fc, too)\n",
    "    c,fc = abs(fa) < abs(fb) ? (a,fa) : (b,fb)\n",
    "    \n",
    "    floattype = typeof((a+b)/2) # this'll eval to bigger of a,b types\n",
    "    ϵ = eps(floattype)  # machine epsilon\n",
    "    n = 1               # iteration count\n",
    "    N = 1024            # max iterations\n",
    "    \n",
    "    # start vector for storing root iterates xₙ, assign x₁ = c\n",
    "    x = zeros(floattype, 1)\n",
    "    x[1] = c\n",
    "    \n",
    "    # stop when f(c) is small, |b-a| is small relative to a and b\n",
    "    # or when we reach maximium number of iterations\n",
    "    while abs(fc) > 10ϵ && abs(b-a)/(abs(a)+abs(b)) > 10ϵ && n < N\n",
    "    \n",
    "        # compute new values for c, fc by secant approximation\n",
    "        c = a - fa*(b-a)/(fb-fa)\n",
    "        fc = f(c)\n",
    "        \n",
    "        # print diagnostics if they're turned on\n",
    "        diagnostics && println(\"n = $n, c = $c, f(c) = $(f(c))\")\n",
    "        \n",
    "        # error and exit if f(c) is bigger than both f(a) and f(b)\n",
    "        if abs(fc) > abs(fa) && abs(fc) > abs(fb)\n",
    "            println(\"secantsearch error: f(c) is bigger than f(a), f(b)\")\n",
    "            break # break out of while loop, don't store bad value c\n",
    "        end\n",
    "        \n",
    "        # replace whichever of a, b has bigger f(a), f(b)\n",
    "        if abs(fa) > abs(fb)\n",
    "            a,fa = c,fc\n",
    "        else\n",
    "            b,fb = c,fc\n",
    "        end      \n",
    "        \n",
    "       \n",
    "        push!(x,c)   # add another value c to x[n]\n",
    "        n += 1\n",
    "        \n",
    "    end\n",
    "    \n",
    "    # return c and vector x\n",
    "    c, x\n",
    "    \n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Run the search again to get the vector of iterates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n = 1, c = 1.6, f(c) = -1.4399999999999995\n",
      "n = 2, c = 2.1538461538461537, f(c) = 0.6390532544378695\n",
      "n = 3, c = 1.9836065573770492, f(c) = -0.06530502553077122\n",
      "n = 4, c = 1.9993904297470282, f(c) = -0.00243790943599409\n",
      "n = 5, c = 2.000002508903471, f(c) = 1.003562017931614e-5\n",
      "n = 6, c = 1.9999999996176037, f(c) = -1.5295853472707677e-9\n",
      "n = 7, c = 1.9999999999999998, f(c) = -8.881784197001252e-16\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(1.9999999999999998, [1.0, 1.6, 2.15385, 1.98361, 1.99939, 2.0, 2.0, 2.0])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r,x = secantsearch(f,1.0,4.0,true)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot logarithmically,  after removing the last iterate to avoid computing log(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.9999999999999998"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pop!(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"600\" height=\"400\" viewBox=\"0 0 600 400\">\n",
       "<defs>\n",
       "  <clipPath id=\"clip00\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"600\" height=\"400\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip00)\" points=\"\n",
       "0,400 600,400 600,0 0,0 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip01\">\n",
       "    <rect x=\"120\" y=\"0\" width=\"421\" height=\"400\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip00)\" points=\"\n",
       "70.3346,360.065 580.315,360.065 580.315,11.811 70.3346,11.811 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip02\">\n",
       "    <rect x=\"70\" y=\"11\" width=\"511\" height=\"349\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  164.954,354.842 164.954,17.0348 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  325.325,354.842 325.325,17.0348 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  485.696,354.842 485.696,17.0348 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  77.9843,360.065 572.665,360.065 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  77.9843,290.415 572.665,290.415 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  77.9843,220.764 572.665,220.764 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  77.9843,151.113 572.665,151.113 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  77.9843,81.4619 572.665,81.4619 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  77.9843,11.811 572.665,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  70.3346,360.065 580.315,360.065 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  164.954,360.065 164.954,354.842 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  325.325,360.065 325.325,354.842 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  485.696,360.065 485.696,354.842 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  70.3346,360.065 70.3346,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  70.3346,360.065 77.9843,360.065 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  70.3346,290.415 77.9843,290.415 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  70.3346,220.764 77.9843,220.764 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  70.3346,151.113 77.9843,151.113 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  70.3346,81.4619 77.9843,81.4619 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  70.3346,11.811 77.9843,11.811 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 164.954, 373.865)\" x=\"164.954\" y=\"373.865\">2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 325.325, 373.865)\" x=\"325.325\" y=\"373.865\">4</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 485.696, 373.865)\" x=\"485.696\" y=\"373.865\">6</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 64.3346, 364.565)\" x=\"64.3346\" y=\"364.565\">10^-10</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 64.3346, 294.915)\" x=\"64.3346\" y=\"294.915\">10^-8</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 64.3346, 225.264)\" x=\"64.3346\" y=\"225.264\">10^-6</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 64.3346, 155.613)\" x=\"64.3346\" y=\"155.613\">10^-4</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 64.3346, 85.9619)\" x=\"64.3346\" y=\"85.9619\">10^-2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 64.3346, 16.311)\" x=\"64.3346\" y=\"16.311\">10^0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:16; text-anchor:middle;\" transform=\"rotate(0, 325.325, 397.6)\" x=\"325.325\" y=\"397.6\">n</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:16; text-anchor:middle;\" transform=\"rotate(-90, 14.4, 185.938)\" x=\"14.4\" y=\"185.938\">|x_n - r|</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#009af9; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  84.768,11.811 164.954,25.6695 245.139,40.1211 325.325,73.9859 405.51,123.774 485.696,206.851 565.882,339.779 \n",
       "  \"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"84.768\" cy=\"11.811\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"84.768\" cy=\"11.811\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"164.954\" cy=\"25.6695\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"164.954\" cy=\"25.6695\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"245.139\" cy=\"40.1211\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"245.139\" cy=\"40.1211\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"325.325\" cy=\"73.9859\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"325.325\" cy=\"73.9859\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"405.51\" cy=\"123.774\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"405.51\" cy=\"123.774\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"485.696\" cy=\"206.851\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"485.696\" cy=\"206.851\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"565.882\" cy=\"339.779\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"565.882\" cy=\"339.779\" r=\"3\"/>\n",
       "<polygon clip-path=\"url(#clip00)\" points=\"\n",
       "459.589,62.931 562.315,62.931 562.315,32.691 459.589,32.691 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  459.589,62.931 562.315,62.931 562.315,32.691 459.589,32.691 459.589,62.931 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#009af9; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  465.589,47.811 501.589,47.811 \n",
       "  \"/>\n",
       "<circle clip-path=\"url(#clip00)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"486.589\" cy=\"47.811\" r=\"6\"/>\n",
       "<circle clip-path=\"url(#clip00)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"486.589\" cy=\"47.811\" r=\"5\"/>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:start;\" transform=\"rotate(0, 507.589, 52.311)\" x=\"507.589\" y=\"52.311\">|x_n - r|</text>\n",
       "</g>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(abs.(x-r), marker=:circ, yscale=:log10, label=\"|x_n - r|\")\n",
    "plot!(xlabel=\"n\", ylabel=\"|x_n - r|\", ylims=(1e-10, 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Does the convergence match the theoretical estimate of convergence?\n",
    "\n",
    "The secant method converges at a rate between bisection's linear convergence and Newton's quadratic convergence. Specifially, for twice differentiable functions and guesses sufficiently close to a non-repeated root, the secant method converges as\n",
    "\n",
    "\\begin{equation*}\n",
    "\\lim_{n\\rightarrow \\infty} \\frac{e_{n+1}}{e_n^\\alpha} = M\n",
    "\\end{equation*}\n",
    "\n",
    "where $\\alpha = (1+\\sqrt{5})/2 \\approx 1.62$ and $M$ is a finite constant related to $f$ and its derivatives at the root. \n",
    "\n",
    "How can we check the actual convergence at finite $n$ is close to the theoretical estimate in the large-$n$ limit? The easiest way is to compute $e_{n+1}/e_n^\\alpha$ for $n=1,2,\\dots$ and see if the values are roughly constant. I'll name the vector of errors `err` to avoid a clash with Julia's built-in Euler constant `e`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7-element Array{Float64,1}:\n",
       " 1.0        \n",
       " 0.4        \n",
       " 0.153846   \n",
       " 0.0163934  \n",
       " 0.00060957 \n",
       " 2.5089e-6  \n",
       " 3.82396e-10"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "err = abs.(x-r)  # vector of errors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.618033988749895"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "α = (1+sqrt(5))/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6-element Array{Float64,1}:\n",
       " 0.4     \n",
       " 0.67759 \n",
       " 0.338838\n",
       " 0.471789\n",
       " 0.399415\n",
       " 0.440899"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "err[2:end]./(err[1:end-1]).^α"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks like $e_{n+1}\\,/\\,e_n^\\alpha = M \\approx  0.44$ (taking the last value as the closest to the large-$n$ limit). Let's construct a theoretical error curve from this formula and the value $e_1 = 0.4$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7-element Array{Float64,1}:\n",
       " 1.0        \n",
       " 0.44       \n",
       " 0.116559   \n",
       " 0.0135861  \n",
       " 0.000419503\n",
       " 1.50981e-6 \n",
       " 1.67785e-10"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N = length(err)\n",
    "err_theory = zeros(N)\n",
    "err_theory[1] = 1.0\n",
    "for n=2:N\n",
    "    err_theory[n] = 0.44*err_theory[n-1]^α\n",
    "end\n",
    "err_theory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now add the thoeretical error estimate to the plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"600\" height=\"400\" viewBox=\"0 0 600 400\">\n",
       "<defs>\n",
       "  <clipPath id=\"clip00\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"600\" height=\"400\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip00)\" points=\"\n",
       "0,400 600,400 600,0 0,0 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip01\">\n",
       "    <rect x=\"120\" y=\"0\" width=\"421\" height=\"400\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip00)\" points=\"\n",
       "70.3346,360.065 580.315,360.065 580.315,11.811 70.3346,11.811 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip02\">\n",
       "    <rect x=\"70\" y=\"11\" width=\"511\" height=\"349\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  164.954,354.842 164.954,17.0348 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  325.325,354.842 325.325,17.0348 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  485.696,354.842 485.696,17.0348 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  77.9843,360.065 572.665,360.065 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  77.9843,290.415 572.665,290.415 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  77.9843,220.764 572.665,220.764 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  77.9843,151.113 572.665,151.113 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  77.9843,81.4619 572.665,81.4619 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  77.9843,11.811 572.665,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  70.3346,360.065 580.315,360.065 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  164.954,360.065 164.954,354.842 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  325.325,360.065 325.325,354.842 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  485.696,360.065 485.696,354.842 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  70.3346,360.065 70.3346,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  70.3346,360.065 77.9843,360.065 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  70.3346,290.415 77.9843,290.415 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  70.3346,220.764 77.9843,220.764 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  70.3346,151.113 77.9843,151.113 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  70.3346,81.4619 77.9843,81.4619 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  70.3346,11.811 77.9843,11.811 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 164.954, 373.865)\" x=\"164.954\" y=\"373.865\">2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 325.325, 373.865)\" x=\"325.325\" y=\"373.865\">4</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 485.696, 373.865)\" x=\"485.696\" y=\"373.865\">6</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 64.3346, 364.565)\" x=\"64.3346\" y=\"364.565\">10^-10</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 64.3346, 294.915)\" x=\"64.3346\" y=\"294.915\">10^-8</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 64.3346, 225.264)\" x=\"64.3346\" y=\"225.264\">10^-6</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 64.3346, 155.613)\" x=\"64.3346\" y=\"155.613\">10^-4</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 64.3346, 85.9619)\" x=\"64.3346\" y=\"85.9619\">10^-2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 64.3346, 16.311)\" x=\"64.3346\" y=\"16.311\">10^0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:16; text-anchor:middle;\" transform=\"rotate(0, 325.325, 397.6)\" x=\"325.325\" y=\"397.6\">n</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:16; text-anchor:middle;\" transform=\"rotate(-90, 14.4, 185.938)\" x=\"14.4\" y=\"185.938\">|x_n - r|</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#009af9; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  84.768,11.811 164.954,25.6695 245.139,40.1211 325.325,73.9859 405.51,123.774 485.696,206.851 565.882,339.779 \n",
       "  \"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"84.768\" cy=\"11.811\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"84.768\" cy=\"11.811\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"164.954\" cy=\"25.6695\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"164.954\" cy=\"25.6695\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"245.139\" cy=\"40.1211\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"245.139\" cy=\"40.1211\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"325.325\" cy=\"73.9859\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"325.325\" cy=\"73.9859\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"405.51\" cy=\"123.774\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"405.51\" cy=\"123.774\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"485.696\" cy=\"206.851\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"485.696\" cy=\"206.851\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"565.882\" cy=\"339.779\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip02)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"565.882\" cy=\"339.779\" r=\"3\"/>\n",
       "<polyline clip-path=\"url(#clip02)\" style=\"stroke:#e26f46; stroke-width:1; stroke-opacity:1; fill:none\" stroke-dasharray=\"8, 5\" points=\"\n",
       "  84.768,11.811 164.954,24.2279 245.139,44.3189 325.325,76.8268 405.51,129.426 485.696,214.533 565.882,352.238 \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip00)\" points=\"\n",
       "415.358,78.051 562.315,78.051 562.315,32.691 415.358,32.691 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  415.358,78.051 562.315,78.051 562.315,32.691 415.358,32.691 415.358,78.051 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#009af9; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  421.358,47.811 457.358,47.811 \n",
       "  \"/>\n",
       "<circle clip-path=\"url(#clip00)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"442.358\" cy=\"47.811\" r=\"6\"/>\n",
       "<circle clip-path=\"url(#clip00)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"442.358\" cy=\"47.811\" r=\"5\"/>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:start;\" transform=\"rotate(0, 463.358, 52.311)\" x=\"463.358\" y=\"52.311\">|x_n - r|</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip00)\" style=\"stroke:#e26f46; stroke-width:1; stroke-opacity:1; fill:none\" stroke-dasharray=\"8, 5\" points=\"\n",
       "  421.358,62.931 457.358,62.931 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip00)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:start;\" transform=\"rotate(0, 463.358, 67.431)\" x=\"463.358\" y=\"67.431\">e_n, theorectical</text>\n",
       "</g>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot!(err_theory, yscale=:log10, linestyle=:dash, label=\"e_n, theorectical\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Not bad."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.0",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}