lectures/Lecture04.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Fixed-point iteration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "#### Definition\n",
    "\n",
    "A function $g(x)$ has a fixed point at $x = p$ if $p = g(p)$.<br><br>\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* This is called a fixed point because $g(g(p)) = g(p) = p$, or more generally $g^{(k)}(p) = p$ (the $k$th composition of $g$ with itself)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* If $g(x)$ has a fixed point at $x = p$ then $f(x) = x-g(x)$ has a root at $x = p$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* If $f(x)$ has a zero at $x = p$ then the function $g(x) = x - f(x)$ has a fixed point at $x = p$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Example\n",
    "\n",
    "Find the fixed point(s) of $g(x) = 1 + \\sqrt{x}$, $x \\geq 0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Existence and uniqueness of fixed points\n",
    "#### Theorem 1\n",
    "\n",
    "1. If $g \\in C[a,b]$ and $g(x) \\in [a,b]$ for all $x \\in [a,b]$, then $g$ has at least one fixed point. \n",
    "2. If, in addition, $g'(x)$ exists on $(a,b)$ and there is a positive constant $K < 1$ with\n",
    "$$ |g'(x)| \\leq K, \\quad \\text{for all} \\quad x \\in (a,b),$$\n",
    "then there is exactly one fixed point in $[a,b]$.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Proof\n",
    "\n",
    "__Part 1__: Let $f(x) = x - g(x)$.  Then $f(a) = a - g(a) \\leq 0$ because $g(a) \\in [a,b]$ or $g(a) \\geq a$.  Also, $f(b) = b - g(b) \\geq 0$ because $g(b) \\in [a,b]$ or $g(b) \\leq b$.  Then because $f(x)$ is continuous, the Intermediate Value Theorem guarantees that there is one root of $f$ in $[a,b]$ and therefore $g$ has a fixed point in $[a,b]$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "__Part 2__ By Part 1, we know that $g(x)$ has at least one fixed point $p \\in [a,b]$.  Let us assume there is another fixed point $p^* \\in [a,b]$, $p^* \\neq p$.  Then $f(x) = x - g(x)$ has two roots in $[a,b]$.  Recall, Rolle's Theorem. It states that $f'(c) = 0$ for some point $c$ that lies between $p$ and $p^*$:\n",
    "\n",
    "$$ 0 = f'(c) = 1 - g'(c) \\Longrightarrow g'(c) = 1.$$\n",
    "\n",
    "This violates our assumptions and $p^*$ cannot exist."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Example\n",
    "\n",
    "Using this theorem, show that $g(x) = 1 + \\sqrt{x}$ has one and only one fixed point on $[1,4]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### This theorem gives a sufficient condition for a unique fixed point to exist.  It does not give a necessary condition. \n",
    "\n",
    "Consider\n",
    "\n",
    "$$g(x) = e^{-x}, \\quad x \\in [0,1].$$\n",
    "\n",
    "Because $g$ is decreasing, $g(0) = 1 \\in [0,1]$, and $g(1) > 0$ we are guaranteed a fixed point on $[0,1]$.  To apply the theorem we look at $g'(x) = - e^{-x}$.\n",
    "\n",
    "On the interval $[0,1]$, $|g'(x)| \\leq 1$ and it attains the value one at $x = 0$.  We cannot apply the theorem (need a strict inequality) to state uniqueness of the fixed point."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "But if we plot $g(x)$ and the function $h(x) = x$, it is clear that they cross only once."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfkAAAFwCAYAAAC7CQL0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd0VVX+/vH3TiMEQihJgEBCJ/SO9FCkKtJRsdOtY9eZ\nsY86OhZULCD2NshIV+m9Se+9lxBSIAEC6bnn98fB+fJzEAIkOffePK+1XNFLbu6js4aHc87en20s\ny0JERES8j4/TAURERKRgqORFRES8lEpeRETES6nkRUREvJRKXkRExEup5EVERLxUnkreGFPZGPOh\nMeY3Y0yaMcYyxlTN43sDjTFvG2NOGGPSL/yMmOsJLSIiIleW1yv5msCtQAqw/Co/4wtgJPAi0Bs4\nAcw1xjS5yp8jIiIiV8HkZRiOMcbHsizXhb8fAXwGVLMs6/AV3tcY2AwMsyzrqwuv+QE7gD2WZfW5\nvvgiIiLyZ/J0Jf97wV+DPkA2MOmin5UD/Aj0MMYUu8afKyIiIldQ0Avv6gOHLMtK+8PrO4AA7McA\nIiIiUgAKuuTLYj/H/6Pki35dRERECoCf0wEuxRgzChgF4FO8VPOgchWIKhtEcX9fh5OJiIgUAMsF\nxr7u3rBhw0nLssLy48cWdMmnAFUu8frvV/DJl/g1LMuaAEwAqNuwiVVqyDuknM/m8Z7RDGtXDR8f\nUzBpRURECtPZOJg6CgJKwJAfwRiMMUfy68cX9O36HUA1Y0zQH16vB2QB+6/0A0oU82P2ozHE1A7l\ntV93cd/X60hMzSiIrCIiIoVnzxwY1w6Ob4C6BbPZrKBL/mfAHxj8+wsXttDdBsyzLCszLz+kbIkA\nPrunBa/2a8Cag6fo9f5yFu1OKJjEIiIiBSknE2b/FSbeBiGVYPQyaHonmPy/S53nkjfGDDLGDAKa\nX3ip14XXOl749SrGmBxjzIu/v8eyrE3Y2+feN8aMMMbciL19rhrw0tUENcZwd+sq/PxIe8KCizHs\n6/W8NGM7Gdm5V/NjREREnHNyP3zeFdaMg1b3w/AFEFqrwD7uap7J//SHf/7kwtelQCfAAL787x8c\nhgKvA68BpYEtQE/LsjZebViA2uWDmf5QO96as4cvVx5i9cFkxg5pSnSF4Gv5cSIiIoVj80T49Unw\nK2Y/f4/uVeAfmaeJd05q0aKFtX79+kv+2pI9iTz101bOZmTz9151uLdtVUwB3O4QERG5Zpmpdrlv\nnQRV2sPAz6BUxJ9+uzFmg2VZLfLjoz36FLpO0eHMeawD7WqU4+WfdzLs63WcPJenx/wiIiIF7/hG\nGN8Btv0EnZ+De2detuDzm0eXPEBoyWJ8eV9LXulTn5UHTtHz/eUs2ZPodCwRESnKXC5Y9SF80R1y\ns+G+WdDxGfAp3HkvHl/yYC/Ku7dtVWY+3I5yJQK476t1vPLzDi3KExGRwncuCf49GOY9D7V7wP3L\noUobR6J4Rcn/rk6FUsx4uB33ta3KVysP0+ejFeyMO+t0LBERKSoOLoHx7eDQcrj5XbjtewhyboK7\nV5U8QKC/Ly/3qc/XQ1uSkpZN349XMH7pAXJd7r3AUEREPFhuNix4Bb7tB4GlYeQiaDmiQPa+Xw2v\nK/nfdYoOZ+5jMdxYpzxvzt7NkM9Wcyz5j4fhiYiIXKeUI/BVL1gxBprdA6MWQ4UGTqcCvLjkwZ6U\nN+6uZrwzuDE7487S64PlTNkQi7tvGxQREQ+xY5q9ej5pDwz6CvqMtefQuwmvLnmwF+UNal6Z2Y92\noG7FYJ78aQsP/rCRlPNZTkcTERFPlZUGM/8CP90HYbXtxXUNBjid6n94fcn/LrJsED+OasOzPeuw\nYFcCPd5fpq12IiJy9RJ2wGedYeO30P4JGDobylR1OtUlFZmSB/D1MTzQqQbTHmxHSHF/7vtqHS9M\n3056lrbaiYjIFVgWrPscPusC6Slw9zTo+hL4+jud7E8VqZL/XYNKIfz8SHuGtavGd6uPcPPY5Ww5\ndtrpWCIi4q7SkmHSXfZ42qrt4f6VUKOz06muqEiWPNhb7V68pR4/jGhFWlYuA8atYsz8vWTluJyO\nJiIi7uTIb/biur1zoftrcMdPUDLM6VR5UmRL/nftaoYy97EY+jSOYOzCffT/ZCW74zVAR0SkyHPl\nwtK34Oub7Fvyw+dB20fAx3Oq03OSFqCQIH/eu60J4+9qTvyZDG75cAUfL95PTq6u6kVEiqSzcfBt\nX1j8OjQYBKOXQaVmTqe6aldznrzX69mgAi2rluGFGdt5e+4e5u1M4N3BjakZXtLpaCIiUlj2zIHp\nD0BOJvQbB42HOD657lrpSv4PypUsxsd3NGPskKYcOXWem8cu5/PlB3FpLK6IiHfLyYTZz8LE2yCk\nMoxeCk3u8NiCB5X8JRlj6NM4gnmPxdC+Ziiv/bqL2yes5sip805HExGRgnByH3x+I6wZD60egBEL\nILSW06mum0r+MsJLBfL5vS14e1Ajdp04S8/3l/Pdb4d1VS8i4i0sCzb9AJ92hDPHYcgk6PUm+BVz\nOlm+UMlfgTGGwS0imft4DC2qluGFGTu458u1HD+d7nQ0ERG5HhlnYepImPGgvajugZUQ3dPpVPlK\nJZ9HEaWL8+2wG3i9fwM2Hk2h53vL+M+6YzrsRkTEEx3fCJ/GwPYp0Pl5uGcGlIpwOlW+U8lfBWMM\nd7aqwpxHY6gXUYpnpmxl2NfrSDib4XQ0ERHJC5cLVn0IX3S3z4C/bxZ0fBp8fJ1OViBU8tcgqlwQ\nE0e25sXe9Vh14BTdxixlso6wFRFxb+eS4N+DYd7zULuHfXJclTZOpypQKvlr5ONjGNa+GrMf7UDt\n8sE89dMW7vtqnZ7Vi4i4owOLYXw7OLQcbh4Dt30PQWWdTlXgVPLXqXpYSf4zug0v31KPtYeS6T5m\nKd+vPqIV+CIi7iA3Gxa8DN/1h+JlYNRiaDnco/e+Xw2VfD7w8THc164a8x6PoUlUaZ6fvp0hn63m\n8EntqxcRcUzKYfiqF6x4D5rfCyMXQ/n6TqcqVCr5fBRZNojvh7fiXwMbsjPuLD0/WMbnyw+Sq6t6\nEZHCtX2qfXJc0l4Y/DXc8gEEBDmdqtCp5POZMYbbWkYx/4mOtKthT8sbNH4V+xJSnY4mIuL9stJg\n5iMweSiERduL6+r3dzqVY1TyBaRCiD0t74Pbm3D45HluHmufbJetk+1ERApGwg6Y0Ak2fgftn4Ch\ns6FMFadTOUolX4CMMfRtUon5T3SkW/3yvD13D30/WsmOuDNORxMR8R6WBWs/gwmdIeM03DMdur5k\nnwFfxKnkC0HohZPtxt/VnMTUTPp+tJJ35u4hMyfX6WgiIp4tLRkm3QWznoJqMXD/SqjeyelUbkPn\nyReing0q0Lp6WV79ZRcfLd7P3B3x/GtQI5pFlXE6moiI5znyG0wZAecSoPvr0PpB8NG168X0X6OQ\nlQ4K4N1bG/PV0Jacz8xh4LhVvDxzB+cyc5yOJiLiGVy5sORf8PVN4BcAw+dB24dV8Jeg/yIO6Rwd\nztzHY7i3TVW++e0w3cYsZf7OBKdjiYi4tzPH4Zs+sOSf0HAwjF5mnyAnl6SSd1BwoD8v96nP1Afa\nElLcn5HfrueB7zfowBsRkUvZPcseTRu3CfqNhwEToFiw06ncmkreDTSNKsPPj7Tn6R7RLNydSNd3\nNRpXROS/sjNg1jPw4xAIibSv3psMcTqVR1DJuwl/Xx8e6lyTuY/F0LByCM9P386tn/6mIToiUrSd\n3AdfdIW1n9oL60YsgNCaTqfyGCp5N1MttAQ/jGjFO4Mbsz/pHDeNXc6YeXvIyNZ2OxEpQiwLNv0A\nn3a0n8MPmQQ93wC/Yk4n8ygqeTdkjGFQ88osfKIjvRtFMHbRfm4au5zVB085HU1EpOBlnIWpI2HG\ng/aiugdWQnRPp1N5JJW8GytXshjv3daEb4fdQHaui9snrOavU7ZyJi3b6WgiIgXj+Ab4NMY+YKbL\n83DPDCgV4XQqj6WS9wAxtcOY91hHRneszk8bYrlxzFJ+3hKHZWlhnoh4CZcLVo6FL7qDKweGzoKY\np8HH1+lkHk0l7yGKB/jyt151mflwOyJKB/LIxE0M/XodR0+lOR1NROT6nEuEHwbB/Bcgupd9clxU\na6dTeQWVvIepHxHCtAfb8WLveqw7lEy395YyduE+zcEXEc90YBGMawdHVsLNY+DW76C4Rn3nF5W8\nB/L1MQxrX42FT3aia93yjJm/l17vL2fl/pNORxMRyZvcbJj/EnzXH4LKwsjF0HI4GON0Mq+ikvdg\nFUIC+fjOZnwz7AZyLYs7P1/DXyZuIjFVE/NExI2lHIYve8LK96H5fXbBl6/ndCqvpJL3Ah1rhzH3\nsRgevbEWc7bHc+M7S/lm1WFyNTFPRNzN9ikwvoM95Gbw13DLBxAQ5HQqr6WS9xKB/r483q02cx+P\noUlUaV6auYN+H69ky7HTTkcTEYGs8zDjYZg8DMLq2Ivr6vd3OpXXU8l7mWqhJfh22A18dEdTEs5m\n0O+TlTw/fRtn0rW3XkQcEr8dJnSCTd9Dhyft7XFlqjidqkjwczqA5D9jDL0bRdCxdhhj5u/lm1WH\nmbM9nudurku/JpUwWtgiIoXBsmDd5zD3OSheGu6ZDtU7OZ2qSNGVvBcLDvTnpVvqM/Ph9lQuE8Tj\nk7Yw5LPV7E/UoTciUsDSkmHSXTDrKajeER5YpYJ3gEq+CGhQKYSpD7Tln/0bsutEKr0+WM6/5uzm\nfGaO09FExBsdWQXj28PeudDjn/bhMiVCnU5VJKnkiwgfH8MdraJY+GRH+japxLglB7jx3aXM1Hhc\nEckvrlxY8iZ8fbN9WtyI+dDmIfBR1ThF/+WLmNCSxXhncGOmPNCW0OAA/jJxE7dPWM3u+LNORxMR\nT3bmOHzTB5a8AQ1vhdHLIKKp06mKvDyVvDEm0hgz2Rhzxhhz1hgz1RgTlcf3RhljvjHGHDXGpBtj\n9hpjXjPGlLi+6HI9mlcpw4yH2vPP/g3Zk5DKzWNX8PLMHVqFLyJXb/csGN8O4jZBv/Ew4FMoFux0\nKgHMlW7VGmOCgC1AJvA8YAGvAUFAI8uyzl/mvSWATYA/8DJwFGgJvALMtCzrtisFbNGihbV+/fq8\n/LvINUo5n8W78/fw7zVHKRMUwLM96zCoeWV8fLQKX0QuIzsD5r8Iaz+FCo1g0FcQWtPpVB7PGLPB\nsqwW+fGz8rKFbiRQHYi2LGv/hQBbgX3AaGDMZd7bDqgF9LQsa+6F1xYbY8oCTxljgizL0jFqDitT\nIoDX+jXk9pZRvDxzB89M2coPa4/yjz71aRxZ2ul4IuKOTu6Dn4ZCwjZo/SB0fdl+Di9uJS+36/sA\nq38veADLsg4BK4G+V3hvwIWvfxy7dvrCZ+tS0Y00qBTCT/e34b3bGhN3Op1+n6zkr1O2cupcptPR\nRMRdWJY91ObTGEiNgzv+Az3fUMG7qbyUfH1g+yVe3wFc6USBBdhX/G8ZY+oZY0oaY7oAjwLjL3er\nX5xhjKF/08oserIjI9pXY/KGWDq/s4RvVh0mJ9fldDwRcVLGWZgyAmY8BJWaw/0roXYPp1PJZeSl\n5MsCKZd4PRm47KG/lmVlAO0vfM4OIBVYCPwCPHxVSaVQBQf689zN9ZjzWAcaVbZn4ff+cAVrDp5y\nOpqIOCF2A3zaAXZMgy7Pwz0zoFRFp1PJFRToFjpjTCAwCSgP3A10BJ4GbgM+vsz7Rhlj1htj1icl\nJRVkRLmCmuHBfDf8Bsbd2YzUjBxum7Cah/+9keOn052OJiKFweWClR/Al93tffBDZ0PM0+Dj63Qy\nyYO8LLxL4dJX7H92hX+x4UAnoNZFz/SXGWPOABOMMeMty9ryxzdZljUBmAD26vo8ZJQCZIyhV8OK\ndIoOZ9zSA3y69ADzdyYwOqY693eqQVCAjkAQ8UrnEmHa/XBgIdTtA33GQvHL3sAVN5OXK/kd2M/l\n/6gesPMK720InL540d4Fay98rZuHzxc3UTzAlye61WbRU53oXr8CYxftp8s7S5m2KRaXzq4X8S77\nF8K4dnBkJfR+D279VgXvgfJS8jOB1saY6r+/YIypir09buYV3hsPlDbG/HHjZKsLX4/nLaa4k0ql\ni/PhkKZMvr8N4aWK8fikLQwYt4qNR690Y0dE3F5utr33/fsBEFQORi6GFsNAp1d6pLwMwymBPQwn\nnf8bhvMqEIw9DOfche+rAhwA/mFZ1j8uvFYV2Ipd9q9jD8NpAbwA7AVusCzrsku2NQzHvblcFlM3\nHeetObtJTM2kb5MInu1Zh4jSxZ2OJiJXK/kQTBkOxzdA86H24TIBQU6nKnIKdRiOZVnnL2x7ew/4\nDntv+0Lgsd8L/vdcgC8X3R2wLOuwMaY19rS714BQ4Bj28/bXr1Tw4v58fAyDmlemV4MKjFtygAnL\nDzJ3Rzz3d6zB6JgaFA/Q4hwRj7BtMvzyOGBg8DdQv5/TiSQfXPFK3mm6kvcsx5LTeHP2bn7ddoKK\nIYH8tVcd+jSOwOhWn4h7yjoPs5+xB9xEtoKBn0PpPB1NIgUkP6/kdQqd5KvIskF8fGczJo1qTdkS\nATz642YGjFvF5mN/HHooIo6L3wYTOsGmH6DDU3DfLBW8l1HJS4FoVb0cMx9uz1sDG3EsOZ1+H6/k\niUmbOXFG++tFHGdZsGYCfHajPcXunhlw4wvgq+2w3ka366XAncvM4ePF+/li+SF8fGBEe3t/fcli\n+g1FpNClJcOMh2HPr1CrO/QbByVCnU4lF8nP2/UqeSk0x5LTeHvuHmZuiSO0ZACPda3N7S0j8fPV\nDSWRQnF4JUwdaQ+56faKfXqc1su4HT2TF48UWTaIsUOaMv2hdlQPLcnz07fT4/1lLNyVgLv/YVPE\no7lyYcmb8E1v+7S4EfOhzUMq+CJAJS+FrklkaSaNbs2ndzfHZcHwb9Yz5LPVbIs943Q0Ee9z5jh8\ncwsseQMa3gqjl0FEU6dTSSHR7XpxVHaui3+vOcoHC/eRfD6L/k0r8VSPaCppmI7I9dv9q30sbE4W\n9B4DjW93OpHkgZ7Ji9c5m5HNuCUH+GLFIQCGt6/GA51qUCrQ3+FkIh4oOwPmvwBrJ0DFxjDoKyhX\nw+lUkkcqefFasSlpvDtvL9M2Hbf32d9YiztaReGvxXkieZO0FyYPg4Rt0Poh6PqS/RxePIZKXrze\nttgzvD5rJ6sPJlM9tATP9qpD93rlNTlP5M9Ylj21bvYz4F8c+o2H2t2dTiXXQKvrxes1rBzCxJGt\n+fyeFhgDo7/bwKDxv7HucLLT0UTcT8YZ+2CZmQ9D5RZw/0oVvAC6khcPkJPrYtL6Y3ywYB+JqZnc\nWCecZ3rWIbpCsNPRRJwXuwEmD4UzsdD579D+cfDRwVCeTLfrpUhKz8rly5WHGL/0AOcyc+jftBJP\ndKtN5TI6ClOKIJcLVo2FRa9CcIR9sExUK6dTST5QyUuRlnI+i3FLD/D1qsNgwd1tqvBQ55qULRHg\ndDSRwnEuEaaNhgOLoF5fuGUsFC/tdCrJJyp5ESDudDrvL9jL5A2xlAjwY1RMdYZ3qEZQgGbiixfb\nv9Au+MxU6PkmNL9Pk+u8jEpe5CL7ElJ5a+4e5u9MICy4GH+5sRa3t4zUtjvxLjlZ9q35VWMhrC4M\n/grC6zqdSgqASl7kEjYcSeZfs/ew9nAyVcsF8WT3aG5uWBEfH13liIdLPmSvnj++AVoMgx7/tLfJ\niVdSyYv8CcuyWLwnkX/N3sOehFQaVgrh2Z51aF9LR2mKh9o2GX5+DHx8oM+H9jN48WoqeZEryHVZ\nTN90nDHz93L8dDptqpfjqR7RNK9SxuloInmTdd4ebLPpe4hsZa+eLx3ldCopBCp5kTzKzMnl32uO\n8vHiA5w8l0nn6DCe7B5Ng0ohTkcT+XMnttqjaU/thw5PQqe/ga8WlBYVKnmRq5SWlcM3q44wfukB\nzqRnc1PDCjzetTa1ymugjrgRy4K1n8G856B4WRgwAap3dDqVFDKVvMg1OpuRzefLD/HF8oOkZ+fS\nr0klHutam6hyGqgjDktLhhkPw55foVYP6PcJlNBakqJIJS9ynZLPZ/Hp0gN889thcnItBreI5C83\n1qRiiFYsiwMOr4SpI+0hN93+Aa0f0N73IkwlL5JPEs9m8NHi/UxcexRjDHe1qsIDnWoQFqyjOaUQ\n5ObAsrdh2VtQphoM+hIimjidShymkhfJZ7EpaYxduI/JG2Ip5ufL0HZVGR1Tg5Agf6ejibc6EwtT\nRsLRVdB4CNz0NhTTGhFRyYsUmINJ53hvwT5+3hJHcKAfIztUZ2i7qgQHquwlH+3+FaY/CK4cuHkM\nNL7N6UTiRlTyIgVs14mzjJm/l/k7Ewgp7s/IDtW4t63KXq5TdgbMex7WfQYVm9i358vVcDqVuBmV\nvEgh2RZ7hg8W7mXBrkRKB/kzskN17mlTRWUvVy9pj733PWE7tHkYbnwJ/HRyovwvlbxIIdsae5oP\nFuxj4e7/K/t721alZDENKJErsCzY9B3Mfhb8g6DfOKjd3elU4sZU8iIO2Rp7mvcX7GORyl7yIuOM\nPXd+x1So1tEebhNcwelU4uZU8iIO23LsNB8stMu+TJA/I2Oqc08blb1cJHa9fXv+TCx0eQ7aPQY+\nvk6nEg+gkhdxE5uPneaDBXtZvCfpv2V/b5uqlFDZF10uF6z6ABa9BsERMOgLiLzB6VTiQVTyIm5m\n87HTvL9gL0sulP2omBrc06aKyr6oSU2AaaPh4GKo1w9u+QCKl3Y6lXgYlbyIm9p0NIUPFu77b9mP\n6FCdu9tUoZRW43u//Qtg2v2QmQo934Tm92k0rVwTlbyIm9t4NIWPFu1n0e5EggP9GNq2KkPbVaNM\nCW2Z8jo5WbDoVVg1FsLr2Xvfw+s6nUo8mEpexENsP36GjxbtZ86OeIICfLm7dRWGd6hGeHCg09Ek\nPyQfhMnDIW4jtBgGPf4J/jrkSK6PSl7Ew+xNSOXjxfv5eUsc/r4+DLkhitEdq+vUO0+2bbK9Pc7H\nB/p8CPX6Op1IvIRKXsRDHTp5nnFL9jN143GMgUHNI3mgYw2dZ+9Jss7DrGdg8/cQ2RoGfgalo5xO\nJV5EJS/i4WJT0vh06UEmrTtGrmXRt0kED3aqSc3wkk5Hk8s5sdXe+35qP8Q8BR3/Cr7aQSH5SyUv\n4iUSzmYwYdlBflhzhMwcFzc1rMjDnWtSt2Ipp6PJxSwL1k6wD5cJKmdPrqsW43Qq8VIqeREvc+pc\nJl+sOMS3vx3hXGYOXeuG80CnmjSvUsbpaJKWDDMegj2zoHZP6PsJlCjndCrxYip5ES91Ji2br1Yd\n4utVhzmdls0N1cryQKcadKodhtGe68J3eAVMGQlpJ6HbP6DV/dr7LgVOJS/i5dKycvhx7TE+W36Q\nE2cyqFuxFPd3rM7NDSvi5+vjdDzvl5sDy96CZW9DmWr23veIJk6nkiJCJS9SRGTluJi5JY7xSw+w\nP/EckWWLMyqmBoObVybQX4edFIgzsfbV+9FV0PgOuOltKKYFkVJ4VPIiRYzLZbFgVwKfLDnA5mOn\nCS0ZwNB21birdRVCimtkbr7Z9Yv9/N2VA73fg0a3Op1IiiCVvEgRZVkWqw8mM27pAZbtTSK4mB93\ntq7CsHZVCS+lKXrXLDsD5j0H6z6Hik3s2/PlajidSooolbyIsP34GcYvPcCsbSfw8/VhUPPKjOpQ\nnaqhJZyO5lmS9sBPQyFxB7R5GG58Cfx0xoA4RyUvIv91+OR5Jiw/yOT1seS4XPRqUJFRMdVpHKkj\nTi/LsmDjtzD7WQgoAf3HQ61uTqcSUcmLyP9KPJvBlysP88OaI6Rm5HBDtbKM6lCdLnXC8fHRtq//\nT8YZe+78jqlQraM93Ca4gtOpRACVvIhcRmpGNpPWHeOrlYc5fjqd6mElGNmhOv2bVtKKfIDY9fZo\n2jOx0OU5aPe4fciMiJtQyYvIFWXnupi17QQTlh1kR9xZQksGcE+bqtzdukrRPNfe5YJVH8Ci1yA4\nAgZ9AZE3OJ1K5H+o5EUkzyzL4reDp5iw7CBL9iQR6O/DrS0iGd6+GlXKFZFFeqkJMG0UHFwC9frB\nLR9Aca1ZEPeUnyWv45NEvJwxhrY1QmlbI5Q98al8vvwgE9ce5bvVR+hZvwIjY6rTLMqLZ+TvWwDT\nRttHxN7yATS7V6NppcjI04MoY0ykMWayMeaMMeasMWaqMSbPBygbY+oaY34yxpw0xqQbY/YYYx69\n9tgici2iKwTz9uDGrHi2Cw90rMHK/ScZ8MkqBo1bxdwd8bhc7n1n76rkZNmnxv0wEEqGw6gl0Pw+\nFbwUKVe8XW+MCQK2AJnA84AFvAYEAY0syzp/hfe3ABYBS4AvgTNALaCkZVljrhRQt+tFCs75zBz+\ns/4YX6w4RGxKOlXLBXFf26oMbhFJiWIefKMv+SBMHg5xG6HFcOjxOvgXdzqVSJ4U6jP5C1fcY4Bo\ny7L2X3itGrAPeOZyRW2M8QG2A3ssy+p/LQFV8iIFLyfXxezt8Xy58hCbjp4mONCP21tGcm/bqlQu\nE+R0vKuz9Sf45cKK+T4fQb0+TicSuSqFXfILgUDLstr94fWlAJZldbzMe7sAC4EYy7KWX0tAlbxI\n4dp4NIUvVxxi9vZ4LMuiZ4MKDGtXjeZVyrj3cbeZ52D2M7D5B4hsDQM/h9KRTqcSuWqFvfCuPjDj\nEq/vAAZf4b3tL3wNNMasBpoDKcCPwLOWZaXnNaiIFI5mUWVodkcZ4k6n8+1vR5i49iiztsXTqHII\nw9tXo1eDigT4udm+8hNb7b3vp/ZDzDPQ8Vnw9eDHDSL5JC//Ty2LXcx/lAxcaUluxIWvk4B5QDfg\nLWAE8O8/e5MxZpQxZr0xZn1SUlIeIopIfosoXZy/9qrDb3/rwqv9GnAuM4dHf9xMh7cW8fHi/aSc\nz3I6oj2advV4+PxGyDoH9/5sD7hRwYsABb+F7vc/RHxvWdaLF/5+iTHGF3jTGFPXsqxdf3yTZVkT\ngAlg364v4IwichlBAX7c3boKd94QxdJ9SXy54hBvz93Dh4v20b9pZYa3r0rN8ODCD3b+lH0s7N7Z\nULsn9P31YIsFAAAgAElEQVQESpQr/BwibiwvJZ/Cpa/Y/+wK/2KnLnyd/4fX5wFvAk2A/yl5EXE/\nPj6GztHhdI4OZ29CKl+tPMTUjbFMXHuUmNphDG1blY61wwpnTv6h5TB1JKSdgp7/glajtTVO5BLy\nUvI7sJ/L/1E9YGce3isiXqZ2+WDeGNCIp3vU4d9rjvDtb0cY+vU6qpYL4u42VRnUvDIhxf3z/4Nz\nc2Dpv2DZ2/Z573dMgoqN8/9zRLxEXp7JzwRaG2Oq//6CMaYq0O7Cr13ObOz99T3+8HrPC1/X5Sml\niLilsiUCeLhLLVY824UPhzQltGQxXv1lJ23eWMhz07axJz41/z7s9DH4pjcsewua3AGjlqrgRa4g\nL1voSmAPw0nn/4bhvAoEYw/DOXfh+6oAB4B/WJb1j4ve/xLwAvaCu0VAC+AlYJJlWfddKaC20Il4\nlu3Hz/Dtb4eZsTmOzBwXbaqX4962Vehatzx+vte4Kn/XzzDjYXDlQu8x0OjWfM0s4k4K/YCaCyNs\n38NeHW+w974/ZlnW4Yu+pypwCHjFsqyXL3rdAI8DDwJRwAngG+BVy7Kyr/TZKnkRz5RyPosf1x3j\n+9VHOH46nUqli3Nn6yhubxlF2byegpedbo+mXfc5RDSFQV9C2epXfp+IB9MpdCLiMXJdFgt2JfDN\nqsOsOnCKAD8f+jSO4N42VWlYOeTP35i42977nrgD2jwMN74EfkXwiFwpcnQKnYh4DF8fQ4/6FehR\nvwJ7E1L5ZtVhpm48zuQNsTSLKs29bav+/wN2LAs2fguzn4WAEnDnZKjVzdl/CREPpSt5ESl0Z9Kz\nmbwhlu9+O8zhU2mElgzg1haR3NE4hMor/gY7pkH1TtD/Uwiu4HBakcKl2/Ui4hVcLotl+5L4fvVR\nkves4AP/j6hokjnc+HGq9fk7vr6+TkcUKXS6XS8iXsHHx9CpViidEr/HOvwaZwPCGZHzGkvWVKHy\n3qXc0SqKW1tEElqymNNRRTySruRFxDmp8TBtNBxcAvX7Q+/3yfIvxbyd8Xy/+girDybj72vo1aAi\nd7WuQsuqbn4Snkg+0JW8iHi+fQvsgs86D7eMhWb3gDEEAL0bRdC7UQT7E1P5fvVRpmyMZeaWOGqX\nL8ldravQv2klggMLYKKeiJfRlbyIFK6cLFj4Cvz2EYTXt/e+h9e57FvSsnL4eUsc368+yrbjZwgK\n8KVvk0rc1TqK+hGX2YYn4oG08E5EPNOpAzBlOMRtgpYjoPtr4F/8qn7ElmOn+X71EWZusSfqNY4s\nzZCWkdzSOIISxXRzUjyfSl5EPM/W/8Avj4OPH/T9COrecl0/7kxaNlM2xvLjuqPsTThHiQBf+jSJ\nYMgNUTSsFKJn9+KxVPIi4jkyz8Gsp2HLvyGqDQz4DEpH5tuPtyyLjUdPM3HtUX7ZGkdGtot6FUsx\n5IZI+jatRCk9uxcPo5IXEc9wYos9mjb5IMQ8DTHPgG/B3VI/m5HNjM1xTFxzlJ0nzhLo70PvRhEM\nuSGSZlFamS+eQSUvIu7NsmDNeJj/IgSFwsDPoGr7Qvx4i23HzzBx7TFmbj7O+axcapcvye0toxjQ\nrBKlgzQDX9yXSl5E3Nf5UzDjQdg7B2r3gr4fQ4lyzsXJtFfmT1x3jC3HThPg50OvBhUYckMUraqV\n1dW9uB2VvIi4p0PLYepISDtlr5y/YRS4UYnujDvLj+uOMm3TcVIzcqgeWoLbWkYysHllTdUTt6GS\nFxH3kpsDS/8Fy96GcjXtve8VGzmd6k+lZ+Uya9sJJq49yvojKfj5GG6sG87g5pF0ig7Dz9fH6YhS\nhKnkRcR9nD4GU0bAsdXQ5E7o9RYUK+l0qjzbl5DKpHXHmL75OCfPZREWXIwBTSsxuEVlaoYHOx1P\niiCVvIi4h50zYebD4HJB7/eg0WCnE12z7FwXi3cn8tOGWBbtTiTXZdE0qjSDm0fSu3FFbcWTQqOS\nFxFnZafD3Odg/RcQ0dS+PV+2utOp8k1SaibTNx3nP+uPsS/xHIH+PvRqUJHBzSvTuno5fHzcZ52B\neB+VvIg4J3GXvfc9cSe0fQS6vAh+3rklzbIstsSe4af1x5i5JY7UjBwqlynOoOaVGdisMpFlg5yO\nKF5IJS8ihc+yYOM3MPuv9jP3fuOhVlenUxWajOxc5u6I56f1saw8cBLLgrY1ynFri0h61K9A8QBf\npyOKl1DJi0jhSj8NPz8KO6dD9U7QfwIEl3c6lWNiU9KYsuE4kzce41hyOsHF/OjdOILBLSrTNLK0\n9t7LdVHJi0jhObYWJg+H1Djo8jy0fRR8tMUMwOWyWHMomZ/WH2PW9hNkZLuoFlqC/k0r0b9pJd3O\nl2uikheRgudywcr3YNHrEFIJBn0FlfPl9x2vlJqRzext8UzdFMvqg8kA3FC1LAOaVeKmRlqdL3mn\nkheRgpUaD1NHwaGlUH8A3PI+BIY4ncpjxKakMWNzHFM2xnIw6TwBfj50q1eegc0q0aFWGP4atiOX\noZIXkYKzbz5Mux+yzsNNb0HTu91qNK0nsSyLrbFnmLoxlplb4khJy6ZciQD6NIlgQNPKNKhUSs/v\n5X+o5EUk/+VkwcJX4LePoHwDe+97WLTTqbxGVo6LpXuTmLYplgU7E8nKdVErvCT9m1WiX5NKRJQu\n7nREcRMqeRHJX6cO2HvfT2yGliPtw2X8A51O5bXOpGXz67YTTN0Yy/ojKRgDbaqXY0CzyvRsUIGS\nxfycjigOUsmLSP7ZMgl+fQJ8/OxjYev2djpRkXLk1HmmbTrOtE3HOXIqjUB/H3rUr0DfJhF6fl9E\nqeRF5PplnoNZT8GWiRDVFgZ+BiGVnU5VZFmWxcajKUzdeJxft53gdFo2ZYL86dWwIn0bR9CyalmN\n0y0iVPIicn3iNtu351MOQcwzEPM0+OoWsbvIynGxfF8SM7fEMW9HAunZuVQMCaRP4wj6NImgXkUt\n2PNmKnkRuTaWBWvGw/wXISjUvnqv2t7pVHIZaVk5zN+ZwMzNcSzdm0SOy6JmeEm78BtHUDW0hNMR\nJZ+p5EXk6p0/CdMfhH1zoXYv6PcJBJV1OpVchZTzWczafoIZm+NYe8geuNM4sjR9G0fQu1FFwktp\nsaQ3UMmLyNU5tAymjIT0ZHvl/A2jtPfdw8WdTueXrXHM2BzHjriz+BhoWyOUPk0i6FG/AiHFNWHP\nU6nkRSRvcnNg6Zuw7B0oV9Pe+16xkdOpJJ/tT0xl5uY4ZmyJ48ipNAJ8fehcJ4xbGkfQpU44QQFa\nb+FJVPIicmWnj8KUEXBsDTS5y55eF6Dnt97Msiy2xJ5h5uY4ft4aR1JqJsX9felSN5xbGlWkU3Q4\ngf46EtfdqeRF5PJ2zoSZD9uHzNzyPjQc5HQiKWS5Lou1h5L5dVscs7fFc+p8FiUCfLmxbnl6N6pI\nTO0wFb6bUsmLyKVlp8Pcv8P6LyGiGQz6AspWdzqVOCwn18WaQ8n8svUEc7afICUtm5LF/OhWzy78\n9rVCKeanwncXKnkR+V+Ju+y974k7oe1foMsL4BfgdCpxM9m5Ln47cIpftsYxd0cCZ9KzCQ70o3u9\nCvRuXJF2NUIJ8NOUPSep5EXk/1gWbPga5vwNipWE/uOhZlenU4kHyMpxsfLASX7ZcoJ5O+NJzcgh\npLg/PetX4OZGFWlTo5zG6jpAJS8itvTT8POjsHM6VO8M/T+F4PJOpxIPlJmTy4p9J/ll6wnm70zg\nXGYOZYL86dmgAjc1rEjr6ir8wqKSFxE4thYmD4fUOPvWfNu/gI9+E5brl5Gdy7K9Sfyy9QQLdiWQ\nlpVLSHF/utYtT68GFWhfK1SL9gpQfpa8Nk+KeBpXLqx4Dxb/0z5QZthcqJwvvx+IABDo70v3+hXo\nXr/Cfwt/zvZ45u2MZ8rGWEoW86NznXB6NahAp+gw7cN3Y/pfRsSTnD0B00bZE+zqD7C3xwWGOJ1K\nvNjFhZ+V42LVgZMXCj+Bn7fEEejvQ8faYfRqUJEudcMpFahJe+5Et+tFPMXeeTD9fnubXK+3oOld\nGk0rjsnJdbH2cDJztsczZ3s8iamZBPj60K5mOXo1qEi3euUpU0K7O66FnsmLFCU5mbDgFVj9MZRv\nYI+mDYt2OpXIf7lcFpuOpTBnezyzt8cTm5KOr4+hVbWy9GpQgR71K+jwnKugkhcpKk4dgMlD4cQW\n+1CZbq+Cv36zFPdlWRY74s4ye/sJZm+P52DSeYyB5lFl6FG/At3qldfxuFegkhcpCrZMgl+fAF9/\n6Psx1LnZ6UQiV8WyLPYlnmP2tnjm7Ihn14mzAESXD6ZbvfJ0r1+ehpVCMHrs9P9RyYt4s8xUmPU0\nbJkIUW1h4Gf2KnoRD3csOY35OxOYtzOetYeScVlQMSSQrnXtwm9VrZym7aGSF/FecZvt0bQph6Dj\ns9DhKfDVJhjxPsnns1i0O5H5O+NZujeJjGwXwYF+dI4Op3v98nSsHUZwEV2pr5IX8TaWBavHwfwX\noUQYDPwcqrZzOpVIoUjPymXF/pPM3xnPgl2JJJ/PIsDXhzY1ytG9fnm61S1fpBbuqeRFvMn5kzD9\nAdg3D6Jvsp+/B5V1OpWII3JdFhuOpDB/ZzxzdyRwNDkNgCaRpelevzzd65WnRlhJr36Or5IX8RaH\nlsGUkZCeAt1fgxtGau+7yAWWZbE34RzzdsQzf1cCW2PPAFC1XBBd6pSna91wWlYr63Uz9Qu95I0x\nkcB7QDfAAAuAxyzLOnpVH2bMX4E3gJWWZbXPy3tU8uKVcnNgyRuw/F0oVxMGfwUVGjqdSsStxZ1O\nZ+GuBBbsSuS3A6fIynURXMyPmOgwutYNp1PtcK8YwFOoJW+MCQK2AJnA84AFvAYEAY0syzqfpw8y\npjqwFTgP7FPJS5F1+ihMGQHH1thT63q9BQHaNyxyNc5n5rBi/0kW7Upk4e5ETp7LxMdA8ypl/nuV\nXzPcM2/rF3bJPwqMAaIty9p/4bVqwD7gGcuyxuTpg4yZCxwGogE/lbwUSTtnwMxHwOWy5843HOR0\nIhGP53JZbDt+hoW7Eli4O5EdcfZ+/KiyQXSpE07XuuW5oVpZj9meV9glvxAItCyr3R9eXwpgWVbH\nK36IMXcAH2AX/FRU8lLUZKfDnL/Bhq+gUnMY+AWUreZ0KhGvdOJMOot2J7JwVyIr958kM8dFyWJ+\nxNQOpUud8nSODqNcyWJOx/xThX3UbH1gxiVe3wEMvtKbjTFlsJ/nP2NZVrIn3joRuS6Ju+CnoZC0\nC9o9Cp2fBz/Pf24o4q4qhhTnzlZVuLNVFdKzclm5/yQLdyewcFcis7bFYww0jSxN5+hwOtcJp17F\nUvj4eGc35aXkywIpl3g9GSiTh/e/DewFvs5rKGPMKGAUQFRUVF7fJuJeLMu+cp/zNygWDHdNhZo3\nOp1KpEgpHuBL13rl6Vqv/H/n6i/YlcCi3Ym8O38v787fS2jJYnSKDqNzdDjta4USUtx7hvAU6Cgt\nY0wH4B6gmXUVe/Usy5oATAD7dn0BxRMpOOkpMPMvsGsm1OgC/T+FkuFOpxIp0owxNKgUQoNKITzW\ntTYnz2WybG8Si/ckMX9nApM3xOLrY2geVYaOF0q/bsVgj1y897u8lHwKl75i/7Mr/It9CnwBxBpj\nSl/0mb4X/jndsqzMvIYV8QhH18CU4ZB6Arr9A9o8Aj6eseBHpCgJLVmMAc0qM6BZZXJyXWyJPc2S\nPUks3pPI23P38PbcPZQvVYxOtcPpFB1Gu1qhlPKwUbt5WXi3CAj440I5Y8ySC+//04V3xpgrXYU/\nblnW+5f7Bi28E4/hyoUVY2DxG1A6EgZ+CZWbO51KRK5BYmoGS/cksWRPEsv2JZGakYOfj6F5lTJ0\nrhNO5+hwapcvmC16hb26/jHgHaC2ZVkHL7xWFXsL3V8ty3r3Mu/tdImX3wd8gUeA/ZZlxV7u81Xy\n4hHOnoCpI+HwcmgwCHqPgcAQp1OJSD7IyXWx8ehpluxJZPGepP8emVsxJJBO0WHE1Aqjbc38e5Zf\n2CVfAnsYTjr/NwznVSAYexjOuQvfVwU4APzDsqx/XObnLUFb6MSb7J1rz57PToeb3oYmd2o0rYgX\niz+TwdK9iSzencTK/SdJzczBx9jz9WNqh9GhVhiNK4fgd43jdgt1C51lWeeNMV2wt8F9hz3WdiH2\nWNtzF+fCvkLXw0cpGnIyYcErsPpjKN8ABn0FYbWdTiUiBaxCSCC3tYzitpZRZOe62HLsNMv2JrF0\n30k+WLiP9xfso1SgH+1rhdKhVhgxtcOoVLq4I1l1QI3ItTh1ACYPhRNb4IbR9gI7/6JzFKaIXNrp\ntCxW7D/Jsr1JLNt7kvizGQDUCCtBTG371n6r6mUJCvjza2ydQifipC0/wq9Pgq+/fSxsnZudTiQi\nbsiyLPYnnmPp3iSW7TvJmoOnyMxxEeDrQ8tqZeyr/Fph/7NNTyUv4oTMVLvct06CKu1gwGcQUsnp\nVCLiITKyc1l3OPm/V/l7ElIBCAsuRp/GEbzQux5Q+GNtRSRuE0weBimHodPfIOZp8PF1OpWIeJBA\nf1861LIX5j13MySczbALf99Jcl0Fc8Gtkhe5HMuC1Z/A/JfsiXX3/gJV2135fSIiV1C+VCCDW0Qy\nuEVkgX2GSl7kz5w/aW+N2zcPom+Gvh9BUFmnU4mI5JlKXuRSDi6FqaPsGfQ3vQMtR2jvu4h4HJW8\nyMVys2HJG7B8DITWgrsmQ4WGTqcSEbkmKnmR36UcgSkjIHYtNL0bev0LAko4nUpE5Jqp5EUAdky3\nj4bFgkFfQoOBTicSEbluKnkp2rLSYO7fYMPXUKk5DPwCylZzOpWISL5QyUvRlbDT3vuetAvaPQZd\nnren2ImIeAmVvBQ9lgXrv4S5f4dipeDuaVCji9OpRETynUpeipb0FPvZ+66ZdrH3/9QeciMi4oVU\n8lJ0HF1tr55PPWGfGtfmEfDRycgi4r1U8uL9XLmwYgwsfgNKR8KweVC5udOpREQKnEpevNvZEzB1\nJBxeDg0GQe/3ILCU06lERAqFSl6815459uz5nAzo+wk0uUOjaUWkSFHJi/fJybRPjVszDso3tIfb\nhNV2OpWISKFTyYt3ObkfJg+F+K1ww2h7gZ1/oNOpREQcoZIX77F5Ivz6JPgFwO0Toc5NTicSEXGU\nSl48X2aqXe5bJ0GVdjDgMwip5HQqERHHqeTFs8VtskfTphyGTn+HmKfAx9fpVCIibkElL57J5YLV\nn8CCl+2Jdff9ClXaOp1KRMStqOTF85xLsrfG7Z8PdXpDnw8hqKzTqURE3I5KXjzLwSUwdRSkn4ab\n3oGWI7T3XUTkT6jkxTPkZsPif8KK9yC0Ntw1FSo0cDqViIhbU8mL+0s5AlOGQ+w6aHYP9HwTAko4\nnUpExO2p5MW97ZgGMx8FLHtyXYOBTicSEfEYKnlxT1lpMOevsPEbqNQCBn0BZao6nUpExKOo5MX9\nJOy0R9Mm7YZ2j0GX58HX3+lUIiIeRyUv7sOyYP2XMPfvUKwU3D0NanRxOpWIiMdSyYt7SE+BmY/A\nrp+hxo3Qf7w95EZERK6ZSl6cd+Q3mDICzsVD99eg9UPg4+N0KhERj6eSF+e4cmH5GFjyTygdBcPn\nQaXmTqcSEfEaKnlxxtk4e3Ld4eXQcDDcPAYCSzmdSkTEq6jkpfDtmWPPns/JgL6fQJM7NJpWRKQA\nqOSl8ORkwvyXYM04qNAQBn0FobWcTiUi4rVU8lI4Tu63977Hb4VW90PXV8A/0OlUIiJeTSUvBcuy\nYMtE+PUp8CsGQ36E6F5OpxIRKRJU8lJwMlPhlydg23+gSnsY+BmUinA6lYhIkaGSl4JxfCNMHgan\nj0Dn56DDk+Dj63QqEZEiRSUv+cvlgtUfw4JXoGR5uG8WVGnjdCoRkSJJJS/551wSTL8f9i+AOr2h\nz4cQVNbpVCIiRZZKXvLHgcUwbTSkn4ab34UWw7X3XUTEYSp5uT652bD4dVjxPoTWhrumQoUGTqcS\nERFU8nI9Uo7AlOEQuw6a3Qs934CAEk6nEhGRC1Tycm12TIOZjwKWPbmuwQCnE4mIyB+o5OXqZKXB\nnL/Cxm+gcksY+DmUqep0KhERuQSVvORdwg5773vSHmj/BHT+O/j6O51KRET+hEpersyyYP0XMPc5\nCAyBu6dBjc5OpxIRkStQycvlpSXDzEdg9y9Qsyv0Gw8lw5xOJSIieaCSlz935DeYMgLOJUD316D1\nQ+Dj43QqERHJozz9jm2MiTTGTDbGnDHGnDXGTDXGROXhfS2NMV8YY/YZY9KMMUeNMT8YY6pdf3Qp\nMK5cWPoWfH2T/cx9+Dxo+4gKXkTEw1zxSt4YEwQsAjKBewELeA1YbIxpZFnW+cu8/TagPjAW2AZE\nAC8A640xTSzLOnad+SW/nY2DqaPg8HJoeKs9vS6wlNOpRETkGuTldv1IoDoQbVnWfgBjzFZgHzAa\nGHOZ975lWdZTF79gjFkJHLrwc1+8ltBSQPbMhukPQk4m9BsHjYdoNK2IiAfLy/3XPsDq3wsewLKs\nQ8BKoO/l3mhZVuIlXjsCJAGVri6qFJicTJj9LEy8HUIqw+il0OQOFbyIiIfLy5V8fWDGJV7fAQy+\n2g80xtQFwoFdV/teKQAn98HkoRC/DVo9AN1eAb9iTqcSEZF8kJeSLwukXOL1ZKDM1XyYMcYPGI99\nJf/F1bxX8pllweZ/w6yn7VIfMgmiezqdSkRE8lFhb6H7CGgL3GxZ1qX+4ACAMWYUMAogKuqKi/jl\namWchV+fgG0/QdUOMGAClIpwOpWIiOSzvJR8Cpe+Yv+zK/xLMsa8iV3c91qWNe9y32tZ1gRgAkCL\nFi2svH6G5MHxDTB5OJw+Ap2fhw5PgI+v06lERKQA5KXkd2A/l/+jesDOvHyIMeY54FngEcuyvst7\nPMk3Lhf89hEsfAVKVoD7ZkGVNk6nEhGRApSX1fUzgdbGmOq/v2CMqQq0u/Brl2WM+Qv2vvrnLMv6\n6NpiynU5lwT/HgzzX4DaPeH+5Sp4EZEiIC8l/xlwGJhhjOlrjOmDvdr+GPDp799kjKlijMkxxrx4\n0Wu3A+8Dc4BFxpjWF/1VLz//ReRPHFgM49vBoeVw8xi47XsIKut0KhERKQRXvF1vWdZ5Y0wX4D3g\nO8AAC4HHLMs6d9G3GsCX//8PDj0vvN7zwl8XWwp0uubkcnm52bD4dVjxPoRF2yfHlb/UUxcREfFW\neVpdb1nWUWDgFb7nMHahX/zafcB91xZNrlnKYftgmdh10Pw+6PEGBAQ5nUpERAqZTqHzNtunws+P\nAgYGfw31+zudSEREHKKS9xZZaTDnWdj4LVRuCQO/gDJVnE4lIiIOUsl7g/jtMHkYnNwL7Z+Azn+3\nj4gVEZEiTSXvySwL1n0Oc5+D4qXhnulQvZPTqURExE2o5D1VWjLMfAR2/wI1u9lHw5YMczqViIi4\nEZW8JzqyCqaMhHMJ0P11aP0g+ORl5IGIiBQlKnlP4sqFZe/A0jehdBUYPg8qNXM6lYiIuCmVvKc4\ncxymjoIjK6DhrXDzuxBYyulUIiLixlTynmD3LJjxIORkQb/x0GSI04lERMQDqOTdWXYGzH8R1n4K\nFRrBoK8gtKbTqURExEOo5N3VyX0weSjEb7MX1nV9GfyKOZ1KREQ8iEre3VgWbP4BZj0NfoEwZBJE\n//FsHxERkStTybuTjLPw6xOw7Seo2gEGTIBSEU6nEhERD6WSdxfHN9ijaU8fgy7P2+NpfXydTiUi\nIh5MJe80lwt++wgWvsL/a+/Og6QozziOfx9OFyTFJqjBKCwE1F0q4EGEREE8wqaCIt6EqqAgIh7x\nPktBQwxRMUrQKgHBo9SUB4frRQmCaMQDsdSAiXIIKJSIiBxyL/vmj7dX12Vhh+2ZfWd6fp+qt4bt\nnbd5+qHpZ7qn335p0RoGvQxtuoeOSkREEkBFPqTv1sC0YbB0FhSfBn3vh4LC0FGJiEhCqMiHsnQ2\nTL0Ytm+EPvdC18FgFjoqERFJEBX5+rZrJ8y+A+aOgQOOgIFlcFBJ6KhERCSBVOTr07fLYfKFsGo+\nHHMBlP4dmjQLHZWIiCSUinx9WTgFXrgKMDjnUeh0RuiIREQk4VTkM23HZph+I3zwOBxyLJw1EQrb\nho5KRETygIp8Jq1e6B9Nu3Yx9LgWet0MDRuHjkpERPKEinwmOAfvTYRXboGCljDwOWjfK3RUIiKS\nZ1Tk023LOnj+z/DJi9CxN/R7EJq3Ch2ViIjkIRX5dFrxFkwZ4h9yUzoKul0CDRqEjkpERPKUinw6\nVOyCN0bD63dBYREMmQkHHxU6KhERyXMq8nFtWAVTL4IVc6Fzf+hzDzRtEToqERERFflYPnkZyi6F\n8h3Qbxwc+cfQEYmIiHxPRb4udm6DmSNg3nj4eWc4+xFo1SF0VCIiIj+iIr+vvl7k533/agF0vxRO\nuR0aNQ0dlYiIyG5U5FPlHHz4JLx8PTQugAHPwGGloaMSERHZIxX5VGzbCC9eDQsnQ1EPOPMh+Enr\n0FGJiIjslYp8bVa+D1MGw/ov4KRb4fhroEHD0FGJiIjUSkV+Tyoq4O37YdZIaNEaBk2HNt1CRyUi\nIpIyFfmafLcGpg2DpbOguC/0HQsFhaGjEhER2Scq8tUtmeUL/PaNcOp9cMwgMAsdlYiIyD5Tka+0\nayfM/ivM/SccUAwDy+CgktBRiYiI1JmKPMC6ZTDlQlj1vj9zLx0FTZqFjkpERCQWFfkFk/3wOAzO\neQw69QsdkYiISFrkb5HfsRmm3wAfPAGHdoOzJkLLNqGjEhERSZv8LPKrF/hH065dDD2ug143Q8P8\nTIWIiCRXflU252DeQzDjVj8kbmAZtD8hdFQiIiIZkT9Ffss6KLscPn0JOvaGfg9C81ahoxIREcmY\n/Cjyy+fC1Iv8Q25KR/nZ4zT2XUREEi7ZRb5iF7wxGl6/CwqLYMhMOPio0FGJiIjUi+QW+Q0rYepQ\nWKfrzIoAAAgMSURBVDEXOveHPvdA0xahoxIREak3ySzyn7wEZZdB+Q44Yzx06R86IhERkXqXrCK/\ncxvMHA7zJkDrLnD2I/CzX4aOSkREJIjkFPmvF/mx718tgO6XwSm3QaOmoaMSEREJJveLvHP+qXXT\nb4DGBTDgWTisd+ioREREgsvtIr9tg3/u/MIp0K4nnDEBftI6dFQiIiJZIXeL/Mr3YfIgfxf9ScPh\n+KuhQcPQUYmIiGSNBqm8ycwONbPJZrbBzDaa2VQzS2k2FzPbz8xGm9mXZrbVzN42s551jriiAt4c\nAw/39pfqB02HntepwIuIiFRT65m8mTUDZgPbgfMBB9wBvGZmnZ1zm2tZxSSgD3A98BlwGfCKmf3G\nOffhPkW76St4bhgsnQ0lp8NpY6Gg5T6tQkREJF+kcrn+IqA9cLhzbgmAmf0HWAxcDNy7p45m1gUY\nAAx2zj0SLXsd+BgYCfRNOdIls2DaxbB9E5w6Bo65QI+mFRER2YtULtf3Bd6pLPAAzrllwFzg9BT6\n7gSertK3HHgKKDWz2se4OQczhsMTZ0KzVjB0DnQdpAIvIiJSi1SKfCdgYQ3LPwZKUui7zDm3pYa+\nTYAOtf7t3yyGt8ZC18Ew9DU4sDiFkEVERCSVy/U/Bb6tYfk6oDBG38rf7135Njj3Wf8dvIiIiKQs\nK4fQmdlQYGj043br1K+mKwmSPq2AtaGDSDjlOPOU4/qhPGfe4elaUSpF/ltqPmPf01l69b5t99AX\nfjij/xHn3ARgAoCZzXfOdU0hTqkj5TjzlOPMU47rh/KceWY2P13rSuU7+Y/x361XVwL8N4W+7aJh\neNX77gCW7N5FRERE0iGVIv880N3M2lcuMLMi4Ljod3vzAtAYOKdK30bAecAM59z2fYxXREREUpRK\nkX8IWA6UmdnpZtYXKAO+AMZXvsnM2ppZuZmNqFzmnPsAP3xujJkNMbOT8cPn2gG3pRjjhBTfJ3Wn\nHGeecpx5ynH9UJ4zL205Nudc7W/yj7C9D/gdYMAs4Crn3PIq7ykClgF/cc7dXmV5AfA3/ENxWgIf\nATc65+akZxNERESkJikVeREREck9KU1Qk25ZNeFNQtU1x2b2azObZGaLzWyLmX1uZk+aWbv6iDuX\nxNmPq63nJjNzZvZmJuLMdXHzbGbFZvasma2NjhmfmtmVmYw518Q8Jrcxs8eiY8VWM1tkZneYWfNM\nx51LzOwQM7s/qllbov/zRSn2rXPdq/ciX2XCmyPwE978CeiIn/AmlZ1iEv55+iOAU4Ev8RPeHJmZ\niHNPzByfhx9NMRb4A3ATcDQw38wOzVjQOSYN+3HletoDtwJrMhFnroubZzPrCrwLNAWG4PfpfwCa\ntjISJ8fR718FegLD8fmdCFwLPJzBsHNRB+Bc/NDyf+9j37rXPedcvTbgSmAX0KHKsnZAOXBNLX27\n4GfBG1RlWSPgU+D5+t6WbG0xc3xgDcvaAhXAyNDbli0tTo6rrecV/A2sc4A3Q29XtrWY+3ID/DDf\naaG3I5tbzBz3jo7JpdWW3xn1bxZ6+7KlAQ2q/HlIlLeiFPrFqnshLteHnfAmP9Q5x8653c4onXMr\ngK+BX6Q5zlwWZz8GwMwG4K+S3JyRCJMhTp57AcXsZaZMAeLluEn0ur7a8vX4D1maSSzinKuoY9dY\ndS9EkQ874U1+iJPj3ZhZMXAg8L+YcSVJrBybWSF+xMoNzrkan/woQLw8Hx+97mdm75jZTjNbY2Zj\no1E/4sXJ8av4acfvNrMSM9vfzE7CXx0Y55zbnN5Q81KsuheiyIed8CY/xMnxj0QPLxqHP5OfFD+0\nxIib49HAIuDRNMaURHHyfHD0+jQwAz8E+G78pdJ/pSvABKhzjp1z2/Afphrgi84m/BDrF4HL0xtm\n3opV97JyghrJKg8AvwX6OOdqm6tAUmBmPYCBwNEu+oJNMqLyJOYJ51zlQ7rmmFlD4E4zK3bO6epU\nDGa2H/5D1EH4G/Y+B47F3yBWDlwSLjqBMEW+3ie8yUNxcvw9M7sTPxvg+c65GWmKLSni5Hg8/qrI\nSjNrGS1rBDSMft7q9MjnSnHy/E30OrPa8hn4G8OORF9BQbwcX4i/96Fjle/03zCzDcAEMxvnnPso\nbZHmp1h1L8Tlek14k3lxcgyAmd0C3Ahc4Zx7PI2xJUWcHBcDw/D/eSvbcUD36M86+/lB3OOF1C5O\njn8FrK96015kXvRaHDM2iVn3QhR5TXiTeXFyjJldAdwB3OKceyBDMea6ODk+sYb2Ef7mpxOByekP\nN2fFyfN0YDtQWm3576PX99ITYs6Lk+PVQEszq37zV7fodVWaYsxn8epegLGCzfGfPBbgh2f0xR/g\nPgP2r/K+tvjvdEZU6/8U/mxnCHAy/oC4Df/9ZvCxkNnQ4uQY6I8fEz8df2ZZtZWE3rZsaXH34xrW\nNweNk097nvETYZUDo4BT8A932go8GnrbsqXFPF4UARvxN5Gej/+Qen20bD5VxoarOYCzo/Ygfuz7\nJdHPJ+wpx9HyOte9UBvaBpgS7QibgOeo9lCAaOdxwO3Vlhfgx72ujjbyXaBX6H+8bGt1zTH+bm+3\nhzYn9HZlU4uzH9ewLhX5DOQZP077mqiI7QBWACOBxqG3K5tazByXAM/gZybdGhX8e4DC0NuVba22\nY2sm6p4mqBEREUmoIBPUiIiISOapyIuIiCSUiryIiEhCqciLiIgklIq8iIhIQqnIi4iIJJSKvIiI\nSEKpyIuIiCSUiryIiEhC/R88XbbHQ4XUggAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1112039b0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from draw import graph_two_functions\n",
    "import numpy as np \n",
    "f1 = lambda x: np.exp(-x)\n",
    "f2 = lambda x: x\n",
    "graph_two_functions(f1, f2, x1=0, y1=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Definition\n",
    "\n",
    "Suppose that $g(x)$ has a unique fixed point on the interval $[a,b]$.  Given $p_0 \\in [a,b]$ (an initial guess), define $p_n = g(p_{n-1})$ for $n \\geq 1$.\n",
    "\n",
    "<br>\n",
    "\n",
    "If the sequence $p_n \\to p$ as $n \\to \\infty$, and if $g$ is continuous\n",
    "\n",
    "$$ p = \\lim_{n\\to\\infty} p_n = \\lim_{n \\to \\infty} g(p_{n-1}) = g \\left( \\lim_{n \\to \\infty} p_{n-1} \\right) = g(p).$$\n",
    "\n",
    "The technique of constructing the sequence $\\{g_n\\}_{n = 1}^\\infty$ is called __fixed-point iteration__.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Example\n",
    "\n",
    "Use fixed-point iteration to find the fixed point(s) of $g(x) = 1 + \\sqrt{x}$, $x \\geq 0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p_n\t\t g(p_n)\n",
      "0.100000\t 1.316228\t\n",
      "1.316228\t 2.147270\t\n",
      "2.147270\t 2.465357\t\n",
      "2.465357\t 2.570145\t\n",
      "2.570145\t 2.603167\t\n",
      "2.603167\t 2.613433\t\n",
      "2.613433\t 2.616612\t\n",
      "2.616612\t 2.617594\t\n",
      "2.617594\t 2.617898\t\n",
      "2.617898\t 2.617992\t\n",
      "difference between approximate and exact fixed points =  4.19776218132e-05\n"
     ]
    }
   ],
   "source": [
    "from fixed_point_iteration_helper import print_heading, print_row\n",
    "import numpy as np \n",
    "\n",
    "def fixed_point_iteration(func, p0, n_max):\n",
    "    print_heading()     \n",
    "    p = p0\n",
    "    for n in range(n_max):\n",
    "        print_row(p, func(p))\n",
    "        p = func(p)\n",
    "    return p\n",
    "\n",
    "g = lambda x: 1 + np.sqrt(x)\n",
    "p_approx = fixed_point_iteration(g, p0=0.1, n_max=10)\n",
    "p_exact = ((1 + np.sqrt(5))/2)**2\n",
    "print('difference between approximate and exact fixed points = ', abs(p_approx - p_exact))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Convergence of fixed point iteration\n",
    "\n",
    "When does fixed point iteration converge, and at what rate?\n",
    "\n",
    "#### Theorem 2\n",
    "\n",
    "Let $g \\in C[a,b]$ be such that $g(x) \\in [a,b]$ for all $x \\in [a,b]$.  Suppose, in addition, that $g'$ exists on $(a,b)$ and there exists $0 \\leq K < 1 $ such that\n",
    "\n",
    "$$ |g'(x)| \\leq K, \\quad \\text{ for all } ~~ x \\in [a,b]. $$\n",
    "\n",
    "Then for any number $p_0 \\in [a,b]$, the sequence $p_n \\to p$ at rate $K^n$ where $p$ is the unique fixed point in $[a,b]$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Proof\n",
    "\n",
    "From Theorem 1, we know that the unique fixed point $p$ exists.  Then\n",
    "\n",
    "$$ |p_n - p| = |g(p_{n-1}) - g(p)| \\leq K|p_{n-1} - p| \\leq K^2 |p_{n-2} - p| \\leq K^n |p_0 - p|. $$\n",
    "\n",
    "It also follows that $|p_0 - p| \\leq b-a$, the length of the interval that contains them both.  And so\n",
    "\n",
    "$$|p_n - p| \\leq K^n (b-a), \\quad \\text{or} \\quad p_n = p + O(K^n),$$\n",
    "\n",
    "as $n \\to \\infty$.  Convergence follows because $K < 1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Example \n",
    "Use fixed-point iteration to find the root of $f(x) = x^4 + 2x^2 -x - 3$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p_n\t\t g(p_n)\n",
      "1.000000\t 0.000000\t\n",
      "0.000000\t -3.000000\t\n",
      "-3.000000\t 96.000000\t\n",
      "96.000000\t 84953085.000000\t\n",
      "84953085.000000\t 52085473681492913685918814044160.000000\t\n"
     ]
    }
   ],
   "source": [
    "# Method 1: find fixed point of g_1(x)\n",
    "g1 = lambda x: x**4 + 2*x**2 - 3\n",
    "root = fixed_point_iteration(g1, p0=1, n_max=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p_n\t\t g(p_n)\n",
      "1.000000\t 1.189207\t\n",
      "1.189207\t 1.080058\t\n",
      "1.080058\t 1.149671\t\n",
      "1.149671\t 1.107821\t\n",
      "1.107821\t 1.133932\t\n",
      "1.133932\t 1.118003\t\n",
      "1.118003\t 1.127857\t\n",
      "1.127857\t 1.121813\t\n",
      "1.121813\t 1.125540\t\n",
      "1.125540\t 1.123249\t\n"
     ]
    }
   ],
   "source": [
    "# Method 2: find fixed point of g_2(x)\n",
    "g2 = lambda x: (3 + x - 2*x**2)**0.25\n",
    "root = fixed_point_iteration(g2, p0=1, n_max=10)"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}