{ "metadata": { "name": "", "signature": "sha256:09c3fff5e022268c34b91e705a5b412603ed2b13fca15184bc5e11d8e0a97982" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "code", "collapsed": false, "input": [ "from __future__ import division\n", "\n", "import numpy as np\n", "import matplotlib.pyplot as pt\n", "\n", "from cmath import exp, pi" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 11 }, { "cell_type": "code", "collapsed": false, "input": [ "def approximate_stability_region_1d(step_function, make_k, prec=1e-5):\n", " def is_stable(k):\n", " y = 1\n", " for i in range(20):\n", " if abs(y) > 2:\n", " return False\n", " y = step_function(y, i, 1, lambda t, y: k*y)\n", " return True\n", " \n", " def refine(stable, unstable):\n", " assert is_stable(make_k(stable))\n", " assert not is_stable(make_k(unstable))\n", " while abs(stable-unstable) > prec:\n", " mid = (stable+unstable)/2\n", " if is_stable(make_k(mid)):\n", " stable = mid\n", " else:\n", " unstable = mid\n", " else:\n", " return stable\n", "\n", " mag = 1\n", " if is_stable(make_k(mag)):\n", " mag *= 2\n", " while is_stable(make_k(mag)):\n", " mag *= 2\n", "\n", " if mag > 2**8:\n", " return mag\n", " return refine(mag/2, mag)\n", " else:\n", " mag /= 2\n", " while not is_stable(make_k(mag)):\n", " mag /= 2\n", "\n", " if mag < prec:\n", " return mag\n", " return refine(mag, mag*2)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 18 }, { "cell_type": "code", "collapsed": false, "input": [ "def plot_stability_region(center, stepper):\n", " def make_k(mag):\n", " return center+mag*exp(1j*angle)\n", "\n", " stab_boundary = []\n", " for angle in np.linspace(0, 2*np.pi, 100):\n", " stable_mag = approximate_stability_region_1d(stepper, make_k)\n", " stab_boundary.append(make_k(stable_mag))\n", " \n", " stab_boundary = np.array(stab_boundary)\n", " pt.grid()\n", " pt.axis(\"equal\")\n", " pt.plot(stab_boundary.real, stab_boundary.imag)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 26 }, { "cell_type": "code", "collapsed": false, "input": [ "def fw_euler_step(y, t, h, f):\n", " return y + h * f(t, y)\n", "\n", "plot_stability_region(-1, fw_euler_step)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEACAYAAABI5zaHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtcFmXaB/DfI+AqaSClpECRAYHJAoqSGYkp+gpp1lqi\n76aoW+am5CnTavNQnrK2ttcyO4lZoqkZlEIpNWQeolS0whRdMUAhpcXVzFVh3j/uBUVOw3O655n5\nfT+f+eA8DM9cXQ3XM1xzzz0WVVVVEBGRobSQHQAREdkfizsRkQGxuBMRGRCLOxGRAbG4ExEZEIs7\nEZEB2Vzcx44dC19fX4SHh9f7fUVR4OXlhaioKERFReH555+3dZdERNQEd1vfYMyYMZg0aRJGjRrV\n4DZ9+vRBRkaGrbsiIiKNbD5zj42NRbt27RrdhvdJERE5l8N77haLBTt27EBERAQSEhKQn5/v6F0S\nEZmezW2ZpnTr1g1FRUXw9PREZmYmhg4dikOHDjl6t0RE5qbawdGjR9WuXbtq2jYwMFAtLy+v8/ot\nt9yiAuDChQsXLs1YIiIi6q21Dm/LlJWV1fTcc3NzoaoqfHx86mx35MgRqKpq8zJ79my7vI8ZFuaK\neWKuXD9P+/btq7f22tyWGTFiBHJycnDq1CkEBARg7ty5uHjxIgBg/PjxWL9+PZYtWwZ3d3d4enpi\nzZo1tu6yUYWFhQ59fyNhrrRhnrRjrrRxRp5sLu5paWmNfv+xxx7DY489ZutuiIioGQx3h2pycrLs\nEFwGc6UN86Qdc6WNM/JkUVVVdfheNLBYLNBJKERELqOh2mm4M3dFUWSH4DKYK22YJ+2YK22ckSfD\nFXciImJbhojIpZmmLUNERAYs7uz5acdcacM8acdcacOeOxERWYU9dyIiF8aeOxGRiRiuuLPnpx1z\npQ3zpB1zpQ177kREZBX23ImIXBh77kREJmK44s6en3bMlTbMk3bMlTbsuRMRkVXYcycicmHsuRMR\nmYjhijt7ftoxV9owT9oxV9qw505ERFZhz52IyIWx505EZCKGK+7s+WnHXGnDPGnHXGnDnjsREVmF\nPXciIhfGnjsRkYkYrriz56cdc6UN86Qdc6UNe+5ERGQVm3vuY8eOxaZNm9ChQwd8//339W6TkpKC\nzMxMeHp6IjU1FVFRUXUDYc+diKjZGqqd7ra+8ZgxYzBp0iSMGjWq3u9v3rwZhw8fRkFBAb755htM\nmDABu3btsnW3RHZz+jRQVASUlwNnzwJnzoivZ88CVVWAu3vtxcsLaN8euP56sVx3HeDmJvu/gqg2\nm4t7bGwsCgsLG/x+RkYGRo8eDQCIiYlBRUUFysrK4Ovra+uu66UoCuLi4hzy3kZjplydOQP88APw\n/ffA/v1AQYEo6MXFgKoCAQGiULdtC7Rpc3lxcwMKCxXccEMcLl0CLl4EKiqAU6cuLxUVgJ8fEBQE\nBAeLJSQE6N4d6NhR9n+5c5npmLKFM/Jkc3FvSklJCQICAmrW/f39UVxc7LDiTqSqwE8/AV99JZad\nO4GyMqBLFyA8XCyJiaKgBwQA114LWCwNv5+iAI39Hl66BPz8s/jAqF4++wzYvRto1Qro2VMst98O\n9OoF/OEP9v4vJqrL4cUdQJ1+kKWx3yQb8axBOyPlqrwc+PRT4JNPgJwccQZ+111A//7As8+Ks2pr\nWydN5cndHejcWSwDB15+XVWBo0eB3Fzg22+BJ58EDhwQcQ0cKJbg4MY/WFyNkY4pR3JGnhxe3P38\n/FBUVFSzXlxcDD8/v3q3TU5ORmBgIADA29sbkZGRNUmoHjrEda5Xr5eVASdOxOHjj4Fvv1XQvTsw\nblwcXn4ZOHJEfnzV6507AzfcoGDwYOCPf4zD1q3AypUK5s0DfHzikJQEdO6s4Kab9BEv1/W9rigK\nUlNTAaCmXtZLtYOjR4+qXbt2rfd7mzZtUgcNGqSqqqru3LlTjYmJqXc7O4Wifvnll3Z5HzNwxVyd\nO6eqH3ygqv37q6qPj6qOHauqGRnidUdxVJ6qqlQ1N1dVp0xR1U6dVDUiQlUXLlTV48cdsjuncMVj\nSgZ75qmh2mnzmfuIESOQk5ODU6dOISAgAHPnzsXFixcBAOPHj0dCQgI2b96MoKAgXHPNNVixYoWt\nuyQTys8Hli4F1q4FoqOBv/wFuPde0dN2VRYL0KOHWJYsAb7+Gnj/fXFtYOBAYOJEoHdvY7VtyHk4\ntwzplqqKC6JLlgDffQdMmACMGQPceKPsyByrogJYuRJ47TWgdWvg8ceBhx4CPDxkR0Z61FDtZHEn\n3VFV4OOPgYULxRj0adNEcWvdWnZkzlVVBWzdKj7cCgqAp54CkpOBli1lR0Z6YpqJw6ovPFDT9Jgr\nRRHDBefOBWbNEqNLHnlEbmGXlacWLYABA4AtW4APPgA2bBCja5YvF+Pt9UiPx5QeOSNPhivu5Jr2\n7QMGDQLGjgVSUoA9e4D77hMFjkTv/bPPgDVrgHXrgIgI4PPPZUdFesa2DEl1+jTw9NPA+vXi6/jx\nbDs0RVXFeP6pU8XF15deEmf0ZE6macuQa1BV4MMPRXG6cEGMhpk0iYVdC4sFGDIE+PFH4M47RRtr\n9myRR6Jqhivu7PlpJytXxcXi9v/nnhMF/s03AR8fKaFootdj6g9/AGbMEC2t3bvFkMq9e+XGpNdc\n6Q177mQ469eLCbXuuEP01Xv3lh2R6/PzE22a6dPF+Phnn+VZPLHnTk5y9qy4ULptmxj50bOn7IiM\n6cQJMbrol1/EX0U33SQ7InI09txJmv37gago0Sveu5eF3ZE6dgQyMoDhw4GYGGDTJtkRkSyGK+7s\n+WnnjFxt2AD06wfMmQO8846YI93VuNoxZbGIkTQbNoi7ep96SkxL7AyulitZ2HMnl1VVJXq/U6cC\nWVnA//6v7IjMp3dvcaE1N1fMw3P2rOyIyJnYcye7O3cOGDlSzLG+YQPQoYPsiMzt4kXgr38V8/N8\n+qm4AEvGwZ47OUVFhRix0bYtkJ3Nwq4HHh5iuOnw4WJM/P79siMiZzBccWfPTzt756q0FOjTRwx1\nXLnSODckGeGYsliAmTPFJGT9+wM7djhmP0bIlTOw504u4+hRcbfksGHAyy9zThi9Gj4cWLUKGDpU\nDEsl42LPnWz288/iuaDTp4sHTJD+bd0KjBghJiHjY09dG+dzJ4coLRWF/dFHxcgYch1ffCHO5Fng\nXZtpLqiy56edrbkqLwfi48WDNIxc2I16TN19t7iL9cEHxfw09mDUXNkbe+6kW7/9JuZfT0gAnnlG\ndjRkrb59xbNpExOBwkLZ0ZA9sS1DzVZVBTzwgLjbNDWVD3A2gn/8A1i2DNi+HbjuOtnRUHOw5052\n88wz4nF42dli2lkyhhkzgJ07xf9XowxjNQP23KkOa3L1/vvA6tXAxo3mKexmOaYWLQK8vcWoJ2uZ\nJVe2Ys+ddGXvXmDKFDHrYPv2sqMhe2vRQoyBz8wUH+Lk2tiWIU3OnBF3ns6bByQlyY6GHOn778VI\nmq1bxYO4Sd/YcyerqSowapRow7z9tuxoyBnS0sSsnnl5wDXXyI6GGsOeO9WhNVfvvSemjn31VcfG\no1dmPKZGjBCTjM2Y0byfM2OurMGeO0lXWAhMmwasXQt4esqOhpzp1VfFs1k/+0x2JGQNtmWoQaoq\nblTq0weYNUt2NCRDdjaQnCymCW7XTnY0VB/23KnZVq0CXnoJ+PZbMSc4mVNKirigvmKF7EioPg7r\nuWdlZSE0NBTBwcFYvHhxne8rigIvLy9ERUUhKioKzz//vK27bBR7fto1lqtffhHjnd9+m4Xd7MfU\n/PnAli3a5oA3e660ckae3G354crKSkycOBFbt26Fn58fevTogSFDhiAsLKzWdn369EFGRoZNgZJz\nzZghJgSLjpYdCcnWtq14yMdjj4lH9bm5yY6ItLDpzD03NxdBQUEIDAyEh4cHkpKSkJ6eXmc7Z7Zb\n4jh3qWYN5WrPHnERbfZs58ajVzymxL0NXl7AG280vh1zpY0z8mRTcS8pKUFAQEDNur+/P0pKSmpt\nY7FYsGPHDkRERCAhIQH5+fm27JIcTFXF6Jg5c8QZGxEgJod77TVg7lzg119lR0Na2NSWsWiYDrBb\nt24oKiqCp6cnMjMzMXToUBw6dKjebZOTkxEYGAgA8Pb2RmRkZM0nXHWPqqn16te0bm/m9by8PEye\nPLnW98+ejcMvvwBBQQoURV/xylq/+tiSHY/M9fvui8OSJcDAgfV/v/o1vcSr1/VXXnnFqvpWndvU\n1FQAqKmX9VJtsHPnTnXgwIE16wsWLFAXLVrU6M8EBgaq5eXldV63MZQaX375pV3exwyuztWlS6oa\nFqaqmzbJiUeveExd9vPPqurjo6qlpfV/n7nSxp55aqh22jQU8tKlS7j11luRnZ2NTp06oWfPnkhL\nS6t1QbWsrAwdOnSAxWJBbm4uHnzwQRTW81QADoWUb+1acePK119zjnZq2OOPi0nGXn5ZdiQENFw7\nbWrLuLu7Y+nSpRg4cCAqKysxbtw4hIWFYfny5QCA8ePHY/369Vi2bBnc3d3h6emJNWvW2LJLchBV\nBRYuFMPeWNipMbNmAV26iGsz/v6yo6GGGO4mJkVRavpU1Lgrc7V5s/ilzctjcb8aj6m6pkwR9z+8\n8ELt15krbeyZJ9NMHEbNp6rijP2pp1jYSZuUFODdd8WzdEmfDHfmTs2XmyvGMRcU8AYV0u7++4H4\neGDCBNmRmBvP3KlBb74JPPIICzs1z+TJ4sHaVVWyI6H6GK64XznelhqnKApOnwY2bADGjJEdjX7x\nmKpfbCzQurWYObIac6WNM/JkuOJOzbN6NdC/P+DrKzsScjUWCzB2rHiYC+kPe+4m160bsHix6J0S\nNdfJk0BwMFBcDLRpIzsac2LPneo4eBAoLQX69ZMdCbmq9u1Fe+ajj2RHQlczXHFnz0+7JUsU/OlP\n4m5DahiPqcaNGnW5NcNcacOeOzlUTg4wbJjsKMjVDR4sntZVXi47EroSe+4mdeQI0Ls3UFLCIZBk\nu3vvBYYPB0aOlB2J+bDnTrV8/LH4hWRhJ3tITAQ2bZIdBV3JcMWdPT9tsrOBTp0U2WG4BB5TTUtI\nEE/vys5WZIfiEthzJ4e4cEFM6xsZKTsSMgp/f7HwQWv6wZ67CW3fDkyaJJ6VSmQvU6cCHToAM2fK\njsRc2HOnGtnZwN13y46CjOaOO4AdO2RHQdUMV9zZH23a9u1Anz7MlVbMkza9egE5OQr4B3jT2HMn\nu1NVYO9eMe0AkT35+QGtWompo0k+9txNprhYFPayMj6Yg+wvKUmMnBk1SnYk5sGeOwG4fNbOwk6O\nEB4O/Pij7CgIMGBxZ3+0cXv2AFFR4t/MlTbMk3aVlQoOHJAdhf6x5052d/CgeHI9kSPcdBNY3HWC\nPXeT6dULWLIEuPNO2ZGQEV28CLRtC1RUiIur5HjsuRMAoLAQuPlm2VGQUXl4iLP3f/5TdiRkuOLO\n/mjDfv8d+Ne/gI4dxTpzpQ3zpJ2iKPDzA06ckB2JvrHnTnZ17BgQEMCHc5Bj3XADi7sesOduIl99\nBTz9NLBtm+xIyMimTRN/HU6fLjsSc2DPnVBRAXh7y46CjI5n7vpguOLO/mjDri7uzJU2zJN2iqKg\nfXvg1CnZkeibS/Tcs7KyEBoaiuDgYCxevLjebVJSUhAcHIyIiAjs3bvX1l2Slf71L6BdO9lRkNG1\nbi0u3pNcNhX3yspKTJw4EVlZWcjPz0daWhoOXHUHw+bNm3H48GEUFBTgzTffxIQJE2wKuClxcXEO\nfX9Xdvo04OV1eZ250oZ50i4uLg6tWgHnz8uORN+ccUzZVNxzc3MRFBSEwMBAeHh4ICkpCenp6bW2\nycjIwOjRowEAMTExqKioQFlZmS27JStduiTGIRM5Eou7PthU3EtKShAQEFCz7u/vj5KSkia3KS4u\ntmW3jWJ/VDvmShvmSTtFUfD778CWLbIj0TdnHFPutvywRePUglcP02no55KTkxEYGAgA8Pb2RmRk\nZM2fL9XJaGq9mtbtzbReWAgEBV1ez8vL01V8XHf9dQAoLwcABYoiPx69rufl5Vn984qiIDU1FQBq\n6mV9bBrnvmvXLsyZMwdZWVkAgIULF6JFixZ48skna7Z59NFHERcXh6SkJABAaGgocnJy4OvrWzsQ\njnN3uGefBdzdxVciR8nOBhYsEF/J8Rwyzj06OhoFBQUoLCzEhQsXsHbtWgwZMqTWNkOGDMF7770H\nQHwYeHt71yns5Dz8/CRH+/13ThqmBzYVd3d3dyxduhQDBw5Ely5dMHz4cISFhWH58uVYvnw5ACAh\nIQGdO3dGUFAQxo8fj9dff90ugTfkyj8Pqba2bYEzZy6vM1faME/aKYqC8+dZ3JvijGPKpp47AAwa\nNAiDBg2q9dr48eNrrS9dutTW3ZAdeHuL+dyJHInFXR84t4yJrFsHrF0LrF8vOxIystdeA374AVi2\nTHYk5sC5ZQjt2om7VIkc6cSJy9NKkzyGK+7sjzasXbvqYWoCc6UN86SdoigoLRWTh1HDnHFMGa64\nU8NuvBH4+WfZUZDRlZbyzF0P2HM3EVUVI2ZKSmrPMUNkT927A2+8AfToITsSc2DPnWCxAIGB4jmq\nRI6gquL5qY3cOElOYrjizv5o426++XJxZ660YZ6027hRgZsb0L697Ej0jT13srtbbgEKCmRHQUZ1\n7BgQFiY7CgLYczedFSvEnB/vvy87EjKi118H9u0D/nuDOjkBe+4EAIiKAvbskR0FGVV+Ps/c9cJw\nxZ390cZ16SJ67r/9xlxpxTxpl52toHt32VHoH3vuZHctW4ozq/37ZUdCRvP778DRo0B0tOxICGDP\n3ZQmTACCg4GpU2VHQkaybRswfTrwzTeyIzEX9typRt++wBdfyI6CjGbHDuCOO2RHQdUMV9zZH21a\n377iLGvrVkV2KC6Bx5Q227YBXl6K7DBcAnvu5BDt24s7CDm3O9nL+fPAV18B3brJjoSqseduUlOn\nAj4+wDPPyI6EjCArC5g/X5y9k3Ox5061JCQAGRmyoyCj2LQJSEyUHQVdyXDFnf1Rbfr0AQ4eVHDs\nmOxI9I/HVONU9XJxZ660Yc+dHMbDA7jzTuCjj2RHQq5u3z5R4Lt2lR0JXYk9dxPLygKeew7Yvl12\nJOTKpk0DPD3FsUTO11DtZHE3sQsXgE6dgN27gZtukh0NuaJLl4CAACAnBwgJkR2NOZnmgip7ftrt\n2KFg5Ejg3XdlR6JvPKYatnWreHxjdWFnrrRhz50c7pFHgHfeEWdgRM313nvAQw/JjoLqw7YMoXdv\n4MkngSFDZEdCrqS0VExCd+SIuGeC5DBNW4aab/x4PlyBmm/ZMmDECBZ2vTJccWfPT7vqXD3wAPDd\nd8CBA3Lj0SseU3WdPw+88QaQklL7deZKG/bcySlatxa/pIsXy46EXMXq1UD37kBoqOxIqCFW99x/\n/fVXDB8+HMeOHUNgYCA+/PBDeHt719kuMDAQ1157Ldzc3ODh4YHc3Nz6A2HPXaqKCvHw7N27xaRi\nRA2prBQ3LL36KhAfLzsasnvPfdGiRYiPj8ehQ4fQr18/LFq0qMEdK4qCvXv3NljYST5vb+Dhh4EX\nX5QdCend6tXAddcB/fvLjoQaY3Vxz8jIwOjRowEAo0ePxscff9zgts48I2fPT7urczVlivjFLS6W\nE49e8Zi67OJFYM4c4PnnAYul7veZK2103XMvKyuDr68vAMDX1xdlZWX1bmexWNC/f39ER0fjrbfe\nsnZ35AS+vuIRfJwGmBqyciVw881AXJzsSKgp7o19Mz4+HqWlpXVenz9/fq11i8UCS30f4wC2b9+O\njh074uTJk4iPj0doaChiY2Pr3TY5ORmB/234ent7IzIyEnH/PYqqP+m4bt/1atXrTz4Zh5AQ4K23\nFAQHy49PD+txcXG6ikfW+vnzwLx5cfjwQ33E48rr1a9Z8/OKoiA1NRUAauplfay+oBoaGgpFUXDD\nDTfgxIkT6Nu3L3766adGf2bu3Llo06YNpk2bVjcQXlDVjWXLgPXrxa3lDXxmkwk9+yxw6BCwZo3s\nSOhKdr+gOmTIEKxcuRIAsHLlSgwdOrTONufOncOZM2cAAL/99hs+//xzhIeHW7tLTa4+I6WGNZSr\nhx8Gjh8HPvnEufHoFY8p4PBh4PXXm77gzlxp44w8WV3cZ86ciS1btiAkJARffPEFZs6cCQA4fvw4\nEv/7SJbS0lLExsYiMjISMTExuOeeezBgwAD7RE4O4+4OLF0KTJoE/PezmUxMVYHHHwdmzAD8/WVH\nQ1pxbhlq0NixQJs2YjwzmVd6OjBzpngoR8uWsqOhq3E+d2q2X38VN6ts2AD06iU7GpLh1CkgIgJI\nSwPuukt2NFQf00wcxp6fdk3lyscH+Mc/gHHjxFwiZmXWY0pVgUcfBUaO1F7YzZqr5tJ1z53MYdgw\ncfY+fbrsSMjZVq8GfvqJj89zVWzLUJMqKoBu3cRIifvvlx0NOcPPPwPR0cBnnwFRUbKjocaw5042\n+eYbYPBg4Ntv+bxVozt/XrRhHnyQf7G5AvbcqY7m5ComBnjiCfFwhv/8x3Ex6ZHZjqmUFPEBXs+9\nhk0yW66sxZ476cq0aZfnn+EfWcb09tvA11+Lh6bz7mTXxrYMNcvZs0BsrBhB8cQTsqMhe9q5E7j3\nXmDbNuDWW2VHQ1o1VDsbnTiM6Gpt2gAZGcDtt4sCwIdqG8PBg8B99wGpqSzsRmG4tgx7ftpZm6uA\nAGDjRjH+3QzPXzH6MXXiBPA//wMsWgQkJNj2XkbPlb2w50661bMnsGKFGEGzb5/saMha//63KOh/\n+QuQnCw7GrIn9tzJJh9+CEyeDHzxBR+W7GrOnAESE4HwcDFRHC+guib23MkhHnwQOHdOPCg5Jwfo\n3Fl2RKTFv/8NDBoE3HYb8H//x8JuRIZry7Dnp529cpWcDDz1FNCnD/Djj3Z5S10x2jF1+jQwYICY\nEOyNN4AWdqwCRsuVozgjTzxzJ7uYMAG49lqgXz8xRWxMjOyIqD4nT4pWzO23i0nheMZuXOy5k11t\n2iTO5FevFq0a0o9Dh8TF0+HDgeefZ2E3CtNMP0ByJSaKYZJ//jPwzjuyo6Fq27aJm89mzgTmz2dh\nNwPDFXf2/LRzVK7uvBP46ivghRfEPCUXLzpkN07j6sdUWhrwpz8Bq1aJIY+O5Oq5chaOcyeXdeut\nYibJw4eBgQPFE33Iuf7zH/Ec3GeeAbZuFRdRyTzYcyeHqqwEnn4aWLMG+OADoHdv2RGZw7FjwAMP\nAH5+4mYzb2/ZEZGjsOdOUri5idvaX31VtAZmzwYuXZIdlbFlZIg7iJOSgI8+YmE3K8MVd/b8tHNm\nroYMAfbuBXbtEhf2/vlPp+3aZq5yTP36KzBqFDBlirioPXWq8y+cukquZGPPnQylY0cgM1OcUfbs\nCbz0Es/i7eWTT8Q0Au3aAfv3A3fcITsiko09d5KioAD461/FTTXLl/OmJ2sVFwMzZojZOd99Vzwe\nj8yFPXfSleBg4PPPRWEaOvRyoSdtfv8deO45IDJSzOezbx8LO9VmuOLOnp92snNlsYgnOuXniwuv\noaGiYJ09KzWsOmTn6UpVVWImzrAw0X759ltxt+k118iOTNBTrvSMPXcyhXbtxMyEubnAgQPirP61\n14ALF2RHph9VVcD69UBUlLg5LDUVWLcOuPlm2ZGRXrHnTrqTlydmmdy/X9zhOn484OUlOyo5qqqA\nDRuAefOAVq3EUNLERE4fQJfZvee+bt063HbbbXBzc8OePXsa3C4rKwuhoaEIDg7G4sWLrd0dmUhk\nJLB5M/Dpp6LAd+4MTJ8OFBXJjsx5ysuBJUuAoCAxquiFF8RfNvfcw8JO2lhd3MPDw7Fx40bc1chV\nnMrKSkycOBFZWVnIz89HWloaDhw4YO0uNWHPTzu95yoyEnj/fTE+vqpKzD+emCjOZJ3ZsnFWnlQV\n+O47YOxYUdR/+EHc2btrl3iwhisUdb0fU3qh6557aGgoQkJCGt0mNzcXQUFBCAwMhIeHB5KSkpCe\nnm7tLsmkbrwR+PvfxbC/pCTxSDh/f/F4v9xcUfhd2U8/iXZLaKiYjjckREzPu3KluB+AyBoOvaBa\nUlKCgICAmnV/f3+UlJQ4cpeIi4tz6PsbiavlytMTeOgh4MsvgZ07gbZtxdzx/v7iYSFZWWKyLHuz\nd54qK8WkavPmib9O7r5bPPZu1Sox0drMmUD79nbdpdO42jElizPy1OiTmOLj41FaWlrn9QULFmDw\n4MFNvrnFFf6OJJd0yy1i2ORzz4mz3PR08e/hw8XdmXfdJZYePYCWLeXGWlUlYty+HfjsMyA7G+jU\nSczS+MorYjoGNze5MZLxNFrct2zZYtOb+/n5oeiKq2BFRUXw9/dvcPvk5GQEBgYCALy9vREZGVnz\nCVfdo2pqvfo1rdubeT0vLw+TJ0/WTTzWroeEAD16KOjRA+jaNQ7btgEffKBgxQrgxIk4dOsG+Pgo\n6NwZGDYsDl27Art3a3//q4+txraPiYnDkSPAhg0KDh4EysrisHs30Lq1gi5dgD//OQ4vvwwUFOgn\nf/Zcr35NL/Hodf2VV16xqr5V5zY1NRUAauplfWweCtm3b1+8+OKL6N69e53vXbp0Cbfeeiuys7PR\nqVMn9OzZE2lpaQgLC6sbiJ2GQiqKUpMQapwZcnX6tOjLf//95SU/H/DxAQICxOLvL762bw+0aSPa\nPW3aiMXNDdi5U0FUVBwuXRIPHqmoEPPTnzwpvv7yi5gI7fBh8e+bbxZj9bt3Fz3z6GjXbbM0lxmO\nKXuwZ54aqp1WF/eNGzciJSUFp06dgpeXF6KiopCZmYnjx4/j4YcfxqZNmwAAmZmZmDx5MiorKzFu\n3DjMmjWrWQES2VtlpRhWWVwslqIisZSXi7tjz5y5/FVVAXf32ouXlyjW118vlvbtLxf0G29ki4Wc\ny+7F3d5Y3ImIms80E4dd2fujxjFX2jBP2jFX2jgjT4Yr7kRExLYMEZFLM01bhoiIDFjc2fPTjrnS\nhnnSjrlvOShIAAAEx0lEQVTShj13IiKyCnvuREQujD13IiITMVxxZ89PO+ZKG+ZJO+ZKG/bciYjI\nKuy5ExG5MPbciYhMxHDFnT0/7ZgrbZgn7ZgrbdhzJyIiq7DnTkTkwthzJyIyEcMVd/b8tGOutGGe\ntGOutGHPnYiIrMKeOxGRC2PPnYjIRAxX3Nnz04650oZ50o650oY9dyIisgp77kRELow9dyIiEzFc\ncWfPTzvmShvmSTvmShv23ImIyCrsuRMRuTD23ImITMTq4r5u3TrcdtttcHNzw549exrcLjAwEH/8\n4x8RFRWFnj17Wrs7zdjz04650oZ50o650kbXPffw8HBs3LgRd911V6PbWSwWKIqCvXv3Ijc319rd\naZaXl+fwfRgFc6UN86Qdc6WNM/Lkbu0PhoaGat7Wmb30iooKp+3L1TFX2jBP2jFX2jgjTw7vuVss\nFvTv3x/R0dF46623HL07IiJCE2fu8fHxKC0trfP6ggULMHjwYE072L59Ozp27IiTJ08iPj4eoaGh\niI2NtS5aDQoLCx323kbDXGnDPGnHXGnjlDypNoqLi1N3796tads5c+aoL774Yr3fi4iIUAFw4cKF\nC5dmLBEREfXWVKt77ldSG+ipnzt3DpWVlWjbti1+++03fP7555g9e3a92/JCDBGR/Vjdc9+4cSMC\nAgKwa9cuJCYmYtCgQQCA48ePIzExEQBQWlqK2NhYREZGIiYmBvfccw8GDBhgn8iJiKhBurlDlYiI\n7MeQd6j+7W9/Q0REBCIjI9GvXz8UFRXJDkmXnnjiCYSFhSEiIgL3338/Tp8+LTsk3dJ6055ZZWVl\nITQ0FMHBwVi8eLHscHRr7Nix8PX1RXh4uMP3ZcjiPmPGDOzbtw95eXkYOnQo5s6dKzskXRowYAB+\n/PFH7Nu3DyEhIVi4cKHskHRL6017ZlRZWYmJEyciKysL+fn5SEtLw4EDB2SHpUtjxoxBVlaWU/Zl\nyOLetm3bmn+fPXsW119/vcRo9Cs+Ph4tWohDICYmBsXFxZIj0q/Q0FCEhITIDkOXcnNzERQUhMDA\nQHh4eCApKQnp6emyw9Kl2NhYtGvXzin7sstoGT16+umnsWrVKnh6emLXrl2yw9G9d999FyNGjJAd\nBrmgkpISBAQE1Kz7+/vjm2++kRgRAS5c3Ju6wWr+/PmYP38+Fi1ahClTpmDFihUSopRPy41o8+fP\nR8uWLTFy5Ehnh6cr9rhpz4wsFovsEKgeLlvct2zZomm7kSNHIiEhwcHR6FdTeUpNTcXmzZuRnZ3t\npIj0S+sxRbX5+fnVGrRQVFQEf39/iRERYNCee0FBQc2/09PTERUVJTEa/crKysKSJUuQnp6OVq1a\nyQ7HZXD0cG3R0dEoKChAYWEhLly4gLVr12LIkCGywzI9Q45zHzZsGA4ePAg3NzfccsstWLZsGTp0\n6CA7LN0JDg7GhQsX4OPjAwDo1asXXn/9dclR6dPGjRuRkpKCU6dOwcvLC1FRUcjMzJQdlm5kZmZi\n8uTJqKysxLhx4zBr1izZIenSiBEjkJOTg/LycnTo0AHz5s3DmDFjHLIvQxZ3IiKzM2RbhojI7Fjc\niYgMiMWdiMiAWNyJiAyIxZ2IyIBY3ImIDIjFnYjIgFjciYgM6P8BD9w3cks19L4AAAAASUVORK5C\nYII=\n", "text": [ "" ] } ], "prompt_number": 27 }, { "cell_type": "code", "collapsed": false, "input": [ "def heun_step(y, t, h, f):\n", " yp1_fw_euler = y + h * f(t, y)\n", " return y + 0.5*h*(f(t, y) + f(t+h, yp1_fw_euler))\n", "\n", "plot_stability_region(-1, heun_step)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEACAYAAAC9Gb03AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl0VFW2wOFfZQCEIAGRMCQYIMEECKkIGOSJhpaASSCi\ntgJqSx62HWZQluDrbhVfy+TYauR1dDUg0h0RuxVETAPahaiEKARRGQSaQAaIA6LMQ6j3x5EYSAiV\nmk7dk/2tVYvcqhvv3uvEXbd2nXuuzel0OhFCCGGsIN0BCCGE8C0p9EIIYTgp9EIIYTgp9EIIYTgp\n9EIIYTgp9EIIYTiPCn1JSQkDBgyge/fu9OjRgxdeeKHW/SZNmkRsbCyJiYkUFRV5ckghhBD1FOLJ\nL4eGhvLcc89ht9s5cuQIvXr1IjU1lfj4+Kp9Vq5cya5du9i5cycbNmxg7NixFBQUeBy4EEII13h0\nRt+2bVvsdjsAYWFhxMfHU15eft4+y5cvZ9SoUQAkJydz6NAhKioqPDmsEEKIevBaj764uJiioiKS\nk5PPe76srIyoqKiq7cjISEpLS711WCGEEJfglUJ/5MgRfv3rX/P8888TFhZW4/ULV1mw2WzeOKwQ\nQggXeNSjBzh9+jS3334799xzD8OGDavxeocOHSgpKanaLi0tpUOHDrXud2HbRwghRN26dOnCrl27\n6tzHozN6p9PJfffdR7du3ZgyZUqt+2RmZrJo0SIACgoKCA8PJyIiosZ+5eXlOJ1OYx+PPfaY9hgk\nN8lP8jPvsXv37kvWao/O6D/++GMWL15Mz549SUpKAmDWrFns27cPgOzsbNLT01m5ciUxMTE0a9aM\nBQsWeHJIyyouLtYdgs+YnBtIflZnen6u8KjQX3/99Zw9e/aS++Xk5HhyGCGEEB6QK2P9JCsrS3cI\nPmNybiD5WZ3p+bnC5nQ6A+LGIzabjQAJRQghLMOV2iln9H7icDh0h+AzJucGkp/VmZ6fK6TQCyGE\n4aR1I4QQFiatGyGEEFLo/cXkPqHJuYHkZ3Wm5+cKKfRCCGE46dELUYtTp+C77yAkBJo0UY/QUJD1\n+ESgcaV2eryomRBW5HTCnj2wfj1s3gzl5XDgwC+Pw4ehVSs4exZOnFCPM2egcWO47DLo2BHi4s5/\ndO0KTZvqzkyImuSM3k8cDgcpKSm6w/AJK+R2+jR88okq7OvXQ0GBOlu/7jro1QsiI6FtW/Vo104V\n+aCfG5vn8qushJMn4dgx2LsXtm9Xj23b1L+7d8PVV0N6unr07auOEeisMH6eMD0/OaMXDdrZs6qo\n/+1vsHQpREdD//5w992Qk6OKe31aMcHB6oy9aVNo3Vq9QVR3+jRs2ADvvQcTJ6o3g4EDVdG/5RZo\n2dKr6QnhMjmjF8bZuhUWL4a//10V5bvvhrvugk6d/BvH/v2Qnw8rVsAHH8Add8CECdCzp3/jEGZz\npXZKoRfG+PRTePxx2LQJ7rlHFffExMD4ArWiAl5+Gf7yF4iJUQV/2DD1Ba8QnpALpgKIyXN5dee2\nYYNqj9x2G6SlwX/+A08+CXa7d4q8N/KLiIBHHoHiYlXkc3Kgc2dYsEC1mHTSPX6+Znp+rpBCLyxr\n40a4+WbVEhk6FHbtgvHj1VTIQBUaquJduxbefBNeeQWuvRY++kh3ZMJk0roRlnP8ODz6KLz2GsyY\nAf/932raoxU5nfD66zB9upoB9OSTcNVVuqMSViKtG2GcDz9UX2aWlMCWLTBmjHWLPKjW0siRanpm\n9+5qJs/cufrbOcIsUuj9xOQ+oT9yO3xYtWVGjoSnn1ZnwW3a+PywgH/ya9pUfUrZtAneeUd911BR\n4fPDAmb/bYL5+bnC40I/evRoIiIiSEhIqPV1h8NBixYtSEpKIikpiSeeeMLTQ4oG5quvIClJXZ36\n5ZdqTrqpOnYEhwP69IFrroH339cdkTCBxz36devWERYWxr333ssXX3xR43WHw8Gzzz7L8uXL6w5E\nevSiFu++q3rwzzwDv/mN7mj8a80aGDVK5T9jhjWushX+55ceff/+/Wl5iUv+pICL+nI6VYvmd7+D\nZcsaXpEHdVXtpk1q+ugtt6ilF4Rwh8979DabjU8++YTExETS09PZunWrrw8ZkEzuE3o7t5MnYfRo\ntXTB+vVqNopOOscuIgJWrlRr7wweDIcOef8YJv9tgvn5ucLnhf6aa66hpKSEzz//nIkTJzJs2DBf\nH1JY2IkT6uz1xx/V3PKOHXVHpF9oKLz6qvqeIiXFf1/SCnN4ZR59cXExQ4cOrbVHf6FOnTqxceNG\nWrVqdX4gNhujRo0iOjoagPDwcOx2e9Wqc+felWXb3O3Tp+HFF1No0gSysx0EBwdWfLq3nU748MMU\nFi+GP/3JQdu2gRWfbPtn2+FwsHDhQgCio6N5/PHH/bPWTV2FvqKigjZt2mCz2SgsLOTOO++kuLi4\nZiDyZWyDduYMjBih/l26VNaAqcuLL8JTT8HHH0NUlO5ohG5++TJ25MiR9OvXjx07dhAVFcX8+fPJ\nzc0lNzcXgDfffJOEhATsdjtTpkzh9ddf9/SQlnTuHdlEnuZWWQn33gtHj8KSJYFX5ANt7CZOVOvl\nZGSoFpenAi0/bzM9P1d4PGErLy+vztfHjx/P+PHjPT2MMJTTCWPHqr7zihXWvsrVnx56SC2Qdscd\nagpqoL05isAia90Irf7yF3jpJTW7JixMdzTWcuaMWuq4TRv4618DYzlm4X+yHr0IaAUFkJmpes2x\nsbqjsaYjR9RMnFtuUcsgi4ZHFjULICb3Cd3JraJCtR3++tfAL/KBPHZhYarllZurrqR1RyDn5w2m\n5+cKKfTC786cgeHD1aX9Q4fqjsb62rZVNzDJyoLvv9cdjQhE0roRfvfII/DZZ+pMNDhYdzTmeOAB\n2LdP3dBE+vUNh/ToRcDZskWt4bJlizoTFd5z4oS6W9WUKWoJCdEwSI8+gJjcJ3Q1t8pK+O1vYdYs\naxV5q4xdkybw97/DtGnqtoquskp+7jI9P1dIoRd+8+KL0KwZ3Hef7kjM1aOHKvSTJ+uORAQSad0I\nvyguht691Xz5QJ9lY3WnTqmC//zz6k5VwmzSuhEBY9IkmDpVirw/NGoEzz6rvpw9fVp3NCIQSKH3\nE5P7hJfKbf16+PxzePBB/8TjbVYcu4wMiI5WVx1fihXzqw/T83OFFHrhU04n/P736sbXso6N/9hs\n8NxzMHMmfPut7miEbtKjFz61Zg2MGwdbt8o9T3UYPx6aN4c5c3RHInxF5tELrZxO6NtXzeseOVJ3\nNA1TcTH06gX/+Q+0aKE7GuEL8mVsADG5T3ix3N59F44fV8sdWJmVxy46GtLT1SqhF2Pl/Fxhen6u\nkEIvfObPf4bp0yFI/sq0mjZNjcWJE7ojEbpI60b4xPbtavncvXvlS9hAMHQoDBkC2dm6IxHeJq0b\noc28eWq5AynygWHaNDW3Xs6lGiYp9H5icp/wwtwOH4bFi805ezRh7K6/Xk25LCio+ZoJ+dXF9Pxc\n4XGhHz16NBERESQkJFx0n0mTJhEbG0tiYiJFRUWeHlIEuMWLYcAAiIrSHYk4x2ZT69UvWKA7EqGD\nxz36devWERYWxr333ssXX3xR4/WVK1eSk5PDypUr2bBhA5MnT6agltMK6dGbo29fmDEDbr5ZdySi\nurIySEiA0lJo2lR3NMJb/NKj79+/Py1btrzo68uXL2fUqFEAJCcnc+jQISoqKjw9rAhQxcVqidyb\nbtIdibhQhw5qvfq339YdifA3n/foy8rKiKr2GT4yMpLS0lJfHzbgmNwnrJ7bm2/CrbdCaKi+eLzN\npLHLyoJXXz3/OZPyq43p+bnCLxelX/ixwnaR+5xlZWURHR0NQHh4OHa7nZSUFOCXwbLq9ubNmwMq\nHl9tv/FGCjNnBk48sn3+9pAhKdx/P6xY4SAsTH88sl3/bYfDwcKFCwGq6uWleGUefXFxMUOHDq21\nRz9mzBhSUlIYMWIEAHFxcaxdu5aIiIjzA5EeveXt2aNaA/v3y7o2gSwtTd1q8I47dEcivCEg5tFn\nZmayaNEiAAoKCggPD69R5IUZ3n4bhg2TIh/ohg6Fd97RHYXwJ48L/ciRI+nXrx87duwgKiqK+fPn\nk5ubS25uLgDp6el07tyZmJgYsrOzmTdvnsdBW9G5j14mOpfbmjUweLDeWHzBtLEbMgTee0/dwxfM\ny+9CpufnCo/PvfLy8i65T05OjqeHEQHu1ClYtw5+/vAmAljHjmoGTkEB/Nd/6Y5G+IOsdSO8Yt06\ndeu6zz7THYlwxbRpap36Rx7RHYnwVED06EXDsGYNDByoOwrhquuvh48+0h2F8Bcp9H5icp/Q4XDw\nwQfmXiRl4tj166daN2fOmJlfdabn5wop9MJjlZVQVKSmVgpraN1a9elrmREtDCQ9euGxLVvgzjvV\nGvTCOn73O+jRAyZN0h2J8IT06IVffPop9OmjOwpRX337qrET5pNC7ycm9wmXLXMYXehNHbuePdWn\nMVPzO8f0/FwhhV54bMcO6NVLdxSivrp3h6+/htOndUcifE169MIjZ8/C5Zertc5btNAdjaiv+HhY\nskSd3Qtrkh698LnSUlXopchbU2Kiat8Is0mh9xNT+4Tbt0Pbtg7dYfiUqWMH0K0b5Oc7dIfhUyaP\nn6uk0AuPbNum1k4R1tSlC5SX645C+Jr06IVHxo+HuDiYOFF3JMIdBQVq7GSapXVJj174XHExuHiT\nGxGAunSB3bt1RyF8TQq9n5jaJ9y3D775xqE7DJ8ydexALYVw8qSDH37QHYnvmDx+rpJCLzyybx/I\nDcOsy2aDNm2gpER3JMKXpEcv3PbjjxAVpf69yP3ehQUMHqzuJXDzzbojEe6QHr3wqdJSiIyUIm91\n7drJzBvTSaH3ExP7hN98o9o2JuZWnen5nT7tMLrQmz5+rvC40Ofn5xMXF0dsbCxz586t8brD4aBF\nixYkJSWRlJTEE0884ekhRYCoqIArr9QdhfDUFVfA/v26oxC+5NHNwSsrK5kwYQJr1qyhQ4cO9OnT\nh8zMTOLj48/b78Ybb2T58uUeBWp1KSkpukPwum++UV/kmZhbdabnd911KSxbpjsK3zF9/Fzh0Rl9\nYWEhMTExREdHExoayogRI1hWy1+MfMlqpm++kTN6E1xxBXz/ve4ohC95VOjLysqIioqq2o6MjKSs\nrOy8fWw2G5988gmJiYmkp6ezdetWTw5pWSb2CX/4QRUJE3OrzvT89uxxGF3oTR8/V3jUurG5MN3i\nmmuuoaSkhKZNm/Lee+8xbNgwvv7661r3zcrKIvrnyyzDw8Ox2+1VH7vODZZVtzdv3hxQ8Xhje8cO\n6Ns3cOKRbfe21TLTDhyOwIhHtuvedjgcLFy4EKCqXl6KR/PoCwoKmDFjBvn5+QDMnj2boKAgpk+f\nftHf6dSpExs3bqRVq1bnByLz6C1nyBDIzoahQ3VHIjzxww/QqRMcOqQ7EuEOn8+j7927Nzt37qS4\nuJhTp06xZMkSMjMzz9unoqKiKojCwkKcTmeNIi+s6dAhCA/XHYXwVPPmcOQIyHmWuTwq9CEhIeTk\n5DB48GC6devG8OHDiY+PJzc3l9zcXADefPNNEhISsNvtTJkyhddff90rgVvNuY9eJvnpJ1UkTMyt\nOtPz++gjB6GhcPy47kh8w/Txc4VHPXqAtLQ00tLSznsuOzu76ufx48czfvx4Tw8jAtDRo9CsmXzk\nN0Hz5nD4MDRtqjsS4Quy1o1wW7t2sHEjtG+vOxLhqc6dYfVqtWyxsBZZ60b41LkzemF9l11mbutG\nSKH3GxP7hMeOqQJhYm7VNYT8mjSBEyd0R+Ibpo+fK6TQC7dUVsLZsxAaqjsS4Q1yRm826dELtxw7\npq6KleJghtRUeOghGDRIdySivqRHL3zm5Elo3Fh3FMJbQkPh9GndUQhfkULvJ6b1CU+dgkaN1M+m\n5XahhpBfo0bmFnrTx88VUuiFW86ckf68SUJD1Zu3MJP06IVbiovhxhth717dkQhvuOsuyMiAu+/W\nHYmoL+nRC585cwZCPL6uWgSKkBA1k0qYSQq9n5jWJ6xe6E3L7UINIb/gYDWmJjJ9/FwhhV645exZ\nCA7WHYXwluBgOaM3mfTohVu++AJGjoQvv9QdifCG7GxISoIxY3RHIupLevTCZyor5YzeJEFB6lOa\nMJMUej8xrU949qwqDmBebhdqCPkFB5tb6E0fP1dIoRducTp/KfTC+oKCpEdvMunRC7d8+imMHQuf\nfaY7EuENkyerNeknT9Ydiagv6dELn3E6wWbTHYXwJjnPMpcUej8xrU9YvXVjWm4Xagj52WzmFnrT\nx88VHhf6/Px84uLiiI2NZe7cubXuM2nSJGJjY0lMTKSoqMjTQ4oAYGpRaKhMLvTCw0JfWVnJhAkT\nyM/PZ+vWreTl5bFt27bz9lm5ciW7du1i586dvPzyy4wdO9ajgK0qJSVFdwhed651Y2Ju1TWE/Exu\nw5k+fq7wqNAXFhYSExNDdHQ0oaGhjBgxgmXLlp23z/Llyxk1ahQAycnJHDp0iIqKCk8OK4TwATmj\nN5dHhb6srIyoqKiq7cjISMrKyi65T2lpqSeHtSST+4Qm5wYNI7+tW2H3bt2R+Ibp4+cKj9YftLn4\nee/CqT8X+72srCyio6MBCA8Px263V33sOjdYVt3evHlzQMXj6famTQ5++gkgMOKRbc+2P/rIgTpH\nC4x4ZPvi2w6Hg4ULFwJU1ctL8WgefUFBATNmzCA/Px+A2bNnExQUxPTp06v2GTNmDCkpKYwYMQKA\nuLg41q5dS0RExPmByDx6SykogClT1L/C+h58EDp0gKlTdUci6svn8+h79+7Nzp07KS4u5tSpUyxZ\nsoTMzMzz9snMzGTRokWAemMIDw+vUeSFEPqZ/IVsQ+dRoQ8JCSEnJ4fBgwfTrVs3hg8fTnx8PLm5\nueTm5gKQnp5O586diYmJITs7m3nz5nklcKs599HLJOdOIkzMrbqGkJ/JH6ZNHz9XeHyPoLS0NNLS\n0s57Ljs7+7ztnJwcTw8jAoyc/ZlHxtRcstaNcMuGDTBxIhQW6o5EeMOUKXDVVfDAA7ojEfUla90I\nnwkKknnXJpG1i8wmhd5PTOsT2my/rF9uWm4Xagj5mVzoTR8/V0ihF26RM3qzVL+RjDCP9OiFWzZv\nhlGj4PPPdUcivGHcOOjRQ/0rrEV69MJnguQeo0aRM3qzydD6iWl9wur3GDUttws1hPwqK6VHbzIp\n9MItQXKPUaNUVkKIx1fViEAlPXrhlq+/hvR02LVLdyTCG0aNggEDICtLdySivqRHL3wmNFTO6E1y\n5oyc0ZtMCr2fmNYnDAmB06fVz6bldqGGkJ/Jhd708XOFFHrhlpAQdRYozHD6NDRqpDsK4SvSoxdu\n+f57iI2Fgwd1RyK8YcgQyM6GoUN1RyLqS3r0wmcaNYKTJ3VHIbxFzujNJoXeT0zrEzZu/EuhNy23\nCzWE/E6eNLfQmz5+rpBCL9xybtaNzLwxw/HjcNlluqMQviI9euG2pk3h22+hWTPdkQhP9ewJr70G\niYm6IxH1JT164VNNm8LRo7qjEN4gZ/Rmk0LvJyb2CZs1g2PHzMytuoaQ37Fj5hZ608fPFW5fInHw\n4EGGDx/O3r17iY6O5o033iA8PLzGftHR0Vx++eUEBwcTGhpKodx7zhjNmskZvSkOH4bLL9cdhfAV\nt3v006ZNo3Xr1kybNo25c+fyww8/MGfOnBr7derUiY0bN9KqVau6A5EeveVcey28+CIkJ+uORHjC\n6VQXwJ06pVYlFdbi0x798uXLGTVqFACjRo3i7bffvui+UsDN1KIF/Pij7iiEp44ehSZNpMibzO1C\nX1FRQUREBAARERFUVFTUup/NZmPgwIH07t2bV155xd3DWZ6JfcLwcDh0yMzcqjM9v3/9y0Hz5rqj\n8B3Tx88VdfboU1NTOXDgQI3nZ86ced62zWbDdpG7Fnz88ce0a9eOb7/9ltTUVOLi4ujfv3+t+2Zl\nZREdHQ1AeHg4drudlJQU4JfBsur25s2bAyoeb2wfOwaHDqXQpk1gxCPb7m3/9BM0aeLA4QiMeGS7\n7m2Hw8HChQsBqurlpbjdo4+Li8PhcNC2bVv279/PgAED2L59e52/8/jjjxMWFsbUqVNrBiI9est5\n6CFo3RqmT9cdifDE2rXwxz/CunW6IxHu8GmPPjMzk1dffRWAV199lWHDhtXY59ixYxw+fBiAo0eP\nsmrVKhISEtw9pAgwrVvDd9/pjkJ46uBBuOIK3VEIX3K70D/88MOsXr2arl278sEHH/Dwww8DUF5e\nTkZGBgAHDhygf//+2O12kpOTGTJkCIMGDfJO5BZz7qOXSdq0gW++MTO36kzP75NPHEYXetPHzxVu\nz6Nv1aoVa9asqfF8+/bteffddwHo3LlzVW9amCciQhV6YW0HD6qxFOaStW6E2z77TK1hvnGj7kiE\nJ8aPh/h4mDBBdyTCHbLWjfCpdu2gvFx3FMJT5eXQvr3uKIQvSaH3ExP7hG3bqo/9q1Y5dIfiUyaO\nXXXbtzuMLvSmj58rpNALtwUHqzNBmXljbd99J2f0ppMevfDIjTfCjBkwYIDuSIQ7jh+Hli3VMgiy\nBII1SY9e+FzHjrB3r+4ohLv27IGrrpIibzop9H5iap8wJgbef9+hOwyfMnXsAHbvhvBwh+4wfMrk\n8XOVFHrhkfh4OaO3st27pT/fEEiPXnhkyxYYORK++kp3JMId48ZBXBxMmqQ7EuEu6dELn4uNhf/8\nB86c0R2JcMeWLerG4MJsUuj9xNQ+4WWXQcuWDnbt0h2J75g6dmfPqkL/448O3aH4lKnjVx9S6IXH\nYmKgqEh3FKK+9u5V94lt0UJ3JMLXpEcvPDZnjlrc7NlndUci6mPZMnj5Zfh5DUJhUdKjF37Rpw98\n+qnuKER9bdwIdrvuKIQ/SKH3E5P7hMePOygqMvcLWVPH7qOP4Prrzc3vHNPzc4UUeuGxsDA1F3vr\nVt2RCFedPq0+hV13ne5IhD9Ij154xW9/C4mJMHGi7kiEKwoL4f774fPPdUciPCU9euE3AwdCLTcc\nEwHqXNtGNAxS6P3E5D6hw+HgV7+CtWvN7NObOHb//jfccIP62cT8qjM9P1e4XeiXLl1K9+7dCQ4O\nZtOmTRfdLz8/n7i4OGJjY5k7d667hxMBrk0biI6W2TdWcPy4elMeNEh3JMJf3O7Rb9++naCgILKz\ns3nmmWe45pprauxTWVnJ1VdfzZo1a+jQoQN9+vQhLy+P+Pj4moFIj97ypk5VF988+qjuSERd3n0X\nnnoK5ETXDD7t0cfFxdG1a9c69yksLCQmJobo6GhCQ0MZMWIEy5Ytc/eQIsClp8Py5bqjEJfyzjsw\nZIjuKIQ/+bRHX1ZWRlRUVNV2ZGQkZWVlvjxkwDK5T3gutxtvhH371NK3JjFp7JxOWLEChg795TmT\n8quN6fm5IqSuF1NTUzlw4ECN52fNmsXQ6n8pF2Gz2eoVTFZWFtHR0QCEh4djt9tJSUkBfhksq25v\n3rw5oOLx1fbtt6ewdCn07RsY8cj2+dthYSlcdhmUlzvYv19/PLJd/22Hw8HChQsBqurlpXg8j37A\ngAEX7dEXFBQwY8YM8vPzAZg9ezZBQUFMnz69ZiDSozfCv/+tevV1fD8vNJo8GcLD4fHHdUcivMVv\n8+gvdpDevXuzc+dOiouLOXXqFEuWLCEzM9MbhxQB6oYboLwcdu7UHYm40KlTkJcHo0bpjkT4m9uF\n/q233iIqKoqCggIyMjJIS0sDoLy8nIyMDABCQkLIyclh8ODBdOvWjeHDh9c646YhOPfRy0TVcwsO\nhjvvhEWL9MXjbaaM3YoV0K0bdO58/vOm5Hcxpufnijp79HW59dZbufXWW2s83759e96ttu5pWlpa\n1ZuAaBjGjIGbboJHHoFGjXRHI85ZuBCysnRHIXSQtW6ET9x0k1r/ZuRI3ZEIgAMH1L1hS0vVInTC\nHLLWjdBm/HjIydEdhTgnJwdGjJAi31BJofcTk/uEteWWmanm1P88q9TSrD52hw9Dbi489FDtr1s9\nv0sxPT9XSKEXPhESAmPHwnPP6Y5EvPyyaqV16aI7EqGL9OiFzxw6BLGxsG6d6g8L/zt5Us2yWbEC\nkpJ0RyN8QXr0QqvwcHjwQXjsMd2RNFyLF0NCghT5hk4KvZ+Y3CesK7dJk+DDD63dq7fq2B0/Dn/6\nE/zhD3XvZ9X8XGV6fq6QQi98qlkz+J//UXPqhX89/TT06QP9++uOROgmPXrhcydPql793/4mRcdf\nSkrAboeNG9UNYYS5pEcvAkLjxurscswYVfSF7z38MIwbJ0VeKFLo/cTkPqErud1xh5reZ8W7SVpt\n7D75RN0qsJZFYmtltfzqy/T8XCGFXviFzQbz5sGLL8LWrbqjMdeJE5CdDU8+KVfBil9Ij1741bx5\nqle/bh0EyWmG102ZAmVl8MYb6s1VmE969CLgjBmj/pV1cLzvX/+Cf/xDLXcgRV5UJ4XeT0zuE9Yn\nt6AgePVVeOIJKCjwXUzeZIWx+/ZbGD1a3QegVav6/a4V8vOE6fm5Qgq98LuYGPjrX9UXtBUVuqOx\nPqcT7r8f7r4bBgzQHY0IRNKjF9o8+qi6anbNGrUImnDPn/4E77yjvvdo3Fh3NMLfpEcvAtpjj8Fl\nl7k+DVDUtGgRzJ8Py5dLkRcX53ahX7p0Kd27dyc4OJhNmzZddL/o6Gh69uxJUlIS1157rbuHszyT\n+4Tu5hYcrGbgvP02vPKKd2PypkAdu/ffV2vMr1wJbdu6/98J1Py8xfT8XOH2B+aEhATeeustsrOz\n69zPZrPhcDhoVd9viESD0KoV5Oer3nKTJvCb3+iOyBq+/FLdpnHpUoiP1x2NCHRuF/q4eiwwLr13\nSElJ0R2Cz3iaW2wsrFqlbo7RuDHcead34vKWQBu7PXsgIwP+/Ge48UbP/3uBlp+3mZ6fK3zeo7fZ\nbAwcOJDevXvzSiB/Phdadeum5oFPmqRaOaJ2X34JN9ygvte46y7d0QirqLPQp6amkpCQUOPxzjvv\nuHyAjz9Znu8oAAAJQ0lEQVT+mKKiIt577z1eeukl1q1b53HQVmRyn9BbufXsCe++qy7hX77cK/9J\nrwiUsSsogIED1fIG48Z5778bKPn5iun5uaLO1s3q1as9PkC7du0AuPLKK7n11lspLCyk/0XWqs3K\nyiL65+X2wsPDsdvtVR+7zg2WVbc3/3znjUCJJ5C3V6yAtDQHd9wB8+alYLMFVnw6tp9+2sETT8Df\n/55Cerr+eGRb37bD4WDhwoUAVfXyUjyeRz9gwACefvppevXqVeO1Y8eOUVlZSfPmzTl69CiDBg3i\nscceY9CgQTUDkXn0opp9++CWWyAxUV3S35CnDublqTVs/vEPuP563dGIQOPTefRvvfUWUVFRFBQU\nkJGRQVpaGgDl5eVkZGQAcODAAfr374/dbic5OZkhQ4bUWuSFuFDHjvDRR3DkCPzqVw3zCtrjx2Hs\nWHV3rtWrpcgL98mVsX7icDiqPoaZxpe5nT0LM2ao9XHy8qBfP58cpk46xm7bNhg+HHr0gL/8BS6/\n3HfHMvlvE8zPT66MFZYXFAT/+79qKuHtt8Pkyeos31ROp7rS9YYbVK5/+5tvi7xoGOSMXljG99/D\nAw+ols7LL6sZKCbZswemToWvv4YlS6B7d90RCSuQM3phlCuuUGu75OSoJXl/+1v44QfdUXnu8GH4\n/e+hd2+45hr49FMp8sK7pND7ybnpUSbyd27p6erCocaN1VW1jz4KBw/67ni+yu/sWdWmufpqdVeo\nLVvgj39UC735k8l/m2B+fq6QQi8s6fLL4aWXoLAQystVwf/jH1V7J9CdPAmLF0OvXmpd/mXL1JfN\nHTrojkyYSnr0wgh79sDs2Wqu+f33w333qeIfSEpL1TUBr7yirgKeOBGGDJHb/gnPSI9eNBidOqkv\naDdtghMnoH9/6NNHzdbZv19fXCdPqgXb7rhDFfdDh8DhUM8NHSpFXviHFHo/MblPGEi5XXWVKu6l\npTBzJmzerBZMS01VZ9LbtqneeH3UN7/iYvi//4PMTLjySnUdwA03qOdffBHqsfCrXwTS+PmC6fm5\nQm7gJowUEgKDBqnH8eOwYgW89ZZq7/zwAyQnw3XXqUevXmpd/PqeXZ89q95Qtm9Xj23bYO1a+O47\nuPlmtV78ggVqtpAQOkmPXjQ4Bw7Ahg1qNcj166GoCE6dUndpqv648kqorFStoOqP48fV2fmOHdCi\nhbrxR1ycevTtq944guSzsvATV2qnFHohgKNH1Xo6Bw788vj2W/XJoEkT9Wjc+JefO3ZU0yJbtNAd\nuWjopNAHEJPX2zA5N5D8rM70/GTWjRBCCDmjF0IIK5MzeiGEEFLo/cXkubwm5waSn9WZnp8rpNAL\nIYThpEcvhBAWJj16IYQQ7hf6hx56iPj4eBITE7ntttv48ccfa90vPz+fuLg4YmNjmTt3rtuBWp3J\nfUKTcwPJz+pMz88Vbhf6QYMG8dVXX/H555/TtWtXZs+eXWOfyspKJkyYQH5+Plu3biUvL49t27Z5\nFLBVbd68WXcIPmNybiD5WZ3p+bnC7UKfmppK0M8LeiQnJ1NaWlpjn8LCQmJiYoiOjiY0NJQRI0aw\nbNky96O1sEOHDukOwWdMzg0kP6szPT9XeKVHP3/+fNLT02s8X1ZWRlRUVNV2ZGQkZWVl3jikEEII\nF9W5THFqaioHDhyo8fysWbMYOnQoADNnzqRRo0bcddddNfazyV0VqhQXF+sOwWdMzg0kP6szPT+X\nOD2wYMECZ79+/ZzHjx+v9fX169c7Bw8eXLU9a9Ys55w5c2rdt0uXLk5AHvKQhzzkUY9Hly5dLlmr\n3Z5Hn5+fz9SpU1m7di2tW7eudZ8zZ85w9dVX8/7779O+fXuuvfZa8vLyiI+Pd+eQQggh3OB2j37i\nxIkcOXKE1NRUkpKSGDduHADl5eVkZGQAEBISQk5ODoMHD6Zbt24MHz5cirwQQvhZwFwZK4QQwjcC\n7srYZ555hqCgIA4ePKg7FK965JFHSExMxG63c9NNN1FSUqI7JK9y9QI6q1q6dCndu3cnODiYTZs2\n6Q7HK0y/mHH06NFERESQkJCgOxSvKykpYcCAAXTv3p0ePXrwwgsv1P0L9f8K1nf27dvnHDx4sDM6\nOtr5/fff6w7Hq3766aeqn1944QXnfffdpzEa71u1apWzsrLS6XQ6ndOnT3dOnz5dc0TetW3bNueO\nHTucKSkpzo0bN+oOx2NnzpxxdunSxblnzx7nqVOnnImJic6tW7fqDsurPvzwQ+emTZucPXr00B2K\n1+3fv99ZVFTkdDqdzsOHDzu7du1a5/gF1Bn9gw8+yJNPPqk7DJ9o3rx51c9Hjhy56BfYVuXKBXRW\nFhcXR9euXXWH4TUN4WLG/v3707JlS91h+ETbtm2x2+0AhIWFER8fT3l5+UX3r3MevT8tW7aMyMhI\nevbsqTsUn/nDH/7Aa6+9RtOmTSkoKNAdjs/Mnz+fkSNH6g5D1KG2ixk3bNigMSLhruLiYoqKikhO\nTr7oPn4t9Be7AGvmzJnMnj2bVatWVT3ntOB3xJe6wGzmzJnMnDmTOXPm8MADD7BgwQINUbrP0wvo\nAp0r+ZlCLmY0w5EjR/j1r3/N888/T1hY2EX382uhX716da3Pf/nll+zZs4fExEQASktL6dWrF4WF\nhbRp08afIXrkYvld6K677qp1yYhAd6n8Fi5cyMqVK3n//ff9FJF3uTp+JujQocN5EwJKSkqIjIzU\nGJGor9OnT3P77bdzzz33MGzYsDr3DYjWTY8ePaioqKja7tSpExs3bqRVq1Yao/KunTt3EhsbC6g2\nVVJSkuaIvCs/P5+nnnqKtWvX0qRJE93h+JQVP21eqHfv3uzcuZPi4mLat2/PkiVLyMvL0x2WcJHT\n6eS+++6jW7duTJkyxaVfCDidOnUybtbN7bff7uzRo4czMTHRedtttzkrKip0h+RVMTExzo4dOzrt\ndrvTbrc7x44dqzskr/rnP//pjIyMdDZp0sQZERHhvPnmm3WH5LGVK1c6u3bt6uzSpYtz1qxZusPx\nuhEjRjjbtWvnbNSokTMyMtI5f/583SF5zbp165w2m82ZmJhY9f/ce++9d9H95YIpIYQwXEBNrxRC\nCOF9UuiFEMJwUuiFEMJwUuiFEMJwUuiFEMJwUuiFEMJwUuiFEMJwUuiFEMJw/w9GDmyVx/PapwAA\nAABJRU5ErkJggg==\n", "text": [ "" ] } ], "prompt_number": 28 }, { "cell_type": "code", "collapsed": false, "input": [ "def rk4_step(y, t, h, f):\n", " k1 = f(t, y)\n", " k2 = f(t+h/2, y + h/2*k1)\n", " k3 = f(t+h/2, y + h/2*k2)\n", " k4 = f(t+h, y + h*k3)\n", " return y + h/6*(k1 + 2*k2 + 2*k3 + k4)\n", "\n", "plot_stability_region(-1, rk4_step)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAW0AAAEACAYAAAB4ayemAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4VdW9xvFvwiggpKABBCQqIEnAJKKmIGIQAyLiFaEX\nVITgVOuIWqgDKt4KlFYKXKEOrTYqXhx6r+WCGBnkgKCIKAGJXsABCFDGCoJMGc79YzXIcBIOOefs\ntfc+7+d58sQT4OzXxc6Pnd9ee62EYDAYREREPCHRdgAREQmfiraIiIeoaIuIeIiKtoiIh6hoi4h4\niIq2iIiHRFS0Dx48SHZ2NpmZmaSlpfHII49EK5eIiISQEOk87f3791OvXj1KS0vp2rUrzzzzDF27\ndo1WPhEROUrE7ZF69eoBcPjwYcrKymjcuHHEoUREJLSIi3Z5eTmZmZk0bdqU7t27k5aWFo1cIiIS\nQsRFOzExkcLCQjZt2sSiRYsIBAJRiCUiIqHUjNYbNWrUiD59+rB8+XJycnKOfL1NmzZ888030TqM\niEhcyMjIoLCw8ISvR1S0d+7cSc2aNUlKSuLAgQPMnTuXJ5988pjf88033xDva1KNHj2a0aNH247h\nen4ep2AQPv8cVq2CHTtg+/afPlf8dzAIl10G3bubj/POg4SE0O/n57GKJi+PU0Ilf/kRFe1//OMf\nDB06lPLycsrLy7n55pvp0aNHJG/pS+vXr7cdwRP8OE779sH06fD88/D996YoJydD06bQoYP57zPP\nNJ9LSmDhQliwAJ58EmrUgCuuMAX8qqugWbOf3tePYxULfhyniIp2x44d+fzzz6OVRcQ3vvgCXngB\n/uu/oFs3GDMGevaExJPcRWrTBm691Vx1r11rCvh778GIEeY9br+98qtviQ8Rz9M+6QESEuK+PRII\nBI7p80toXh+nYBDeeAOmToXvvoPbbjMfrVpF/t5FRTBkiLlC/8tfYO1ab4+VU7x8TlVWO1W0RaJg\n+3YYNgy2bYNHH4W+faFWregeo6QEnn7atFomT4ZBg6L7/uIuldVOrT3iAE2DDI9Xx2nePMjKggsu\ngI8/huuvj37BBvOeTz0Fs2bByJEBBg2CXbuifxw/8eo5VRUVbZFqKimBhx+GoUPhlVdg3LjYFOvj\nXXwx/PnP5sZk587mZqfED7VHRKrh22/hhhvgjDMgP9/MALEhLw/q1jUtE/EXtUdEouSNN+DnPzdF\ne9YsewUbTG+7oABmz7aXQZylou0AP/bVYsEL4zRuHDz2mCmUw4fbm35XMVaNGpkr/dtvh5077WRx\nMy+cU6dKRVskDMGgKdbTpsHixXDhhbYT/SQnx8wk+dWvTE7xN/W0RU4iGDRX1R9+CHPmmD622xw8\nCBddZG6MDh5sO41Eg+Zpi1RDWRn88pfw1Vfw7ruQlGQ7UeVWrIBeveCzz6LzQI/YpRuRFvmxrxYL\nbhunkhJz1frdd/D+++4q2KHGKivLPDX53HPO53Ert51T0aCiLRLCoUPwi1/A3r3mCrtBA9uJwtOr\nl+m5i3+pPSJynB9/hH79zJX1tGlQu7btROHbuxeaNzdPStapYzuNRELtEZEw7NljlkFt0cIsqeql\ngg1w+unQrp1Zu1v8SUXbAX7sq8WC7XHatQt69ICMDHjpJbOetVtVNVZdu6pFUsH2ORULKtoiwNat\nZr7zlVfCs8+efN1rN1PR9jf1tCXubdxoivWQIeYBGq9vMrB5s/lpYccO7/+/xDP1tEVCKCw0W4D9\n6lcwapQ/ilyLFqa3vWaN7SQSCyraDvBjXy0WnB6nt96C3FwYPx4eeMDRQ0fsZGOlFonhx+89FW2J\nO2VlZneZkSPNY+l+3AGma1dYssR2CokF9bQlruzZAzfdZOYz/+1vdpdVjaU5c8xPEPPn204i1aWe\ntsS9wkLIzoaUFLNFmF8LNkBxsdYf8SsVbQf4sa8WC7Eap5IS+O1voWdPMztkyhRntgWLpZON1YYN\n0Lq1M1nczI/fezVtBxCJpS+/NHs4Nm5snhJs2dJ2Imds3AjdutlOIbGgnrb4Ummp2Yrrd7+Dp5+G\nO+7wx3S+cHXvbn6quPJK20mkuiqrnbrSFt9ZsgTuvttcXX/yCZx7ru1EzlN7xL/U03aAH/tqsRDp\nOG3fDsOGwcCBZgeX+fP9W7CrGquyMvNUpG5E+vN7T0VbPK+01NxcTE+HJk3MLjODBsVXO+RoW7ea\nnzLq1rWdRGIhop52cXExQ4YMYfv27SQkJHDHHXdw3333HXsA9bQlhubMMU8zNm8OkyZBhw62E9n3\n0UdmTD75xHYSiURMetq1atVi4sSJZGZmsm/fPjp16kRubi6pqamRvK3ISa1ZAw89ZD5PmAB9+8bv\nlfXx1M/2t4jaI82aNSMzMxOABg0akJqaypYtW6ISzE/82FeLhXDG6fvvzVVk165mhkRREVx7bfwV\n7KrGauNGFe0Kfvzei1pPe/369axYsYLs7OxovaXIEWVl8MILkJoK+/ebYv3QQ97bWcYJutL2t6jM\n0963bx85OTmMGjWK66677tgDqKctEVq4EO67D372M9O3/tcPd1KJPn3gzjtNy0i8K2bztEtKSujf\nvz+DBw8+oWBXyMvLIyUlBYCkpCQyMzPJyckBfvrxRa/1+vjX27bB4MEBVq2CqVNz6N8fFi4MEAi4\nI59bX3/5JZx9tnvy6HV4rwOBAPn5+QBH6mUoEV1pB4NBhg4dSpMmTZg4cWLoA+hKm0AgcOQvSSpX\nMU5lZfDii/DEE2be9RNPQIMGttO5S2XnVDAIDRuaBaOSkpzP5TZe/t6LyZX2kiVLmDZtGhdccAFZ\nWVkAjBs3jquuuiqSt5U4tno13HqrWdDpgw+gY0fbibxlzRozR1sF27+09oi4QlkZPPMM/OEPMHYs\n3HabtzfXteW558z87H/9lC0eprVHxLXWrjUr8Z12Gixfbta7luoJBODqq22nkFjStYwDKm42yLGC\nQdO77tLF7CYzalRABTtMoc6pYJB/3aR1Oo17+fF7T1faYsWPP5od0FesMKvynX++KThSfV9+CfXr\na46236mnLY5bswYGDICsLNODrV/fdiJ/mDLF/CP40ku2k0g0aI9IcYWZM+Gyy+Dee+GVV1Swo2nB\nAvNov/ibirYD/NhXq47nn4df/hJmzQq9k4zGKXzHj1V5uXlyVEX7WH48p9TTlpgLBmHUKHj7bfjw\nQzjvPNuJ/OeLL8z87BYtbCeRWFNPW2KqpMTMuV67Fv73f+HMM20n8qfJk82NyBdesJ1EokU9bXFc\nWRkMGQI7dpitv1SwY0f97Pihou0AP/bVTqa83DyOvnMn/M//QL16J/8z8ThO1XX0WJWVwaJFmp8d\nih/PKfW0JeqCQTMH+7vvYPZs7VUYaytXQrNm5kP8Tz1tibqRI80Nxzlz4PTTbafxvwkT4Jtv4E9/\nsp1Eokk9bXHEK6+YdsisWSrYTvnoI+jWzXYKcYqKtgP82FcLZelSGDHCzBJp0uTU/3y8jFM0HD1W\nq1drCdvK+PGcUtGWqNi0Cfr3h5dfhrQ022nix8GDZiPftm1tJxGnqKctESstNY+m/9u/wcMP204T\nXwoLzQqJRUW2k0i0qactMTN+vFlDZORI20niT1ERpKfbTiFOUtF2gB/7ahVWrDBP4/31r5HvNOPn\ncYq2irFS0a6aH88pFW2ptoMHYfBg+OMfoVUr22nik4p2/FFPW6pt1Cj4v/8zC0Edv2KfOKNNG7Pc\nbWqq7SQSbZXVThVtqZb166FTJ/M0XsuWttPEp/37zdTKH34wu9eLv+hGpEV+7Kv95jdw//3RLdh+\nHKdYCQQCfPMNnHOOCnZV/HhOae0ROWWLF8PHH5ubj2LP1q3QvLntFOI0tUfklASDkJ0Nw4fDjTfa\nThPfpk2D996D11+3nURiQe0RiYq5c00vddAg20lk61Zo2tR2CnGairYD/NRXGz/e9LMjnZMdip/G\nKdYCgQDbtqlon4wfzykVbQnbsmXw9de6ynaLrVu1hnY8Uk9bwta/P1x+Odx3n+0kAtCzJzz4IFx1\nle0kEguV1U7NHpGwbN5s9iF89VXbSaTC99/Dz35mO4U4LeL2yC233ELTpk3pqAV9K+WHvtq0aTBg\ngFkYKlb8ME5OCQQC7N0LDRvaTuJufjynIi7aw4YNo6CgIBpZxKWCQcjPh7w820nkaHv3anegeBSV\nnvb69evp27cvX3zxxYkHUE/b8z75BG6+Gdas0RojbtKwIRQXQ6NGtpNILGietlTba6/BkCEq2G4S\nDMKPP0KDBraTiNMcuRGZl5dHSkoKAElJSWRmZpKTkwP81HPy8+vCwkKGDx/umjyn8nrBggBvvQUL\nF8b+eEf3H93y/+/W1wcPQu3aOXz4oTvyuPX1pEmTPFNvAoEA+fn5AEfqZShqjzggEAgc+UvymlWr\noF8/Mz871lfaXh4np82cGeDmm3PYvdt2Enfz8jml9ohFXj1pAGbNgmuucaY14uVxctrFF+dQp47t\nFO7nx3Mq4qJ9ww030KVLF9auXUurVq34q5Z+85VZs6BPH9sp5HiHDqGiHaciLtrTp09ny5YtHDp0\niOLiYoYNGxaNXL5ydK/WS/buNe2Ryy935nheHScbFi0KqGiHwY/nlNojUqlPPoGsLF3RuVFpqTY/\niFdae0Qq9R//YaaVjR9vO4kcb+VKM3d+1SrbSSRWdCNSTtmSJdCli+0UEkppKdTUykFxSUXbAV7s\nq5WXw9Kl0Lmzc8f04jjZsmxZgBo1bKdwPz+eUyraEtLGjWZdi+Rk20kklPJyXWnHKxVtB3hxrujq\n1ZCe7uwxvThOtmRm5uhKOwx+PKdUtCWkoiLo0MF2CqlMWVlstnwT99NfuwO82Fdbvdr5ou3FcbLl\n888DKtph8OM5pb92CWnNGmjf3nYKqUwwqFUX45XmaUtIyclmDrA2jnWnDz6Ap582n8WfNE9bwrZ/\nv3mEvWlT20mkMuXlutKOVyraDvBaX23DBjj7bOeLgtfGyaaVKwMq2mHw4zmloi0n2LABWre2nUKq\noo5j/FLRdoDX5opu3QrNmzt/XK+Nk00XXJCjK+0w+PGcUtGWE2zfrichvUBFOz6paDvAa301W0Xb\na+Nk08qVAdsRPMGP55SKtpxAV9rup552/FLRdoDX+mq7dkGTJs4f12vjZFNGhnra4fDjOaWiLSf4\n4Qdo1Mh2ChEJRUXbAV7rq+3ZAw0bOn9cr42TTepph8eP55SKtpxAV9oi7qW1R+QEjRvDunV2+toS\nnvffhz/+0XwWf9LaIxK2AwegXj3bKUQkFBVtB3iprxYMwqFDUKeO88f20jjZpp52ePx4TqloyzEO\nHoTatbUriohbqactx/j+ezjnHNi923YSqYp62v6nnraEpaTEXGmLiDupaDvAS3210lKoWdPOsb00\nTrappx0eP55TERftgoIC2rdvT9u2bRk/fnw0MolFNou2iJxcRD3tsrIyzj//fObNm0eLFi24+OKL\nmT59OqmpqT8dQD1tT/n2W7jySvNZ3Es9bf+LSU972bJltGnThpSUFGrVqsWgQYOYMWNGJG8plpWX\na+aIiJtF9O25efNmWrVqdeR1y5Yt2bx5c8Sh/MaPfbVY0DiFb+HCAEuX2k7hfn48pyLqXiaEuTZk\nXl4eKSkpACQlJZGZmXlkycSKQfXz68LCQlflqer1J58EOHAAwB159Dr06x07zBoxbsnj1teFhYWu\nylPV60AgQH5+PsCRehlKRD3tpUuXMnr0aAoKCgAYN24ciYmJ/OY3v/npAOppe8rXX8NVV5nP4l7q\naftfTHraF110EevWrWP9+vUcPnyYN998k2uvvTaStxQX0L+xIu4VUdGuWbMmU6ZMoVevXqSlpTFw\n4MBjZo6IUfEjkBfUrAllZXaO7aVxsk3ztMPjx3Mq4hm5vXv3pnfv3tHIIi5Qs6aZqy0i7qS1R+QY\nW7dCZqb5LO6lnrb/ae0RCUvNmmb9ERFxJxVtB3ipr1a3rllP2wYvjZNt6mmHx4/nlIq2HKNuXbNz\njTpaIu6knracoFYt+PFHLdHqZupp+5962hK2007jX09FiojbqGg7wGt9tQYNzJW207w2Tjappx0e\nP55TKtpygkaNYM8e2ylEJBQVbQdULA7jFQ0bmsWInOa1cbIpIyPHdgRP8OM5paItJ2jYUFfaIm6l\nou0Ar/XVkpLMruxO89o42bRyZUDTMsPgx3NKRVtOkJxs1msW9wpzKXvxIRVtB3itr5acDNu3O39c\nr42TTepph8eP55SKtpzAVtGWU6P2SHxS0XaA1/pqTZvCtm3OH9dr42TTqlUB2xE8wY/nlIq2nKBl\nSygutp1CTkZX2vFJa4/ICbZvh7Q02LnTdhKpzLx5MG4czJ9vO4nEitYekbCdeSbs3w/79tlOIpVJ\nSNCVdrxS0XaA1/pqCQnQujVs2ODscb02TjZpnnZ4/HhOqWhLSOedB+vW2U4hlUlM1JV2vFLRdoAX\n54qmp0NRkbPH9OI42XLhhTkq2mHw4zmloi0hdegAq1fbTiGVSUyE8nLbKcQGFW0HeLGvlp7ufNH2\n4jjZUlgYoKzMdgr38+M5paItIaWmwtdfw+HDtpNIKImJqGjHKc3TlkplZsILL0B2tu0kcrzly+HO\nO81n8SfN05ZT1qULLFliO4WEUrMmlJbaTiE2qGg7wKt9tUsvhY8+cu54Xh0nG1asCKhoh8GP55SK\ntlSq4kpb3S33qVEDSkpspxAbql203377bdLT06lRowaff/55NDP5jlfniqakmB/D16xx5nheHScb\nunbN0U3iMPjxnKp20e7YsSPvvPMO3bp1i2YecZGEBOjTB95913YSOV6dOnDokO0UYkO1i3b79u1p\n165dNLP4lpf7atdcA7NmOXMsL4+T05YvD6hoh8GP55R62lKlK66Azz6D3bttJ5Gj1aqlK+14VbOq\nX8zNzWXr1q0nfH3s2LH07ds37IPk5eWRkpICQFJSEpmZmUd6TRX/Evr9dQW35An39bJlAdLSoKAg\nh0GDYnu8nJwc6/+/Xnnds2cOBw7AggUBEhLs53Hr64qvuSVPVa8DgQD5+fkAR+plKBE/XNO9e3cm\nTJjAhRdeGPoAerjG8159Fd56y7k2iYTntNNg1y6oV892EomFmD5co6JctYp/Tb2qf38z9e8f/4jt\ncbw+Tk4KBAKcfjrs3Ws7ibv58ZyqdtF+5513aNWqFUuXLqVPnz707t07mrnERerXh3794PXXbSeR\no6loxyetPSJhWbQI7roLvvjCTAUU+zIyID8fsrJsJ5FY0NojEpHLLjNP4C1aZDuJVGjUCPbssZ1C\nnKai7QA/9NUSEuDXv4bx42N3DD+Mk1MCgQDJybB9u+0k7ubHc0pFW8I2ZAisXAmrVtlOIgBNm8K2\nbbZTiNPU05ZT8vvfm8Ktm5L2/fa3cPAgjBljO4nEgnraEhV33gnvvw9r19pOIrrSjk8q2g7wU1+t\nYUMYORJGjIj+e/tpnGItEAjQrJmK9sn48ZxS0ZZTdv/9ZtPfefNsJ4lvTZtCiFUmxOfU05Zqeecd\neOIJWLHCrLktztu6FdLTYedOzZ33I/W0Jaquuw7OOANefNF2kvjVtKn5rBZJfFHRdoAf+2oJCTB1\nKjz5JHz7bXTe04/jFCuBgFndLz0diopsp3EvP55TKtpSbWlp8MgjMHQolJXZThOfVLTjj3raEpHy\ncrNRwtVXm1kl4qwpU8xN4eeft51Eok09bYmJxESzaNEf/mBuSoqzdKUdf1S0HeDHvtrRUlLMFV+/\nfpGtheH3cYqmirGqKNr6YTY0P55TKtoSFQMHwk03wYABcPiw7TTxIzkZatTQfO14op62RE15ubna\nTk42UwE1d9gZOTkwahRceaXtJBJN6mlLzCUmwrRp8PHHMGGC7TTx46KLYOlS2ynEKSraDvBjX60y\np58Os2ebOdx/+tOp/dl4GqdIHT1WOTmwYIG1KK7mx3NKDyBL1J19Nsyfb4pJ3bpwyy22E/nbZZfB\nDTfAoUNQp47tNBJr6mlLzKxdC927mzW4b7rJdhp/u+QSeOYZ6NbNdhKJFvW0xXHt2sHcueahm6lT\nbafxt+7d1SKJFyraDvBjXy1caWmweDE8+yw8/LCZYVKZeB6nU3X8WKloh+bHc0pFW2LunHNgyRL4\n8EO4+WbTe5XouvRSWL4cDhywnURiTT1tccyBA6a3vXMnvPkmNG9uO5G/dO5s9ou84grbSSQa1NMW\n6047Dd5+G3r0gE6dwIc/uVqlFkl8UNF2gB/7atVVo4ZZgzs/HwYNgvHjf+pza5zCF2qsunfXP4TH\n8+M5paItVvTsCZ9+Cn//O1xzDWzZYjuR9116qVlpcf9+20kkllS0HZCTk2M7giu1agULF8LFF0NW\nFmzZkqPV6sIU6pyqV8+M45IlzudxKz9+71W7aI8YMYLU1FQyMjK4/vrr2bNnTzRzSZyoXRueegre\nfReefhp+8QvYscN2Ku/SI+3+V+2i3bNnT4qKili5ciXt2rVj3Lhx0czlK37sq0XbRRfBpEkBzj0X\nOnaEV17RGtFVqeyc0s3IY/nxe6/aRTs3N5fERPPHs7Oz2bRpU9RCSXyqXds88j5rlnkYJycHvvzS\ndipv6dwZCgu1prmfRWWedt++fbnhhhu48cYbTzyA5mlLNZSVmX0PR4+GW2+Fxx+H+vVtp/KGlBT4\n4AM491zbSSQS1ZqnnZubS8eOHU/4mDlz5pHfM2bMGGrXrh2yYItUV40acPfdsGoVFBdD+/bw+utq\nmYSjdWvYsMF2ComVKpdmnTt3bpV/OD8/n9mzZzN//vwqf19eXh4pKSkAJCUlkZmZeeSubkXPyc+v\nCwsLGT58uGvyuPX1sWtEm19fsybA7bfDXXflcP/9MHZsgHvvhTvvtJ/X5uuKr4X69dq1YcMGd+W1\n9XrSpEmeqTeBQID8/HyAI/UypGA1vffee8G0tLTgjh07qvx9ERzCNxYsWGA7giecbJzKyoLBl18O\nBps3Dwbz8oLBzZudyeVGVY3VY48Fg6NHO5fFzbz8vVdZ7ax2T7tt27YcPnyYxo0bA9C5c2f+FGKr\nEvW0Jdp++AHGjoW//AUeeAAefNA8Ii/Gn/9sth976SXbSSQSldVOLRglnvXtt2at7uXLzePw//7v\n2kwYYM4cMwtn3jzbSSQSWjDKoqP7kFK5Ux2nc8+Fv/3NzOn+3e/MFMHVq2MSzXWqGquzz9aNyAp+\n/N5T0RbPu/xyc7U9cKBZlvShh0wLJV6dfbaZcVPVhhPiXWqPiK9s3252yHn/fbNn4qBB8dkySU42\n0yWbNbOdRKpL7RGJC8nJ8PLLZt3usWOhXz/Yts12KudprrZ/qWg7wI99tViI5jh16WJaJmlpkJFh\net9+crKxUl/b8OP3noq2+FadOuZq++9/h8cegxtvhH/+03YqZ+hK27/U05a4sH8/PPoo/Pd/m6vu\n7GzbiWJr8mRYtw6mTLGdRKpLPW2Ja/XqwaRJMHUq9O1rtjvzs9atYeNG2ykkFlS0HeDHvlosODFO\n115r9lEcMwaGD4fS0pgfMiZONlbBoFkpMd758XtPRVviTloaLFsGX30FvXrBrl22E0XfRx+ZtbXF\nf9TTlrhVVgaPPGL63AUF0Lat7UTR07kzjBtnnhIVb9LaIyKVePFFs0/lnDmQnm47TeQOHIAzzjB7\nbdarZzuNVJduRFrkx75aLNgapzvuMAss9egBn31mJcIpq2qsPv0UOnRQwQZ/fu9VuQmCSLy46SZT\n5Hr3hnfegUsvtZ2o+hYvhq5dbaeQWFF7ROQo778PgwfDG2+YK28vuvpquP128wi/eJd62iJhWrQI\nBgwwmwj07Ws7zakpL4fGjWHNGmja1HYaiYR62hb5sa8WC24Zp27d4N13zdXqm2/aThNaZWNVVGQW\nzVLBNtxyTkWTetoiIVx8sZlNctVV5hH4YcNsJwqP+tn+p/aISBXWroXcXBgxAu65x3aaqv34I2Rl\nwcSJ0KeP7TQSKfW0Rapp/Xq48kq47TazwYJb3X232bHntddsJ5FoUE/bIj/21WLBreOUkmJuTr76\nKowaZdb1sO34sSoogJkz4dln7eRxK7eeU5FQ0RYJw1lnwcKFMHu2WWjKTYsx7dplfgr4618hKcl2\nGok1tUdETsHu3WY6YEkJTJsGrVrZzRMMmn0wzzrL9LLFP9QeEYmCpCQzq6R3b7joIvP0pE3Tp8MX\nX5gdeiQ+qGg7wI99tVjwyjglJpobkjNmwEMPwV13mUWanBQIBCgsNK2aadPgtNOcPb5XeOWcOhUq\n2iLV9POfw4oVZt/JSy4xD7Y4obTUFOrcXHPj8cILnTmuuIN62iIRCgbN9mUjR5olXm+/HWrVis2x\n1q2DIUOgfn14+WWz67r4k3raIjGSkGCemFy82Gwa3Lo1PP54dPdoLC83m/R27mxWJJwzRwU7XlW7\naD/++ONkZGSQmZlJjx49KC4ujmYuX/FjXy0WvD5O558PH3wA8+aZh1yyssyCU+++W/0pggcOwPz5\nZlu0adPMNmL33AOLFgWimt2vvH5OhVLtoj1y5EhWrlxJYWEh1113HU899VQ0c/lKYWGh7Qie4Jdx\nSkuDyZPNlXa/fqZlct55ZobHkiWmxbF7d+iHdA4dMvPBR4+Gyy+HM8+EJ54ws1UWL4Z27czv88tY\nxZofx6naC0adfvrpR/573759nHHGGVEJ5Ee7d++2HcET/DZO9evDLbeYj88+M9ua/frXZhuw7dtN\ngT7jDLMqX3Kymfv96aeQmgrdu8Ojj5rNGBo0OPG9/TZWseLHcYpolb/HHnuM1157jXr16rF06dJo\nZRLxnU6d4IUXjv3awYM/FfAdO8yVd5cu0KiRnYziDVW2R3Jzc+nYseMJHzNnzgRgzJgxbNy4kby8\nPB544AFHAnvR+vXrbUfwhHgbp7p1zROVnTqZJWB79w6/YMfbWFWXL8cpGAUbNmwIpqenh/y1jIyM\nIKAPfehDH/o4hY+MjIyQNbXa7ZF169bRtm1bAGbMmEFWVlbI3+fHGwEiIrZU++GaAQMGsGbNGmrU\nqMF55529FwQJAAACwElEQVTHc889R3JycrTziYjIUWL+RKSIiESPnoh02IQJE0hMTOSf//yn7Siu\nNGLECFJTU8nIyOD6669nz549tiO5SkFBAe3bt6dt27aMHz/edhzXKi4upnv37qSnp9OhQwf+8z//\n03akqFHRdlBxcTFz586ldevWtqO4Vs+ePSkqKmLlypW0a9eOcePG2Y7kGmVlZdxzzz0UFBTw5Zdf\nMn36dL766ivbsVypVq1aTJw4kaKiIpYuXcrUqVN9M1Yq2g568MEH+f3vf287hqvl5uaSmGhOy+zs\nbDZt2mQ5kXssW7aMNm3akJKSQq1atRg0aBAzZsywHcuVmjVrRmZmJgANGjQgNTWVLVu2WE4VHSra\nDpkxYwYtW7bkggsusB3FM15++WWuvvpq2zFcY/PmzbQ6aqucli1bsnnzZouJvGH9+vWsWLGC7Oxs\n21GiIqInIuVYubm5bN269YSvjxkzhnHjxjFnzpwjX4vn+7+VjdPYsWPp27cvYMasdu3a3HjjjU7H\nc62EhATbETxn3759DBgwgMmTJ9Mg1HoAHqSiHUVz584N+fXVq1fz3XffkZGRAcCmTZvo1KkTy5Yt\ni8tpkpWNU4X8/Hxmz57N/PnzHUrkDS1atDhmNc3i4mJatmxpMZG7lZSU0L9/fwYPHsx1111nO07U\naMqfBeeccw6fffYZjRs3th3FdQoKCnjooYdYuHChFiE7TmlpKeeffz7z58/nrLPO4pJLLmH69Omk\npqbajuY6wWCQoUOH0qRJEyb6bMdj9bQt0I+5lbv33nvZt28fubm5ZGVlcdddd9mO5Bo1a9ZkypQp\n9OrVi7S0NAYOHKiCXYklS5Ywbdo0FixYQFZWFllZWRQUFNiOFRW60hYR8RBdaYuIeIiKtoiIh6ho\ni4h4iIq2iIiHqGiLiHiIiraIiIeoaIuIeIiKtoiIh/w/wGB14zTfgpEAAAAASUVORK5CYII=\n", "text": [ "" ] } ], "prompt_number": 29 }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] } ], "metadata": {} } ] }