{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Foundations of Computational Economics #25\n", "\n", "by Fedor Iskhakov, ANU\n", "\n", "" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "## Newton-Raphson method with bounds\n", "\n", "" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "\n", "\n", "[https://youtu.be/1fD-BSiO1f4](https://youtu.be/1fD-BSiO1f4)\n", "\n", "Description: Robust implementation of Newton method for problems with strict bounds." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "- example of how elements of different methods can be combined\n", " in a poly-algorithm \n", "- deeper merge than a multi-stage poly-algorithm we mentioned before, such as\n", " 1. Robust slow method for initial search of parameter space\n", " 2. Fast accurate method invoked from a very good starting values " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Example\n", "\n", "Solve the following equation\n", "\n", "$$\n", "f(x) = a \\log(x) + b \\log(1-x) + c = 0, \\; ab<0\n", "$$\n", "\n", "This equation arises in the models of discrete choice in IO, for example when computing a mixed\n", "strategy equilibrium in a two players game with binary actions." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Theoretical properties and numerical difficulties\n", "\n", "- $ x \\in (0,1) $ are strict limits, any algorithm will break down if stepping outside \n", "- there is exactly one solution for any values of parameters $ a,b,c $ (where $ a $ and $ b $ have opposite signs, without loss of generality assume $ a>0 $ and $ b<0 $)\n", " - $ f(x) $ is continuous in its domain $ x \\in (0,1) $\n", " - $ f'(x) \\ge 0 $ everywhere in the domain, so the function is monotonically increasing, and vise versa\n", " - when $ x \\rightarrow 0 $ from the right $ f(x) \\rightarrow -\\infty $, and when $ x \\rightarrow 1 $ from the left $ f(x) \\rightarrow \\infty $, and vise versa \n", "- but depending on the value of $ c $ the solution may be arbitrary close to 0 or 1! " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "- Newton — can easily overshoot to the outside of domain \n", "- bisections — may take forever to converge \n", "- successive approximations — ? " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Newton-bisection hybrid algorithm\n", "\n", "1. Initialize with interval $ (a,b) $ such that the function\n", " has opposite signs at the ends, and starting value $ x_0 $ \n", "1. Compute the Newton step and check if it remains within the bounds \n", "1. If yes, make Newton step and continue to next iteration \n", "1. If not, make a bisection step and continue to next iteration \n", "1. Iterate until convergence, or until maximum number of iterations is exceeded " ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "hide-output": false, "slideshow": { "slide_type": "slide" } }, "outputs": [ { "data": { "text/plain": [ "{'divide': 'warn', 'over': 'warn', 'under': 'ignore', 'invalid': 'warn'}" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline\n", "plt.rcParams['figure.figsize'] = [12, 8]\n", "np.seterr(all=None, divide='ignore', over=None, under=None, invalid='ignore') # turn off log(0) warning" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "hide-output": false, "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "# code up the function\n", "# make a plot of the function\n", "# code up the solver" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "hide-output": false, "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "#\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "# solution below\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "#" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "hide-output": false, "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "def newton_bounds(fun,grad,x0=None,bounds=(0,1),\n", " tol=1e-6,maxiter=100,callback=None):\n", " '''Polyalgorithm that combines bisections and Newton-Raphson\n", " to solve fun(x)=0 within given lower and upper bounds.\n", " Callback function is invoked at each iteration if given.\n", " '''\n", " a,b = bounds\n", " sa,sb = np.sign(fun(a)),np.sign(fun(b)) # a and b signs, never change\n", " if sa*sb > 0:\n", " raise(ValueError('Function has the same signs at the initial lower and upper limits'))\n", " x0 = (a+b)/2 if x0==None else x0 # midpoint is default x0\n", " for i in range(maxiter):\n", " f0 = fun(x0)\n", " newt = x0 - f0/grad(x0) # Newton step\n", " if not a < newt < b:\n", " a,b = (x0,b) if np.sign(f0)*sa > 0 else (a,x0) # update limits\n", " x1 = (a+b)/2 # bisections step\n", " step_type = 'bisection'\n", " else:\n", " x1 = newt\n", " step_type = 'newton'\n", " err = np.abs(x0-x1) # save error for both steps\n", " if callback:\n", " callback(iter=i,type=step_type,err=err,x0=x0,x1=x1,bounds=(a,b))\n", " if err < tol:\n", " break\n", " x0 = x1\n", " else:\n", " raise(RuntimeError('Failed to converge in %d iterations'%maxiter))\n", " return x1" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "hide-output": false, "scrolled": false, "slideshow": { "slide_type": "slide" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsEAAAHSCAYAAAANGxbcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deZiU1Z328e8BGhARBQFxQAQccAGRGESMGo1iVBLjHs1EYjIqMY7rmLwzyZhMlpmJWRyzGaNGY+Ikxj0uo7gjxh2QHY2KCgIKyKLsNJz3j9M93SBLNV1dT1U938911fV0V5Xdd/JEc3v41TkhxogkSZKUJ62yDiBJkiSVmiVYkiRJuWMJliRJUu5YgiVJkpQ7lmBJkiTljiVYkiRJudMmi1/atWvX2KdPnyx+tSRJknJkwoQJi2KM3TZ9PpMS3KdPH8aPH5/Fr5YkSVKOhBDe3tzzjkNIkiQpdyzBkiRJyh1LsCRJknLHEixJkqTcsQRLkiQpdyzBkiRJyh1LsCRJknLHEixJkqTcsQRLkiQpdyzBkiRJyh1LsCRJknLHEixJkqTcsQRLkiQpdyzBkiRJyh1LsCRJknLHEixJkqTcsQRLkiSpZaxbB0uWQG1t1kk+whIsSZKklvHss9ClC4wbl3WSj7AES5IkqWWsW5euNTXZ5tgMS7AkSZJahiVYkiRJuVNfgtu0yTbHZliCJUmS1DLqPxDnSrAkSZJyw3EISZIk5Y4lWJIkSbljCZYkSVLuWIIlSZKUO5ZgSZIk5Y4lWJIkSbljCZYkSVLuWIIlSZKUO5ZgSZIk5c66ddCqVXqUmfJLJEmSpOqwbl1ZrgKDJViSJEktZd06aNMm6xSbZQmWJElSy6itdSVYkiRJOeM4hCRJknLHEixJkqTcsQRLkiQpdyzBkiRJyh1LsCRJknLHEixJkqTcsQRLkiQpdyzBkiRJyh1LsCRJknLHEixJkqTcsQRLkiQpd9atgzZtsk6xWZZgSZIktYzaWleCJUmSlDOOQ0iSJCl31q6t/BIcQtgjhPBkCGFmCGF6COGSuue/G0KYG0KYVPcY2XJxJUmSVDFWr4Yddsg6xWY1ZVK5Frg8xjgxhLATMCGE8Gjda1fHGH9a/HiSJEmqWKtXQ/v2WafYrIJLcIxxPjC/7usPQwgzgZ4tFUySJEkVbtWqsl0J3q6Z4BBCH+BjwAt1T10YQpgSQrgphNB5C3/N6BDC+BDC+IULF25XWEmSJFWIGMt6JbjJJTiE0BG4C7g0xvgBcC2wFzCEtFJ81eb+uhjj9THGoTHGod26dWtGZEmSJJW9detSEa6GleAQQg2pAP8xxng3QIzxvRjj+hjjBuAGYFjxY0qSJKmirFqVrpW+EhxCCMCNwMwY4383en73Rm87GZhWvHiSJEmqSKtXp2uZluCm7A5xKDAKmBpCmFT33LeAL4QQhgAReAv4alETSpIkqfLUrwSX6ThEU3aH+CsQNvPSg8WLI0mSpKpQ5ivBnhgnSZKk4qsvwWW6EmwJliRJUvFVywfjJEmSpII5DiFJkqTcKfMPxlmCJUmSVHyuBEuSJCl3LMGSJEnKHcchJEmSlDuuBEuSJCl3XAmWJElS7qxcma6WYEmSJOXGihXQrh20aZN1ks2yBEuSJKn4li+Hjh2zTrFFlmBJkiQVnyVYkiRJuWMJliRJUu5YgiVJkpQ7lmBJkiTlzooVsOOOWafYIkuwJEmSis+VYEmSJOWOJViSJEm5YwmWJElSrmzYkGaCLcGSJEnKjZUr09UPxkmSJCk3PvwwXXfaKdscW2EJliRJUnEtW5auO++cbY6tsARLkiSpuCzBkiRJyp36ErzLLtnm2ApLsCRJkopr6dJ0dSVYkiRJueE4hCRJknLHEixJkqTcWbYMQvCwDEmSJOXIsmVpFbhV+VbN8k0mSZKkylRfgsuYJViSJEnFtXRpWW+PBpZgSZIkFduiRdC1a9YptsoSLEmSpOKyBEuSJCl3Fi60BEuSJClHamthyRJLsCRJknJk8eJ07dYt2xzbYAmWJElS8SxalK6uBEuSJCk3LMGSJEnKnYUL09VxCEmSJOWGK8GSJEnKnfoSvOuu2ebYBkuwJEmSimfRIthpJ2jXLuskW2UJliRJUvEsXFj288BgCZYkSVIxVcCRyWAJliRJUjFZgiVJkpQ7CxdWVwkOIewRQngyhDAzhDA9hHBJ3fNdQgiPhhBeq7t2brm4kiRJKlsxwrvvQo8eWSfZpqasBNcCl8cY9wWGA/8UQtgP+Ffg8Rhjf+Dxuu8lSZKUN4sWwdq10LNn1km2qeASHGOcH2OcWPf1h8BMoCdwIvD7urf9Hjip2CElSZJUAebOTddqKsGNhRD6AB8DXgB2izHOh1SUge7FCidJkqQKUs0lOITQEbgLuDTG+EET/rrRIYTxIYTxC+vPlJYkSVL1qNYSHEKoIRXgP8YY7657+r0Qwu51r+8OLNjcXxtjvD7GODTGOLRbBWygLEmSpCaaOxdCqK4PxoUQAnAjMDPG+N+NXroPOLvu67OBe4sXT5IkSRVj7lzYbTeoqck6yTa1acJ7DwVGAVNDCJPqnvsWcCVwewjhHGA2cHpxI0qSJKkizJ1bEaMQ0IQSHGP8KxC28PLRxYkjSZKkijV3LvTtm3WKgnhinCRJkoqjglaCLcGSJElqvlWrYPFiS7AkSZJypH57tF69ss1RIEuwJEmSmm/WrHR1JliSJEm58eab6dqvX7Y5CmQJliRJUvPNmgVt28Lf/V3WSQpiCZYkSVLzzZqVRiFaVUa9rIyUkiRJKm+zZlXMKARYgiVJklQMlmBJkiTlypIlsHSpJViSJEk5Ur89miVYkiRJufHGG+laIXsEgyVYkiRJzfXKKxAC9O+fdZKCWYIlSZLUPK+8AnvuCR06ZJ2kYJZgSZIkNc/MmbDvvlmnaBJLsCRJkrbfhg3w6quwzz5ZJ2kSS7AkSZK23+zZsGqVK8GSJEnKkZkz09USLEmSpNywBEuSJCl3pk+Hbt1g112zTtIklmBJkiRtv8mT4YADsk7RZJZgSZIkbZ/aWpg2zRIsSZKkHHn1VVizBoYMyTpJk1mCJUmStH0mT05XV4IlSZKUG5MnQ9u2FXdQBliCJUmStL0mTYKBA6GmJuskTWYJliRJUtPFCOPHw4EHZp1ku1iCJUmS1HRvvAGLF8PBB2edZLtYgiVJktR0L76YrsOGZZtjO1mCJUmS1HQvvggdOqSZ4ApkCZYkSVLTvfgifPzj0KZN1km2iyVYkiRJTbNmDUycWLGjEGAJliRJUlONH5+K8GGHZZ1ku1mCJUmS1DRPPZWuhx+ebY5msARLkiSpaZ56CgYNgl13zTrJdrMES5IkqXC1tfDMM3DEEVknaRZLsCRJkgo3cSKsWAGf/GTWSZrFEixJkqTCjRuXrpZgSZIk5cZTT8GAAdCjR9ZJmsUSLEmSpMKsXZtK8Kc+lXWSZrMES5IkqTDPPgsffgjHHZd1kmazBEuSJKkwDz0ENTVw9NFZJ2k2S7AkSZIKM2ZMOiVup52yTtJslmBJkiRt29y5MGUKHH981kmKwhIsSZKkbRszJl2rYB4YLMGSJEkqxP/+L/TsmY5LrgKWYEmSJG3dihVpJfikkyCErNMUhSVYkiRJWzdmDKxaBaeemnWSoim4BIcQbgohLAghTGv03HdDCHNDCJPqHiNbJqYkSZIyc9dd0LUrHH541kmKpikrwTcDm5uEvjrGOKTu8WBxYkmSJKksrFkDDzyQRiHatMk6TdEUXIJjjOOAxS2YRZIkSeXmscfSKXGnnJJ1kqIqxkzwhSGEKXXjEp2L8PMkSZJULm69FXbZpSpOiWusuSX4WmAvYAgwH7hqS28MIYwOIYwPIYxfuHBhM3+tJEmSWtzy5XDPPXDGGdC2bdZpiqpZJTjG+F6McX2McQNwAzBsK++9PsY4NMY4tFu3bs35tZIkSSqFu++GlSth1KiskxRds0pwCGH3Rt+eDEzb0nslSZJUYW65Bfr2hU98IuskRVfwR/xCCLcCRwJdQwjvAP8OHBlCGAJE4C3gqy2QUZIkSaU2dy48/jhccUXVHJDRWMElOMb4hc08fWMRs0iSJKlc3HwzxFiVoxDgiXGSJEna1Pr1cP31aUeI/v2zTtMiLMGSJEna2EMPwezZcP75WSdpMZZgSZIkbezaa6FHDzjxxKyTtBhLsCRJkhq8+WZaCT73XKipyTpNi7EES5IkqcENN6TdIEaPzjpJi7IES5IkKVmxIn0g7oQTYI89sk7ToizBkiRJSn73O3j/ffjGN7JO0uIswZIkSYLaWrjqqnQ63KGHZp2mxRV8WIYkSZKq2J13wltvwc9+lnWSknAlWJIkKe9ihB//GPbZJ80D54ArwZIkSXk3Zgy8/DL89rfQKh9rpPn4TylJkqTNixG+/W3o0wdGjco6Tcm4EixJkpRn994LEyaknSHats06Tcm4EixJkpRXGzbAd74DAwbAWWdlnaakXAmWJEnKqzvugKlT4Y9/hDb5qoWuBEuSJOXRmjXwrW/BwIFwxhlZpym5fFV+SZIkJb/8JcyalXaGaN066zQl50qwJElS3ixcCD/4ARx/PBx7bNZpMmEJliRJypvvfhdWrICf/jTrJJmxBEuSJOXJ5Mlw3XXw1a/CfvtlnSYzlmBJkqS8WL8+ld8uXdI4RI75wThJkqS8uO46eOEFuOWWVIRzzJVgSZKkPJg/H775TTj6aPjiF7NOkzlLsCRJUh5cemnaG/jaayGErNNkzhIsSZJU7e65B26/Ha64Avr3zzpNWbAES5IkVbMFC9KH4Q48EP7lX7JOUzYswZIkSdUqRjjvPPjgg/RhuJqarBOVDXeHkCRJqlY33wz33Qf//d+53hN4c1wJliRJqkZvvAGXXAJHHpmu2oglWJIkqdqsXg2f/zy0aZNWg1tZ+TblOIQkSVK1+frXYeLENAqx555ZpylL/muBJElSNbn9drjmGrj8cjjhhKzTlC1LsCRJUrV49VU491wYPhx++MOs05Q1S7AkSVI1WLIkrfy2bw+33eZ2aNvgTLAkSVKlq62FM86At96CJ56A3r2zTlT2LMGSJEmV7utfh0cfhd/+Fg47LOs0FcFxCEmSpEp2443w85/DpZfCOedknaZiWIIlSZIq1ZgxcP758OlPw09+knWaimIJliRJqkQvvQSnnQaDBsEdd6SDMVQwS7AkSVKlee01GDkSuneHhx6CTp2yTlRxLMGSJEmV5N134dhj09cPPww9emSbp0K5bi5JklQpFi6EESNgwQJ48kno3z/rRBXLEixJklQJ3n8/FeA33oAHH4SDDso6UUWzBEuSJJW7JUvgmGPSscj33w+f+lTWiSqeJViSJKmcLV0Kxx0H06fDX/6SyrCazRIsSZJUrhYsSHsAz5gBd94Jxx+fdaKqYQmWJEkqR3PmpFXf2bPTCET9jhAqCkuwJElSuXn9dTj66DQK8cgjcNhhWSeqOpZgSZKkcjJlSlr1XbcubYN24IFZJ6pKBR+WEUK4KYSwIIQwrdFzXUIIj4YQXqu7dm6ZmJIkSTlQv+rbujWMG2cBbkFNOTHuZuC4TZ77V+DxGGN/4PG67yVJktRUN90En/kM9O0Lzz8P++2XdaKqVnAJjjGOAxZv8vSJwO/rvv49cFKRckmSJOVDjPDtb8M558BRR8HTT0OvXlmnqnpNWQnenN1ijPMB6q7dt/TGEMLoEML4EML4hQsXNvPXSpIkVYGVK+Gss+A//iOV4AcegE6dsk6VC80twQWLMV4fYxwaYxzarVu3Uv1aSZKk8vT222n+99Zb4Yc/hBtugJqarFPlRnN3h3gvhLB7jHF+CGF3YEExQkmSJFW1sWPh9NPTDhAPPAAjR2adKHeauxJ8H3B23ddnA/c28+dJkiRVrxjhl7+EESOga1d48UULcEaaskXarcBzwN4hhHdCCOcAVwLHhBBeA46p+16SJEmbWr4cvvQluPjitAvECy/AgAFZp8qtgschYoxf2MJLRxcpiyRJUnWaPBk+//l0Etz3vgdXXAGtSvbRLG2GJ8ZJkiS1lBjhuuvg0kuhSxd4/HE48sisU4kS7g4hSZKUK8uWwZlnwte+lorvpEkW4DJiCZYkSSq2J56A/feHu+6CK6+EBx+E7ls8TkEZsARLkiQVy6pVcNllcPTRsMMO8Mwz8C//4vxvGXImWJIkqRgmTIBRo2DmTLjwQvjRj6BDh6xTaQv81xJJkqTmWLs27fgwfDh88AE8/HDaC9gCXNZcCZYkSdpeL7wA554L06bBF74A11wDnTtnnUoFcCVYkiSpqZYvh0sugUMOgaVL4b774E9/sgBXEFeCJUmSmuKhh+D882HOHLjgAviv/4JOnbJOpSZyJViSJKkQc+bAGWfAyJGw447w9NPwq19ZgCuUJViSJGlr1qyBH/4Q9tknjT1873vw8stw6KFZJ1MzOA4hSZK0JWPGwMUXw2uvwUknwdVXQ58+WadSEbgSLEmStKk334STT4bjj0/fP/QQ3HOPBbiKWIIlSZLqLV4Ml1+eRh8eeSSNQUydCscdl3UyFZnjEJIkSWvWwK9/DT/4Qdry7Ctfge9/H3r2zDqZWogrwZIkKb9ihNtvh/32g3/+Zxg2DCZNghtvtABXOUuwJEnKpyefhE98Im171rFjOu54zBgYPDjrZCoBS7AkScqXZ5+Fo45Kjzlz4KabYOJE+PSns06mErIES5KkfJgwIR10ceihMGMG/Pzn8Prraf63deus06nELMGSJKm6TZ0Kp5wCQ4fCCy/Aj34Eb7yR9v9t3z7rdMqIu0NIkqTq9NJL8J//Cffem442/t734NJLPeZYgCVYkiRVm3HjUvl95BHYZRf4znfgkkugS5esk6mMWIIlSVLlizHt7vCf/wl//St07w5XXglf+5orv9osS7AkSapctbVw113wk5+kD7716gW/+AWccw506JB1OpUxS7AkSao8H3wAv/1t2uFh9mz4+7+HG26AL30J2rbNOp0qgCVYkiRVjjlzUvG94YZUhA8/PK38nnACtHLTKxXOEixJksrfhAlw1VXpiGOA006Dyy+Hgw7KNpcqliVYkiSVpzVr4M474Zpr4LnnYKed0i4PF18Me+6ZdTpVOEuwJEkqL3PmwHXXpZGHBQugf3+4+up0stvOO2edTlXCEixJkrIXIzz5ZFr1vfde2LABPvtZuPBCGDHCeV8VnSVYkiRlZ/Fi+J//gd/8BmbOhF13ha9/Hc4/H/r0yTqdqpglWJIkldaGDfDUU2mLs7vuSrO/Bx0EN98MZ5wB7dtnnVA5YAmWJEmlMX9+Kro33ghvvJGOND7vvHSwxZAhWadTzliCJUlSy6mthTFj0qrvAw/A+vVwxBHw3e/CqafCDjtknVA5ZQmWJEnFN2UK3HIL/PGPaQV4t93SrO8//iMMGJB1OskSLEmSiuTdd+FPf4I//AEmT4Y2beAzn4Gzz047PdTUZJ1Q+j+WYEmStP1Wrkxbmv3hD/DII+lDb8OGwa9+lT7k1rVr1gmlzbIES5Kkplm/HsaOTaMOd94JH34IvXvDN78JZ50F++yTdUJpmyzBkiRp2zZsgGefhT//ORXf995LxxiffjqMGgWf/KQHWqiiWIIlSdLmxQjjx8Ntt6XHO++kPXw/+1k480wYOdLdHVSxLMGSJKlBjDB1aiq9f/4zzJqVPtB23HFw5ZXwuc+lFWCpwlmCJUnKuxjh5Zfh7rvTCW6vvAKtW8NRR8G//RucfDJ07px1SqmoLMGSJOVR/Yzv3Xenx9tvp5neI46Aiy6C006D7t2zTim1GEuwJEl5sW5d2tXh7rvhL39J+/q2bQvHHAPf+U4adXBLM+WEJViSpGq2ahU8+mgac7j/fliyBHbcMX2o7ZRT0rVTp6xTSiVnCZYkqdrMmwcPPJBK7+OPpyLcuXNa6T3llLTy664OyjlLsCRJla7+g233358eEyak5/v0gXPPTeX3iCM8tlhqxBIsSVIlWrUqrfI+8EB6zJ0LIcAhh8APfwgnnAD77Zeek/QRRSnBIYS3gA+B9UBtjHFoMX6uJElq5O234eGHU+l97LFUhDt2hGOPTaV35Ejo1i3rlFJFKOZK8KdijIuK+PMkScq31ath3DgYMyY9Zs5Mz++5J5xzTiq+RxwB7dplm1OqQI5DSJJULmKE115rKL1jx6bV3nbtUtkdPTqd3Lb33o45SM1UrBIcgUdCCBG4LsZ4/aZvCCGMBkYD9O7du0i/VpKkCrd8OTzxREPxffPN9PyAAXDeean0HnEEdOiQbU6pyhSrBB8aY5wXQugOPBpCeCXGOK7xG+qK8fUAQ4cOjUX6vZIkVZb169PuDY89lh5//Ws6xGLHHeHoo+Eb30gzvv36ZZ1UqmpFKcExxnl11wUhhHuAYcC4rf9VkiTlQIzwt781lN6xY2Hp0vTa4MFw6aVptffQQ53tlUqo2SU4hLAj0CrG+GHd158Gvt/sZJIkVar589P2ZfXFd+7c9Pyee8Kpp8KIEXDUUdC9e7Y5pRwrxkrwbsA9IQ3otwH+FGMcU4SfK0lSZfjgA3jqqYbSO2NGer5LlzTiMGJEuvbr5wfapDLR7BIcY5wFHFCELJIkVYYPPkizvGPHpseECbBhQzqK+PDD4ctfTsX3gAOgVauMw0raHLdIkyRpW5Ytg6efTqu9Y8fCxImp9NbUwPDh8K1vpZXeQw5xrleqEJZgSZI2tXTpxqX35ZdT6W3bNpXeK65I25YNH+7WZVKFsgRLkvT++/DMM6nwPvVUKr0xplXd4cPh29+GI4+Egw9OIw+SKp4lWJKULzHCW2+lmd76R/0H2dq1SyMN//7vDaW3ffss00pqIZZgSVJ1W78epkzZuPTOm5de23ln+MQn4ItfhMMOg2HDLL1STliCJUnVZcUKePHFhsL73HPw4YfptT32SLO8hx2WHgMHQuvW2eaVlAlLsCSpss2bl4rus8+m0jtxItTWpv14Bw2Cs85qKL29e2edVlKZsARLkirHmjXpQ2vPP5+K7/PPw+zZ6bV27dI4wze+kQrvIYdA587Z5pVUtizBkqTyNWfOxoV3wgRYuza9tsceqehedlnaweFjH3OPXkkFswRLksrD6tVplKG+8D73HMydm15r1w6GDoWLL06Fd/hw6Nkz27ySKpolWJJUejHCm2+mD7DVF96XX4Z169LrffrAJz+Zyu4hh6Tjh9u2zTSypOpiCZYktbx334WXXkqPF1+E8ePTARWQDp846KA01nDIIan49uiRbV5JVc8SLEkqrg8+SCW3cemdMye91qpV2pbspJNS8T3oINh/f6ipyTazpNyxBEuStt+aNTB5ciq69YX31VfTuANAv37pMIphw1LhPfBA2HHHbDNLEpZgSVKh1q+HV17ZuPBOmdIwx7vbbqno/sM/pOvQodC1a7aZJWkLLMGSpI+qrYUZM9KWZBMnpuukSbBqVXp9p51Syb3ssoZV3j32SAdUSFIFsARLUt6tXQvTpzeU3QkT0grv6tXp9R13THvwjh6dxhkOOgj23jvN90pShbIES1KerFkDU6duXHinTm04gGKnnVLRveCCdP34x6F/f2jdOtvcklRklmBJqlarVqUV3frCO3EiTJvWMMO7yy6p6F5ySUPh3WsvV3gl5YIlWJKqwdKlaZeGSZPSY+LENOKwfn16vUuXVHIvv7yh8Pbt6wyvpNyyBEtSJYkRZs9uKLv1j7feanhPjx4wZAiccEJD4e3d28IrSY1YgiWpXK1dCzNnfrTwLl2aXg8BBgyAgw+G889PxfeAAzxtTZIKYAmWpHKw6TjDpElpnKF+fneHHWDwYDjjjFR2hwxJJ6158IQkbRdLsCSVUozpCOH6ovvyyx8dZ+jePW1JduyxDYXXHRokqagswZLUUlasSLsxTJmStiGbMiU9lixJrzceZ/jqVxsKr+MMktTiLMGS1FwbNsCsWQ0lt77wvvFGWvkF6NgxjS+cfvrG4wwdO2abXZJyyhIsSU2xePHGq7pTp6bHypXp9RDS6MKQITBqVJrjHTwY+vRx/11JKiOWYEnanHXr4NVXG8pufeF9552G93TpknZjOPfchrI7cCB06JBdbklSQSzBkvItRnj33Y3L7pQpaWuy+p0Zampg333hyCNT0d1//3TdfXf33pWkCmUJlpQfH34IM2akD6s1Hml4//2G9/TsmQru8cc3lN2994a2bbPLLUkqOkuwpOqzenUaZZg2beNH423IOnSAQYPg5JMbVnf33x923TWz2JKk0rEES6pctbVpB4ZNy+5rr8H69ek9bdrAPvvA8OFpdnfgwFR++/Z1311JyjFLsKTyt2EDzJ790bI7c2Y6WhjSbO5ee6WCe9pp6TpoUNqpwVEGSdImLMGSykeM8N57Hy2706fD8uUN79tjj1Rwjzmmoezuu6+7MkiSCmYJlpSNJUtSud208Db+kFrXrmlO9ytfaSi7++0Hu+ySXW5JUlWwBEtqWStWNOzI0Pgxb17Dezp1SgX31FMbZnYHDYLu3bPLLUmqapZgScWxdu3md2R4882Go4Pbt08ruSNGNBTdQYOgVy/325UklZQlWFLTrF8Ps2Z9tOz+7W9ptwZIOzLsvTcMHQpf/nJD2e3Xzx0ZJEllwRIsafPqd2Son9utv86cmfbhhbR6269fw3679WV3wAB3ZJAklTVLsJR3McL8+RvvxDBtWprjbbwjQ69eqeAeddTGOzLsuGN22SVJ2k6WYClPFi7cuOjWX5cubXjPbrulD6f94z82fEjNHRkkSVXGEixVo6VLP1p0p01LJbhe586p4J55ZroOHJge3bpll1uSpBKxBEuVbPnyNLbQuOhOnw5z5za8p2PHVG4/97mNtx/r0cMdGSRJuWUJlirBqlXwyisfXd19662G99RvP3b00Q1ld+BA6N3bsitJ0iYswVI5WbcubTW26dzu66+n3RoAampgn31g+HA499yGwtu3r9uPSZJUIEuwlIVN99qtL7yvvtqw126rVtC/fzo2uPHcbv/+qQhLkqTtZgmWWtp778HUqRs/pk9PIw71+vX76Nzu3nunEQdJklR0lmCpWFasSOV208LbeGGSvL8AAA5pSURBVEeG3XZLK7vnn5+u9duPudeuJEklVZQSHEI4Dvg50Br4bYzxymL8XKksrV+fZnQ3LbtvvJEOngDo0CGt6J5wQiq79Y/u3bPNLkmSgCKU4BBCa+Aa4BjgHeClEMJ9McYZzf3ZUqZibBhlmDKloezOmNFwbHD93O6QITBqVEPZ7dcvvSZJkspSMVaChwGvxxhnAYQQ/gycCFiCVTlWr06jDJMmbVx4Fy1qeE+PHqngXnBBQ9ndbz/YYYfsckuSpO1SjBLcE5jT6Pt3gIOL8HOL6o477uCdd97JOobKwdq16ZCJ5cvTHO/y5bByZcMoQ+vW0LMnDBiQZnU7dkzXxjsyLFkC48alhyRJ2qpevXpx+umnZx1jI8UowZvbhT9+5E0hjAZGA/Tu3bsIv1bahhhTuW1cdpcvTyW4Xvv2qeR27ZquHTum5zxcQpKkqlaMEvwOsEej73sB8zZ9U4zxeuB6gKFDh36kJLe0cvu3DxXZsmUweXLDY9KkNN5QP7vbtm0aXRgyBA44oOHRpUu2uSVJUiaKUYJfAvqHEPoCc4EzgX8ows+VNu+992DCBJg4MT1efnnj44O7dk1l95/+KRXdIUPSCWseMCFJkuo0uwTHGGtDCBcCD5O2SLspxji92cmkGGHu3IayW1985zX6g4YBA2DYMBg9umGVd/fdHWeQJElbVZR9gmOMDwIPFuNnKadihLff3niFd+JEWLAgvd6qVVrNPeoo+PjH4cADU+nt1Cnb3JIkqSJ5YpxKL0aYMwdefBFeeqmh+C5Zkl5v3TodNPGZz6Sye+CBaYXXU9UkSVKRWILV8pYsSWX3xRcbHu+9l16rqUn77Z52WkPhHTw47dAgSZLUQizBKq7Vq9PuDI0L79/+1vD6PvvAscemOd5hw1Lhbdcuu7ySJCmXLMHafjGmXRmeeQaeey4V3smTYd269Pruu8PBB8OXv5wK78c/DrvskmViSZIkwBKsplizJm1H9uyzqfg++yy8+256rWPHVHQvv7xhlbdnz2zzSpIkbYElWFu2cGEquvWPl15KRRigXz8YMQI+8Yn0GDQofaBNkiSpAliC1WDuXHjqqfQYO7ZhlremBoYOhQsvTIX3kEPSqIMkSVKFsgTn2Zw5DYX3qafg9dfT8zvvDIcfDuecA4cemmZ53a1BkiRVEUtwnsybB4891lB6Z81Kz++yC3zyk3DBBXDEEWlPXkcbJElSFbMEV7OVK2HcOHjkkfSYXneadZcuqfRefHEqvfvvb+mVJEm5YgmuJhs2wKRJqfA++ij89a+wdm3ah/fww+Hss+GYY9LevK1aZZ1WkiQpM5bgSrd8eSq8998P//u/sGBBen7w4LTSe8wxqQDvsEO2OSVJksqIJbgSvf12Kr0PPABPPplWe3feGY4/HkaOTMW3R4+sU0qSJJUtS3ClmDkT7rgD7rwTpk5Nzw0YABddBJ/9bNrFoaYm24ySJEkVwhJcrmJMH2S7885UfmfMgBBS2b3qqlR8BwzIOqUkSVJFsgSXm1mz4JZb4M9/hldeScX38MPhl7+EU06Bv/u7rBNKkiRVPEtwOVi2LK32/v73aUeHENIWZhddlIqv872SJElFZQnOSozpQ23XXw/33gurV8M++8B//Rd88YvQu3fWCSVJkqqWJbjUlixJK76/+Q28+ip07gznngtf+hIMHZpWgSVJktSiLMGlMm0aXH013HorrFoFw4enMnz66e7hK0mSVGKW4JYUIzz9NPz4x+kgix12gFGj4Pzz4WMfyzqdJElSblmCW0KM8NBD8P3vwwsvQNeu8L3vwQUXpK8lSZKUKUtwsT3xBFxxBTz3HPTtC7/+NZx9NnTokHUySZIk1bEEF8uUKXDZZakE9+oF110HX/mKp7hJkiSVoVZZB6h4ixfDhRemGd9Jk+BnP4PXXoPRoy3AkiRJZcqV4O0VI/zud/D//l/a9uxrX0szwF26ZJ1MkiRJ22AJ3h7z5sF558GDD6YjjX/1Kxg8OOtUkiRJKpDjEE11220wcGA67e0Xv4CxYy3AkiRJFcYSXKi1a+Gii+DMM2HffWHy5PR9K/8rlCRJqjSOQxRiwQI46aS07dlll8GPfuSH3iRJkiqYJXhbXn8djjsuzQHfdht8/vNZJ5IkSVIzWYK3ZuLEVIA3bEj7/w4fnnUiSZIkFYEleEumTIERI6BTJ3j4Ydh776wTSZIkqUgswZvzt7/BMceko46ffDIdfyxJkqSq4dYGm1q6FE44IY1APP64BViSJKkKuRLc2Pr1aQu0WbPSDLAjEJIkSVXJEtzYT3+a5n9/85t0EpwkSZKqkuMQ9aZNg+98B049FUaPzjqNJEmSWpAlGKC2Fs4+G3beGa69FkLIOpEkSZJakOMQADfemPYEvv126NYt6zSSJElqYa4Ef/hhGoM4/HA47bSs00iSJKkEXAm+9lpYsADuv98xCEmSpJzI90rwmjXws5+lk+GGDcs6jSRJkkok3yvBt94K8+fDzTdnnUSSJEkllO+V4N/9DgYMSEckS5IkKTfyW4LffhvGjYNRo5wFliRJypn8luA//Sldv/jFbHNIkiSp5PJbgu+5B4YPh759s04iSZKkEmtWCQ4hfDeEMDeEMKnuMbJYwVrU0qUwYQJ8+tNZJ5EkSVIGirE7xNUxxp8W4eeUzrhxsGEDHHVU1kkkSZKUgXyOQzzxBLRvn8YhJEmSlDvFKMEXhhCmhBBuCiF0LsLPa3lPPAGHHgrt2mWdRJIkSRnYZgkOITwWQpi2mceJwLXAXsAQYD5w1VZ+zugQwvgQwviFCxcW7T9Ak61YAdOmpRIsSZKkXNrmTHCMcUQhPyiEcAPwwFZ+zvXA9QBDhw6NhQYsuunTIUYYMiSzCJIkScpWc3eH2L3RtycD05oXpwQmT07XwYOzzSFJkqTMNHd3iB+HEIYAEXgL+GqzE7W06dOhQwf3B5YkScqxZpXgGOOoYgUpmVmzoF8/aJXPjTEkSZKUxy3S3nwzlWBJkiTlVr5KcIypBDsKIUmSlGv5KsGLFqUt0izBkiRJuZavEjxvXrr26pVtDkmSJGUqXyV40aJ07do12xySJEnKVL5K8Pvvp6slWJIkKdfyVYLrV4J33TXbHJIkScpUvkpw/UqwJViSJCnX8lWCFy2CTp2gpibrJJIkScpQvkrw++87DyxJkqScleBFixyFkCRJEm2yDlBSffvCXntlnUKSJEkZy1cJvvbarBNIkiSpDORrHEKSJEnCEixJkqQcsgRLkiQpdyzBkiRJyh1LsCRJknLHEixJkqTcsQRLkiQpdyzBkiRJyh1LsCRJknLHEixJkqTcsQRLkiQpdyzBkiRJyh1LsCRJknLHEixJkqTcsQRLkiQpdyzBkiRJyh1LsCRJknLHEixJkqTcCTHG0v/SEBYCb5f8F0NXYFEGv1el5X3OB+9zPnifq5/3OB+yvM97xhi7bfpkJiU4KyGE8THGoVnnUMvyPueD9zkfvM/Vz3ucD+V4nx2HkCRJUu5YgiVJkpQ7eSvB12cdQCXhfc4H73M+eJ+rn/c4H8ruPudqJliSJEmC/K0ES5IkSdVZgkMIx4UQXg0hvB5C+NfNvN4uhHBb3esvhBD6lD6lmquA+/zPIYQZIYQpIYTHQwh7ZpFTzbOt+9zofaeFEGIIoaw+faxtK+QehxA+X/f38/QQwp9KnVHNV8A/s3uHEJ4MIbxc98/tkVnk1PYLIdwUQlgQQpi2hddDCOEXdf8bmBJCOLDUGRuruhIcQmgNXAMcD+wHfCGEsN8mbzsHWBJj/HvgauBHpU2p5irwPr8MDI0xDgbuBH5c2pRqrgLvMyGEnYCLgRdKm1DNVcg9DiH0B74JHBpjHAhcWvKgapYC/16+Arg9xvgx4Ezg16VNqSK4GThuK68fD/Sve4wGri1Bpi2quhIMDANejzHOijGuBf4MnLjJe04Efl/39Z3A0SGEUMKMar5t3ucY45MxxpV13z4P9CpxRjVfIX8/A/yA9C85q0sZTkVRyD0+D7gmxrgEIMa4oMQZ1XyF3OcIdKr7emdgXgnzqQhijOOAxVt5y4nAH2LyPLBLCGH30qT7qGoswT2BOY2+f6fuuc2+J8ZYCywDdi1JOhVLIfe5sXOAh1o0kVrCNu9zCOFjwB4xxgdKGUxFU8jfywOAASGEZ0IIz4cQtrbSpPJUyH3+LnBWCOEd4EHgotJEUwk19f+7W1SbrH5xC9rciu6mW2AU8h6Vt4LvYQjhLGAocESLJlJL2Op9DiG0Io00fblUgVR0hfy93Ib0x6dHkv5E5+kQwqAY49IWzqbiKeQ+fwG4OcZ4VQjhEOCWuvu8oeXjqUTKqn9V40rwO8Aejb7vxUf/SOX/3hNCaEP6Y5etLd+r/BRynwkhjAD+DfhcjHFNibKpeLZ1n3cCBgFjQwhvAcOB+/xwXEUp9J/Z98YY18UY3wReJZViVY5C7vM5wO0AMcbngPZA15KkU6kU9P/dpVKNJfgloH8IoW8IoS1puP6+Td5zH3B23denAU9EN0yuNNu8z3V/TH4dqQA7Q1iZtnqfY4zLYoxdY4x9Yox9SLPfn4sxjs8mrrZDIf/M/gvwKYAQQlfSeMSskqZUcxVyn2cDRwOEEPYlleCFJU2plnYf8KW6XSKGA8tijPOzClN14xAxxtoQwoXAw0Br4KYY4/QQwveB8THG+4AbSX/M8jppBfjM7BJrexR4n38CdATuqPvc4+wY4+cyC60mK/A+q4IVeI8fBj4dQpgBrAe+EWN8P7vUaqoC7/PlwA0hhMtIf0T+ZReoKksI4VbS2FLXutnufwdqAGKMvyHNeo8EXgdWAl/JJmniiXGSJEnKnWoch5AkSZK2yhIsSZKk3LEES5IkKXcswZIkScodS7AkSZJyxxIsSZKk3LEES5IkKXcswZIkScqd/w8MdMkURdSNiQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "x0 = 0.5\n", "iter=0 newton x=0.24205584583201645 err=2.5794e-01 bounds=(0, 1)\n", "iter=1 newton x=0.22186227635177588 err=2.0194e-02 bounds=(0, 1)\n", "iter=2 newton x=0.22209978959601762 err=2.3751e-04 bounds=(0, 1)\n", "iter=3 newton x=0.222099829644735 err=4.0049e-08 bounds=(0, 1)\n", "Converged in 4 iterations, x* = 0.2220998296447350\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAHSCAYAAADxDj0WAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3debhdZX0v8O/LFBDCoMxjEIIGMAqmAapivVgFrVBQLCqCilKUW63VWr3YVq/Xq7Uqrdo61LGtlKpIRVFQkGIuQiSIDEkYEmaFMIVowhCSrPvHm8PZgQwnOWfvdc7Zn8/zrGftvddir19cJvnmPb/1vqVpmgAAANVGbRcAAACjiYAMAAAdBGQAAOggIAMAQAcBGQAAOgjIAADQYZO2C+i0/fbbN5MmTWq7DAAAxrmrrrrq/qZpdljdsVEVkCdNmpRZs2a1XQYAAONcKeX2NR3TYgEAAB0EZAAA6CAgAwBABwEZAAA6CMgAANBBQAYAgA4CMgAAdBCQAQCgg4AMAAAdBGQAAOggIAMAQAcBGQAAOgjIAADQQUAGAIAOXQ/IpZQjSyk3llLmlVLe3+3rAQDAcHQ1IJdSNk7yT0mOSrJ/kteVUvbv5jUBAGA4uj2CPD3JvKZpbmmaZmmSs5Mc0+VrAgDABtuky9+/W5I7O97fleSQLl9zvX3729/OXXfd1XYZAAB9Zffdd8/xxx/fdhlP0e0R5LKaz5pVTijl1FLKrFLKrPvuu6/L5QAAwNp1ewT5riR7dLzfPclvOk9omuZLSb6UJNOmTVslPPfKaPyXCwAA7ej2CPKVSSaXUvYupWyW5IQk53X5mgAAsMG6OoLcNM2yUsr/THJhko2TfLVpmtndvCYAAAxHt1ss0jTND5P8sNvXAQCAkWAlPQAA6CAgAwBABwEZAAA6CMgAANBBQAYAgA4CMgAAdBCQAQCgg4AMAAAdBGQAAOggIAMAQAcBGQAAOgjIAADQQUAGAIAOAjIAAHQQkAEAoIOADAAAHQRkAADoICADAEAHARkAADoIyAAA0EFABgCADgIyAAB0EJABAKCDgAwAAB0EZAAA6CAgAwBABwEZAAA6CMgAANBBQAYAgA4CMgAAdBCQAQCgg4AMAAAdBGQAAOggIAMAQAcBGQAAOgjIAADQQUAGAIAOAjIAAHQQkAEAoIOADAAAHQRkAADoICADAEAHARkAADoIyAAA0EFABgCADgIyAAC99/a3J5/9bNtVrNYmbRcAAEAf+t73kmXL2q5itYwgAwDQe7/9bTJxYttVrJaADABAb61YkSxZkmy9dduVrJaADABAby1eXPdGkAEAILW9IhGQAQAgSfK739W9FgsAAIgRZAAAWIURZAAA6GAEGQAAOhhBBgCADgMB2QgyAAAkWbSo7o0gAwBAknvvraPHEya0XclqCcgAAPTWggXJTju1XcUaCcgAAPSWgAwAAB0WLEh23rntKtZIQAYAoLeMIAMAwEpLlyYPPiggAwBAkuSee+peiwUAACS57ba6nzSpzSrWSkAGAKB3BGQAAOgwEJD32qvVMtZGQAYAoHduvTXZdddRu4peIiADANBL8+Yl++zTdhVrJSADANAbTZPMnp1MmdJ2JWslIAMA0Bv33pssXJjsv3/blayVgAwAQG/MmVP3AjIAACS56qq6nzq13TrWoWsBuZTyoVLKr0spv1q5vaJb1wIAYAyYObPOfzyKl5lOkk26/P1nNk3zyS5fAwCAsWDmzOQFL2i7inXSYgEAQPfdfXdy553JIYe0Xck6dTsg/89SyrWllK+WUrbr8rUAABitZs6s+/EekEspF5VSrl/NdkySzyfZJ8nzktyd5FNr+I5TSymzSimz7rvvvuGUAwDAaHXZZcmmmyYHHdR2JetUmqbp/kVKmZTkB03THLi286ZNm9bMmjWr6/UAANBjz3lOsuOOycUXt11JkqSUclXTNNNWd6ybs1js0vH22CTXd+taAACMYnfemVx/ffKKsTGpWTdnsfhEKeV5SZoktyX50y5eCwCA0epHP6r7o45qt44h6lpAbprmjd36bgAAxpDzz0/23DOZMqXtSobENG8AAHTPgw8mF1yQHHdcUkrb1QyJgAwAQPd85zvJ0qXJiSe2XcmQCcgAAHTPv/978uxnJwcf3HYlQyYgAwDQHXPnJjNmJG9845hpr0gEZAAAuuVzn0smTEje9ra2K1kvAjIAACPvoYeSb3wjed3rkh12aLua9SIgAwAw8j7zmWTJkuRd72q7kvUmIAMAMLIeeig588zkmGOS5z2v7WrWm4AMAMDI+sQnakj+m79pu5INIiADADBy5s9PPvWpOnPFGJrarZOADADAyHnPe5JNN00+/vG2K9lgm7RdAAAA48R55yXf+14Nx7vu2nY1G8wIMgAAw3fffXW+4+c+N3n3u9uuZliMIAMAMDxNk5x2Wn0w76KLks02a7uiYRGQAQAYns99Lvnud+vsFc95TtvVDJsWCwAANtyMGclf/EVy9NH1Ab1xQEAGAGDDzJ+fvOY1yd57J//6r8lG4yNajo9fBQAAvbVgQfKylyXLlyff/36yzTZtVzRi9CADALB+fvvb5KijknvuSX760+RZz2q7ohElIAMAMHQLFyZHHplcd12d9/iQQ9quaMQJyAAADM1999W2ijlz6qwVRx3VdkVdISADALBu8+cnr3xlcvvtdeT45S9vu6KuEZABAFi7GTOSY4+tr3/yk+SFL2y3ni4ziwUAAGv29a8nRxyRbL99csUV4z4cJwIyAACr8/DDyVvekrz5zcnhhyeXX57su2/bVfWEgAwAwKrmzEmmT6+jxx/8YHLBBcl227VdVc8IyAAAVMuXJ5/8ZPL859cZKy68MPnIR5JN+uuxtf761QIAsHpz59Z2ipkzk2OOSb7whWTnnduuqhVGkAEA+tmjj9ZR4oMOSm6+OTnrrOTcc/s2HCdGkAEA+lPTJN//fvLudye33JIcf3zy2c8mO+3UdmWtM4IMANBvrr++LvpxzDHJhAnJRRcl3/qWcLySgAwA0C9uuy05+eRk6tTkssuST386ueaaOs8xT9BiAQAw3t1zT/Lxjyef/3yy0UbJe9+b/NVfJc94RtuVjUoCMgDAeHX77cknPpF85SvJ44/XhT/+9m+T3Xdvu7JRTUAGABhvbrihjhh/85tJKbWt4n3vSyZPbruyMUFABgAYD5qmPmz3mc8k55+fbL55cvrptZ3CiPF6EZABAMayxYuTf/u3OkXb3LnJDjvU5aH/7M/qa9abgAwAMBbNnp18+cvJ176WLFqUHHxw8o1vJH/yJ3XqNjaYgAwAMFb87nfJ2WfXh+5mzkw23TQ59tjkXe9KDjus9hszbAIyAMBotmJF8vOfJ1/9al3MY8mSZMqU5FOfSt74Rm0UXSAgAwCMRtddl5x1VvIf/1Gna9tyy+SEE5JTTkkOPdRocRcJyAAAo8Vtt9VAfNZZdTnojTdO/vAPk498JPnjP04mTmy7wr4gIAMAtGnevOS7363bzJn1s9///eRzn0uOPz7Zccd26+tDAjIAQC81TXLttcm559ZQfN119fNp05L/+39rG8Xee7dbY58TkAEAum3ZsuTyy5Pzzquh+JZbag/xi16U/MM/1PaJvfZqu0pWEpABALrhvvuSH/0o+eEPkwsvTB56qE7L9tKXJu9/f3L00clOO7VdJashIAMAjIQVK5Krr67LPJ9/fnLllbWdYqed6lzFr3xlDcfbbNN2payDgAwAsKHuvDO56KK6XXxxsmBBbZ2YPj358IeTV7wiOeigZKON2q6U9SAgAwAM1UMPJZdcMhiKb7qpfr7jjskRRyRHHZUceaTFO8Y4ARkAYE0efTS54orBQHzllbWVYsstkxe/ODnttNo2ceCBFu4YRwRkAIABixfX2SYuvTT52c/qvMRLl9YFOw45JPngB2sgPuSQZLPN2q6WLhGQAYD+tXBhctllNQxfemly1VXJ8uU1ED//+ck735kcfngdLd5667arpUcEZACgf/z613WEeMaMGoivvbbONLHZZvXBuve/vwbiww6zrHMfE5ABgPFp6dLkmmuSn/+8huLLL0/uuKMe22KLupzzhz9cA/H06fUziIAMAIwXCxbUEDwQiGfNqg/ZJckee9RR4Xe/u+4POkgPMWskIAMAY89jj9X2iF/8YjAQ33prPbbpprV/+O1vr2H4sMOS3Xdvt17GFAEZABjdli9P5s6tU6wNbNdckzz+eD2+yy41BJ9+et0ffHCy+ebt1syYJiADAKNH0yTz568ahn/5y+Thh+vxiROTadNqq8Tv/V7d9tzTHMSMKAEZAGhH0yR33VWnVhsIw7Nm1anXkmTChNorfMopg2F4v/0s20zXCcgAQPetWJHMm1dHg6++enC7//56fOON62p0r371YBg+8MDaTww9JiADACNr6dJkzpwagAcC8TXX1FXqkhp6DzwwOfroOkJ88MHJ856XPO1p7dYNKwnIAMCGW7Kkht+BEeFf/jKZPbuG5CTZcssaft/0psEwvP/+plhjVBOQAYB1a5rkzjvr1GoD2zXXJDfeWI8lyfbb1xD8538+GIb33VfPMGOOgAwArGrx4uT661cNw9demyxaNHjO3nsnU6cmJ5wwGIZ3281sEowLAjIA9KsVK+riGk8OwvPnD44KT5xYg/DrX1/3U6fW/uGtt263dugiARkA+sGiRcl1163aHnHddbWHOKkjv5Mn137hk05KnvvcGob32suoMH1HQAaA8WTJkrrq3OzZtU1iYH/nnYPnbLddDcCnnDI4KnzAAWaRgJUEZAAYix57LLnhhqcG4VtvHWyPmDAhefazk8MPrwF4YFRYrzCslYAMAKPZ448nN9/81CA8b16yfHk9Z5NN6gpz06YlJ59ce4QPOCDZZ596DFgvftcAwGiwfHlyyy1PDcI33lhDclKnS9tnnxqAjz9+MAjvt595hWEECcgA0EvLlye3315Xmps9ezAIz52bPPro4HmTJtUA/MpXDgbhZz872WKL1kqHfiEgA0A3LF1a2yDmzq1heGB/442rBuHddqsB+CUvGQzC+++fbLVVe7VDnxOQAWA4Hn64Piw3d+6qYXjevGTZssHz9tqrBt8jjqj7KVPqfttt26sdWC0BGQCG4qGHVg3AA69vv31w1oiNN65LK++/f3LccYMh+FnPSrbcst36gSEbVkAupRyf5ENJpiSZ3jTNrI5jH0hySpLlSd7ZNM2Fw7kWAHRd0yT33rtqS8RAGL777sHzBqZPO/TQ5C1vqUF4ypS60IaH5WDMG+4I8vVJjkvyxc4PSyn7JzkhyQFJdk1yUSllv6Zplg/zegAwfCtW1IUzntwfPHdusnDh4HkTJ9bg+/KXD44GT5lSH6DbeOPWyge6a1gBuWmauUlSnjrZ+DFJzm6a5rEkt5ZS5iWZnuTy4VwPANbLsmXJ/PlP7Q+eO7f2Dg/YYYcafF/72lX7g3fd1YIa0Ie61YO8W5IrOt7ftfIzABh5jz6a3HTTU/uDb765ziYxYPfda/B929tWHRHefvv2agdGnXUG5FLKRUl2Xs2hM5qm+d6a/rPVfNas4ftPTXJqkuy5557rKgeAfrZo0WAA7txuvbW2TSR1MY1nPrMG3z/6o8Eg/Oxn15YJgHVYZ0BumualG/C9dyXZo+P97kl+s4bv/1KSLyXJtGnTVhuiAegjTVMfiHtyCL7hhlUflNtss7qC3MEHJ294w+Bo8H77JZtv3l79wJjXrRaL85KcVUr5dOpDepOT/KJL1wJgLBpYWnkg/HaG4d/+dvC8rbde9UG5KVPqaPDeeyebmK0UGHnDnebt2CSfTbJDkvNLKb9qmublTdPMLqV8K8mcJMuSnG4GC4A+9cgjdfW4J4fgm25atT94l11q+H3jG2sAHgjDu+ziQTmgp0rTjJ6uhmnTpjWzZs1a94kAjD4LF66+P/i22wYX0thoozryOxB+O0eErSgH9FAp5aqmaaat7pifTQEwdE2T/PrXqw/C9947eN6ECXX1uOnTk5NPXnUhDf3BwCgnIAPwVE+eP3igPeKGG5Lf/W7wvG23XXW2iIHRYAtpAGOYgAzQzx5+uPYHP3k0+Oabk8cfHzxvt91q+O0cDZ4yJdlpJ/3BwLgjIAP0g4UL68IZnQtpzJ2b3H774DkbbZTss08Nvq961aojwltv3V7tAD0mIAOMJ/ffPxiEB7bZs5N77hk8Z4stan/w7/9+csopg0F4331r7zBAnxOQAcaapqkPxK0uCN933+B5W21VF8848si6H9j22quOFgOwWgIywGg1sKLck4PwnDnJAw8Mnrf11skBByRHHz0Ygg84INl9d/3BABtAQAZoW9Mkd921+iD80EOD5223XQ2+r371qkHYQhoAI0pABuiVFSuSO+5YfRDunDpthx1q+H3d61YNwjvuKAgD9ICADDDSmiZZsCC5/vrkuuvq/vrra4/wkiWD5+28cw2/J5+8ao/wDju0VzsAAjLAsCxcWIPvQAge2Dp7hHfYIXnOc+qMEQOjwVOmJM94Rnt1A7BGAjLAUCxZUucNfnIQ/vWvB8+ZODE58MDkuOPqfmDbccf26gZgvQnIAJ2WLk1uuumpQfiWW2rrRFLnCt5//+R//I9Vg/Aee+gRBhgHBGSgPzVNfWDu2mvrNtArfOONybJl9ZyNN0722y85+ODkpJMGg/A++9RjAIxLAjIw/i1ZUsPvtdcm11wzGIoXLRo8Z9Kk2id89NGDQfhZz7KyHEAfEpCB8aNpkttue2oQnjdvsD1i4sRk6tTk9a+v+6lTazCeOLHV0gEYPQRkYGz63e/qqHBnEL722sH5hEtJ9t23BuATT0ye+9z62jLLAKyDgAyMbgO9wldfnfzqV4NBeP78wXO22aaG35NOGgzCBxyQbLVVe3UDMGYJyMDosWxZnUHi6qsHt1/9KnnwwXp8o42SyZPrQ3NvfvNgi8See5o9AoARIyAD7XjkkTpzRGcYvu66+nlSH46bOjV5zWuS5z0vOeig+v5pT2u3bgDGPQEZ6L6FC+tIcGcYvuGGZPnyenybbWoAPu20uj/ooOTZz0428UcUAL3nbx9gZN1/f3LVVcmsWXX/y18mt98+eHzXXWsAPvbYwTA8aZIWCQBGDQEZ2HALFw6G4YGtMwxPnpwccsiqI8OWXQZglBOQgaH57W/raHBnGO6cSWKffWoYPv30ZNq0+iDdNtu0Vy8AbCABGXiqxYtrn3Dn6PCNNw4e32uvGoLf+tbBMPz0p7dXLwCMIAEZ+t3y5cncuckVVyQzZ9Zt9uxkxYp6fLfdagg+8cS6f/7zkx12aLdmAOgiARn6zd13DwbhmTOTK6+sI8ZJst12yfTp9QG63/u9GoZ32aXdegGgxwRkGM8eeaT2DXeODt9xRz22ySZ11bmTTkoOPbT2D0+ebDYJAPqegAzjRdMkt96aXHZZcvnlNQxfe21dnS6pfcOHHpq86111f9BByRZbtFszAIxCAjKMVY8/Xh+ku+yywe2ee+qxiRNri8Rf/mUNw9OnJzvv3G69ADBGCMgwVixcmPz854Nh+MorB5dlnjQpOeKI5AUvqNsBByQbb9xquQAwVgnIMBo1TZ1juHN0eM6cemzjjWt7xKmnDgbiXXdtt14AGEcEZBgNmia54Ybk0ksHt7vvrse22SY57LDkda+rYXj69GTLLdutFwDGMQEZ2rBiRZ1reCAM/+xnyb331mO77JK8+MV1e+ELk/33TzbaqN16AaCPCMjQC8uX1xklBgLxjBnJAw/UY3vskbzsZYOheN99TbUGAC0SkKEbmqb2DF98cd0uvTRZtKgee+Yzk6OPHgzEkya1WioAsCoBGUbKHXcMBuKLLx6ccm2ffZLXvraG4cMPryPGAMCoJSDDhnrggeSSSwYD8c0318933LFOuTawGSEGgDFFQIahWrq0zkN8wQXJT35SF+lommSrrero8DveUQPxgQfqIQaAMUxAhrW5/fYaiH/0ozpKvHhxsskmddq1D32oBuLp05NNN227UgBghAjI0OnRR+uUaz/6UQ3GN9xQP99rr+QNb0iOOip5yUuSrbdut04AoGsEZLjlluT882sgvuSSunzzhAm1beJP/zQ58sjkWc/SNgEAfUJApv+sWJH84hfJeefVbfbs+vl++yVve1sNxC9+cfK0p7VbJwDQCgGZ/rBkSX2w7vvfT37wg7pq3cYb1yD81rcmr3pVnY4NAOh7AjLj191310B83nnJRRcljz2WbLNN8opX1EB85JHJdtu1XSUAMMoIyIwvd9yRfPe7yXe+U6dka5pk772T006rq9e96EVmnAAA1kpAZuybPz8555y6/eIX9bOpU5MPfzg59tjkgAM8YAcADJmAzNh04411lPicc+qCHUkybVrysY8lr351Mnlyu/UBAGOWgMzYcccdydlnJ2edlVxzTf3ssMOST36yhmJLOgMAI0BAZnS7//46UnzWWcmMGfWzQw9N/uEfaijeffd26wMAxh0BmdFn8eI688RZZyUXXpgsW5ZMmZL8n/+TnHCC6dgAgK4SkBkdmqYu8fz1ryff/nadt3j33ZN3v7su8Tx1qgftAICeEJBp1x13JN/4Rg3Gt9ySTJxYR4lPOil54QuTjTZqu0IAoM8IyPTeI48k556bfO1rycUX19Hjl7wk+dCHkuOOS7bcsu0KAYA+JiDTO3PnJl/4Qh0xXrQo2Wuv5G/+Jjn55LqYBwDAKCAg011Ll9bR4s9/Prn00rqK3Wtek7z1rckf/IEWCgBg1BGQ6Y7bbku+9KXkK19J7r23jhB//OPJm9+c7Lhj29UBAKyRgMzIGZiJ4swz6zRtpSSvelVy2mnJy15mtBgAGBMEZIbvsceS//zPunjH1Vcn22+f/K//VYOxhTwAgDFGQGbD3Xdffejun/85ueeeZP/9k3/5lzpv8RZbtF0dAMAGEZBZf7ffnnzyk8mXv5w8+mhy1FHJn/958od/aDEPAGDME5AZuhtuSP7u75J///f6/qSTkve+ty4DDQAwTgjIrNsvf5l87GPJOeckm2+evOMdyXvek+y5Z9uVAQCMOAGZNbv66uSv/zo5//xk662TD3ygtlLssEPblQEAdI2AzFPNmVNXuDvnnGS77ZKPfjQ5/fRkm23argwAoOsEZAbNm5d8+MPJN7+ZbLVV8rd/m7z73YIxANBXBGSSBQtqGP7yl5PNNkv+8i/rtv32bVcGANBzAnI/e+SRurjHxz5WX592WnLGGckuu7RdGQBAawTkftQ0ydlnJ+9/f3LHHckxxySf+ESy335tVwYA0LqN2i6AHps1KznssOT1r0+e/vTkpz9N/uu/hGMAgJUE5H6xcGGdiWL69LoS3te+VsPyS17SdmUAAKOKFovxrmnqynfvfW9y//3JO99ZZ6owMwUAwGoJyOPZzTcnb3tbcumlyaGHJhdemDzveW1XBQAwqmmxGI+WL0/OPDOZOjW55prki19MLrtMOAYAGIJhBeRSyvGllNmllBWllGkdn08qpTxSSvnVyu0Lwy+VIbnxxuTww5O/+IvkpS9NZs9OTj012ci/hQAAhmK4LRbXJzkuyRdXc2x+0zSGLHtlxYo6p/EZZyRbbJH8278lb3hDUkrblQEAjCnDCshN08xNkiKEtevuu5OTTkouuig5+ujkC1+w2AcAwAbq5s/d9y6lXF1KubSU8qIuXqe//eAHtdf4ssuSf/mXOqexcAwAsMHWOYJcSrkoyc6rOXRG0zTfW8N/dneSPZumeaCU8vwk/1VKOaBpmt+u5vtPTXJqkuy5555Dr7zfPfZY8r73JZ/5TPLc5yb/8R/JlCltVwUAMOatMyA3TfPS9f3SpmkeS/LYytdXlVLmJ9kvyazVnPulJF9KkmnTpjXre62+dNddyWtek8ycmbzrXcnHP55svnnbVQEAjAtdmQe5lLJDkgebplleSnlmkslJbunGtfrOz36WHH988vDDyTnnJMcd13ZFAADjynCneTu2lHJXksOSnF9KuXDlocOTXFtKuSbJd5Kc1jTNg8Mrtc81TW2nOOKIZNtt6+ixcAwAMOKGO4vFuUnOXc3n5yQ5ZzjfTYfHH09OOy356lfrLBX/+q+WigYA6BKrR4x2ixYlr3hFDcd/8zfJuecKxwAAXdSVHmRGyB131HB8443J17+enHxy2xUBAIx7AvJodd11yctfXh/Gu+CC2nsMAEDXCcij0ZVX1nD8tKcl/+//JQce2HZFAAB9Qw/yaDNjxuBMFTNmCMcAAD0mII8mP/lJHTnebbcajvfeu+2KAAD6joA8Wvz0p8mrXpXst19y6aU1JAMA0HN6kEeDn/+8zm88eXJy8cXJM57RdkUAAH3LCHLbfvnL5Kijkl13rS0WwjEAQKsE5DbNn197jrfbro4c77xz2xUBAPQ9AbktDzxQFwFpmuTHP0722KPtigAAiB7kdjz6aPLHf5zcfnsdOd5vv7YrAgBgJQG515omOeWUugDI2WcnL3hB2xUBANBBi0WvnXlmctZZyUc/mvzJn7RdDQAATyIg99Kllybve19y3HHJBz7QdjUAAKyGgNwrv/lNHTHeZ5/ka19LSmm7IgAAVkMPci8sW5a89rXJ735XH8rbeuu2KwIAYA0E5F746EeTyy5LvvnN5IAD2q4GAIC10GLRbTNnJh/5SPKGNySvf33b1QAAsA4CcjctXpyceGKy227J5z7XdjUAAAyBFotuOuOMupz0JZck227bdjUAAAyBEeRu+cUvks9+Njn99OTFL267GgAAhkhA7obHH09OPTXZddf6gB4AAGOGFotuOPPM5Jprku9+15RuAABjjBHkkfab3yQf/nBy9NHJsce2XQ0AAOtJQB5pf/3XtcXi059uuxIAADaAgDySrrmmLiP9Z39Wl5QGAGDMEZBHStMk73lPst12yQc/2HY1AABsIA/pjZQLL0wuvjj5x3+sIRkAgDHJCPJIaJrkQx9K9torOe20tqsBAGAYjCCPhB//OJk5M/niF5PNNmu7GgAAhsEI8nA1TZ3Wbc89kze9qe1qAAAYJiPIw3Xppcnllyef/7zRYwCAccAI8nB9+tPJDjskJ5/cdiUAAIwAAXk4bms0aoUAAAzYSURBVLop+cEPkre/Pdlii7arAQBgBAjIw/GP/5hsumnyjne0XQkAACNEQN5QixYlX/968oY3JDvt1HY1AACMEAF5Q519dvLww7W9AgCAcUNA3lBf/nLynOck06a1XQkAACNIQN4Q11yTzJqVvPWtSSltVwMAwAgSkDfEV76STJiQnHhi25UAADDCBOT1tWxZctZZyTHHJE9/etvVAAAwwgTk9XXJJckDDySve13blQAA0AUC8vr61reSrbZKjjyy7UoAAOgCAXl9PP54cu65ydFHJ5tv3nY1AAB0gYC8PgbaK1772rYrAQCgSwTk9fGd7yQTJyYvf3nblQAA0CUC8lA1TfLDH9ZwrL0CAGDcEpCHavbs5Ne/9nAeAMA4JyAP1QUX1L32CgCAcU1AHqoLLkgOPDDZffe2KwEAoIsE5KFYvDiZMUN7BQBAHxCQh+K//ztZulRABgDoAwLyUPz0p8mECckLX9h2JQAAdJmAPBQzZiSHHFJDMgAA45qAvC6LFydXX5286EVtVwIAQA8IyOtyxRXJ8uUCMgBAnxCQ12XGjGSjjZLDDmu7EgAAekBAXpfLL0+mTk223rrtSgAA6AEBeW2aJrnqqmTatLYrAQCgRwTktbnttuTBBwVkAIA+IiCvzVVX1f3zn99uHQAA9IyAvDZXXZVsumnynOe0XQkAAD0iIK/NVVclBx5ogRAAgD4iIK/N1VcnBx/cdhUAAPSQgLwm992X3H9/HUEGAKBvCMhrMmdO3e+/f7t1AADQUwLymgjIAAB9SUBekzlzkokTk912a7sSAAB6SEBekzlz6uhxKW1XAgBADwnIazIQkAEA6CsC8uosXJjcc08yZUrblQAA0GMC8urMn1/3kye3WwcAAD0nIK/ObbfV/d57t1oGAAC9N6yAXEr5+1LKDaWUa0sp55ZStu049oFSyrxSyo2llJcPv9QeuvXWup80qdUyAADoveGOIP8kyYFN00xNclOSDyRJKWX/JCckOSDJkUn+uZSy8TCv1Tu33ZZst12yzTZtVwIAQI8NKyA3TfPjpmmWrXx7RZLdV74+JsnZTdM81jTNrUnmJZk+nGv11K23Gj0GAOhTI9mD/JYkP1r5erckd3Ycu2vlZ2PDrbfqPwYA6FPrDMillItKKdevZjum45wzkixL8s2Bj1bzVc0avv/UUsqsUsqs++67b0N+DSOraWqLhRFkAIC+tMm6Tmia5qVrO15KOTnJHyU5ommagRB8V5I9Ok7bPclv1vD9X0rypSSZNm3aakN0Ty1YkDz6qBFkAIA+NdxZLI5M8ldJjm6a5uGOQ+clOaGUMqGUsneSyUl+MZxr9Ywp3gAA+to6R5DX4XNJJiT5SSklSa5omua0pmlml1K+lWROauvF6U3TLB/mtXrj9tvrfq+92q0DAIBWDCsgN02z71qOfTTJR4fz/a24556633nndusAAKAVVtJ7sgULkk02SZ7+9LYrAQCgBQLyky1YkOy4Y7KR/2kAAPqRFPhkCxYkO+3UdhUAALREQH6ye+4RkAEA+piA/GRGkAEA+pqA3KlpknvvFZABAPqYgNzpoYeSpUtN8QYA0McE5E4LFtS9EWQAgL4lIHcSkAEA+p6A3ElABgDoewJyp4FlpgVkAIC+JSB3WrAg2Xjj5BnPaLsSAABaIiB3evDB5OlPt8w0AEAfkwQ7LVyYbLdd21UAANAiAbnTQw8l227bdhUAALRIQO4kIAMA9D0BudPChQIyAECfE5A7GUEGAOh7AvKApqkB2UN6AAB9TUAe8OijydKlRpABAPqcgDxg4cK6F5ABAPqagDzgoYfqfptt2q0DAIBWCcgDFi+u+4kT260DAIBWCcgDliyp+622arcOAABaJSAPGBhB3nLLdusAAKBVAvIAI8gAAERAHmQEGQCACMiDjCADABABedDACLKADADQ1wTkAUuWJJtskmy2WduVAADQIgF5wOLFRo8BABCQn7B4sQf0AAAQkJ+wZIkRZAAABOQnaLEAACAC8qAlS7RYAAAgID/hkUeSLbZouwoAAFomIA947LFkwoS2qwAAoGUC8gABGQCACMiDBGQAACIgD3rsMavoAQAgID9h6VIjyAAACMhP0GIBAEAE5EECMgAAEZCrpqktFnqQAQD6noCcJI8/XvdGkAEA+p6AnNT2ikRABgBAQE4iIAMA8AQBORGQAQB4goCc1Af0Eg/pAQAgICcxggwAwBME5ERABgDgCQJyIiADAPAEATkZDMh6kAEA+p6AnAw+pGcEGQCg7wnIiRYLAACeICAnAjIAAE8QkBMBGQCAJwjIiYVCAAB4goCcJNttl0yfnmy1VduVAADQsk3aLmBUeNWr6gYAQN8zggwAAB0EZAAA6CAgAwBABwEZAAA6CMgAANBBQAYAgA4CMgAAdBCQAQCgg4AMAAAdBGQAAOggIAMAQAcBGQAAOgjIAADQQUAGAIAOwwrIpZS/L6XcUEq5tpRybill25WfTyqlPFJK+dXK7QsjUy4AAHTXcEeQf5LkwKZppia5KckHOo7Nb5rmeSu304Z5HQAA6IlhBeSmaX7cNM2ylW+vSLL78EsCAID2jGQP8luS/Kjj/d6llKtLKZeWUl40gtcBAICu2WRdJ5RSLkqy82oOndE0zfdWnnNGkmVJvrny2N1J9mya5oFSyvOT/Fcp5YCmaX67mu8/NcmpK98uLqXcuAG/jpGwfZL7W7o2veM+j3/ucX9wn/uD+zz+tXmP91rTgdI0zbC+uZRycpLTkhzRNM3Dazjnv5O8t2maWcO6WBeVUmY1TTOt7TroLvd5/HOP+4P73B/c5/FvtN7j4c5icWSSv0pydGc4LqXsUErZeOXrZyaZnOSW4VwLAAB6YZ0tFuvwuSQTkvyklJIkV6ycseLwJP+7lLIsyfIkpzVN8+AwrwUAAF03rIDcNM2+a/j8nCTnDOe7W/CltgugJ9zn8c897g/uc39wn8e/UXmPh92DDAAA44mlpgEAoEPfBeRSypGllBtLKfNKKe9fzfEJpZT/XHl8ZillUu+rZDiGcI//opQyZ+US6ReXUtY4zQuj17ruc8d5rymlNKWUUfeUNOs2lPtcSnntyt/Ts0spZ/W6RoZnCH9m71lKuWTl2grXllJe0UadDE8p5aullHtLKdev4XgppXxm5f8Pri2lHNzrGjv1VUBeObPGPyU5Ksn+SV5XStn/SaedkmThyv7qM5P8XW+rZDiGeI+vTjJt5RLp30nyid5WyXAN8T6nlDIxyTuTzOxthYyEodznUsrkJB9I8oKmaQ5I8uc9L5QNNsTfyx9M8q2maQ5KckKSf+5tlYyQryc5ci3Hj0qd9Wxy6voYn+9BTWvUVwE5yfQk85qmuaVpmqVJzk5yzJPOOSbJN1a+/k6SI8rKKToYE9Z5j5umuaRjWkJLpI9NQ/m9nCQfSf0H0KO9LI4RM5T7/LYk/9Q0zcIkaZrm3h7XyPAM5R43SbZe+XqbJL/pYX2MkKZpfpZkbTOaHZPkX5vqiiTbllJ26U11T9VvAXm3JHd2vL9r5WerPadpmmVJFiV5Rk+qYyQM5R53OiWrLpHO2LDO+1xKOSjJHk3T/KCXhTGihvL7eb8k+5VSLiulXLFyfn7GjqHc4w8lObGUcleSHyb5s96URo+t79/fXTXceZDHmtWNBD95Go+hnMPoNeT7V0o5Mcm0JC/uakV0w1rvcyllo9QWqTf1qiC6Yii/nzdJ/ZHsH6T+NGhGKeXApmke6nJtjIyh3OPXJfl60zSfKqUcluTfVt7jFd0vjx4aVfmr30aQ70qyR8f73fPUH9U8cU4pZZPUH+dY5GTsGMo9TinlpUnOSF0F8rEe1cbIWdd9npjkwCT/XUq5LcmhSc7zoN6YM9Q/s7/XNM3jTdPcmuTG1MDM2DCUe3xKkm8lSdM0lyfZPMn2PamOXhrS39+90m8B+cokk0spe5dSNktt9j/vSeecl+Tkla9fk+Snjcmix5J13uOVP3r/Ymo41q84Nq31PjdNs6hpmu2bppnUNM2k1F7zo5ummdVOuWygofyZ/V9JXpIkpZTtU1subulplQzHUO7xHUmOSJJSypTUgHxfT6ukF85LctLK2SwOTbKoaZq72yqmr1osmqZZVkr5n0kuTLJxkq82TTO7lPK/k8xqmua8JF9J/fHNvNSR4xPaq5j1NcR7/PdJtkry7ZXPX97RNM3RrRXNehvifWaMG+J9vjDJy0opc5IsT/KXTdM80F7VrI8h3uP3JPmXUsq7U3/k/iYDV2NPKeU/Uluhtl/ZT/63STZNkqZpvpDaX/6KJPOSPJzkze1UWllJDwAAOvRbiwUAAKyVgAwAAB0EZAAA6CAgAwBABwEZAAA6CMgAANBBQAYAgA4CMgAAdPj/UaAjVvqeEX0AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "x0 = 0.5\n", "iter=0 bisection x=0.75 err=2.5000e-01 bounds=(0.5, 1)\n", "iter=1 bisection x=0.875 err=1.2500e-01 bounds=(0.75, 1)\n", "iter=2 bisection x=0.9375 err=6.2500e-02 bounds=(0.875, 1)\n", "iter=3 bisection x=0.96875 err=3.1250e-02 bounds=(0.9375, 1)\n", "iter=4 bisection x=0.984375 err=1.5625e-02 bounds=(0.96875, 1)\n", "iter=5 bisection x=0.9921875 err=7.8125e-03 bounds=(0.984375, 1)\n", "iter=6 bisection x=0.99609375 err=3.9062e-03 bounds=(0.9921875, 1)\n", "iter=7 bisection x=0.998046875 err=1.9531e-03 bounds=(0.99609375, 1)\n", "iter=8 bisection x=0.9990234375 err=9.7656e-04 bounds=(0.998046875, 1)\n", "iter=9 bisection x=0.99951171875 err=4.8828e-04 bounds=(0.9990234375, 1)\n", "iter=10 bisection x=0.999755859375 err=2.4414e-04 bounds=(0.99951171875, 1)\n", "iter=11 bisection x=0.9998779296875 err=1.2207e-04 bounds=(0.999755859375, 1)\n", "iter=12 newton x=0.9999986681276718 err=1.2074e-04 bounds=(0.999755859375, 1)\n", "iter=13 newton x=0.9999939680663965 err=4.7001e-06 bounds=(0.999755859375, 1)\n", "iter=14 newton x=0.9999817931722779 err=1.2175e-05 bounds=(0.999755859375, 1)\n", "iter=15 newton x=0.9999651586097119 err=1.6635e-05 bounds=(0.999755859375, 1)\n", "iter=16 newton x=0.9999559390072096 err=9.2196e-06 bounds=(0.999755859375, 1)\n", "iter=17 newton x=0.9999546240099589 err=1.3150e-06 bounds=(0.999755859375, 1)\n", "iter=18 newton x=0.999954604196403 err=1.9814e-08 bounds=(0.999755859375, 1)\n", "Converged in 19 iterations, x* = 0.9999546041964030\n", "x0 = 0.5\n", "iter=0 bisection x=0.25 err=2.5000e-01 bounds=(0, 0.5)\n", "iter=1 bisection x=0.125 err=1.2500e-01 bounds=(0, 0.25)\n", "iter=2 bisection x=0.0625 err=6.2500e-02 bounds=(0, 0.125)\n", "iter=3 bisection x=0.03125 err=3.1250e-02 bounds=(0, 0.0625)\n", "iter=4 bisection x=0.015625 err=1.5625e-02 bounds=(0, 0.03125)\n", "iter=5 bisection x=0.0078125 err=7.8125e-03 bounds=(0, 0.015625)\n", "iter=6 bisection x=0.00390625 err=3.9062e-03 bounds=(0, 0.0078125)\n", "iter=7 bisection x=0.001953125 err=1.9531e-03 bounds=(0, 0.00390625)\n", "iter=8 bisection x=0.0009765625 err=9.7656e-04 bounds=(0, 0.001953125)\n", "iter=9 bisection x=0.00048828125 err=4.8828e-04 bounds=(0, 0.0009765625)\n", "iter=10 bisection x=0.000244140625 err=2.4414e-04 bounds=(0, 0.00048828125)\n", "iter=11 bisection x=0.0001220703125 err=1.2207e-04 bounds=(0, 0.000244140625)\n", "iter=12 bisection x=6.103515625e-05 err=6.1035e-05 bounds=(0, 0.0001220703125)\n", "iter=13 bisection x=3.0517578125e-05 err=3.0518e-05 bounds=(0, 6.103515625e-05)\n", "iter=14 bisection x=1.52587890625e-05 err=1.5259e-05 bounds=(0, 3.0517578125e-05)\n", "iter=15 bisection x=7.62939453125e-06 err=7.6294e-06 bounds=(0, 1.52587890625e-05)\n", "iter=16 bisection x=3.814697265625e-06 err=3.8147e-06 bounds=(0, 7.62939453125e-06)\n", "iter=17 bisection x=1.9073486328125e-06 err=1.9073e-06 bounds=(0, 3.814697265625e-06)\n", "iter=18 bisection x=9.5367431640625e-07 err=9.5367e-07 bounds=(0, 1.9073486328125e-06)\n", "Converged in 19 iterations, x* = 0.0000009536743164\n", "x0 = 0.9\n", "iter=0 bisection x=0.95 err=5.0000e-02 bounds=(0.9, 1)\n", "iter=1 bisection x=0.975 err=2.5000e-02 bounds=(0.95, 1)\n", "iter=2 bisection x=0.9875 err=1.2500e-02 bounds=(0.975, 1)\n", "iter=3 newton x=0.9952833780711102 err=7.7834e-03 bounds=(0.975, 1)\n", "iter=4 newton x=0.9936312647131805 err=1.6521e-03 bounds=(0.975, 1)\n", "iter=5 newton x=0.9933150757720691 err=3.1619e-04 bounds=(0.975, 1)\n", "iter=6 newton x=0.9933071537399688 err=7.9220e-06 bounds=(0.975, 1)\n", "iter=7 newton x=0.9933071490757167 err=4.6643e-09 bounds=(0.975, 1)\n", "iter=8 newton x=0.9933071490757152 err=1.5543e-15 bounds=(0.975, 1)\n", "iter=9 newton x=0.9933071490757152 err=0.0000e+00 bounds=(0.975, 1)\n", "Converged in 10 iterations, x* = 0.9933071490757152\n" ] } ], "source": [ "def run(a,b,c,plot=False,**kwargs):\n", " '''Solves the equation with illustrations'''\n", " assert a*b<0, 'Must have different signs on a and b by the conditions of the problem'\n", " f = lambda x: a*np.log(x) + b*np.log(1-x) + c\n", " g = lambda x: a/x - b/(1-x)\n", " # plot\n", " if plot:\n", " xd = np.linspace(0,1,1000)\n", " plt.plot(xd,f(xd),c='r')\n", " plt.plot([0,1],[0,0],c='dimgrey')\n", " plt.show()\n", " def printiter (**kwargs):\n", " printiter.cout += 1\n", " iter = kwargs['iter']\n", " type = kwargs['type']\n", " x = kwargs['x1']\n", " bounds = kwargs['bounds']\n", " err = kwargs['err']\n", " if iter == 0:\n", " print('x0 = {}'.format(kwargs['x0']))\n", " print('iter={:<3d} {:<9s} x={:<23} err={:1.4e} bounds={:}'.format(iter,type,x,err,bounds))\n", " printiter.cout = 0\n", " xs = newton_bounds(f,g,bounds=(0,1),callback=printiter,**kwargs)\n", " print('Converged in %d iterations, x* = %1.16f' % (printiter.cout,xs))\n", "\n", "run(1,-4,0.5,plot=True)\n", "run(2,-1,-10,plot=True)\n", "run(2,-3,1e25)\n", "run(1,-1,-5,x0=.9,tol=1e-15)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Further learning resources\n", "\n", "- Oscar Veliz video on Newton-bisections hybrid (more general)\n", " [https://www.youtube.com/watch?v=FD3BPTMGJds&t=335s](https://www.youtube.com/watch?v=FD3BPTMGJds&t=335s) " ] } ], "metadata": { "celltoolbar": "Slideshow", "date": 1612589585.79215, "download_nb": false, "filename": "25_newton_bounds.rst", "filename_with_path": "25_newton_bounds", "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.6" }, "title": "Foundations of Computational Economics #25" }, "nbformat": 4, "nbformat_minor": 4 }