{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Plotting Approximate Stability Regions\n",
"\n",
"Copyright (C) 2020 Andreas Kloeckner\n",
"\n",
"\n",
"MIT License
\n",
"Permission is hereby granted, free of charge, to any person obtaining a copy\n",
"of this software and associated documentation files (the \"Software\"), to deal\n",
"in the Software without restriction, including without limitation the rights\n",
"to use, copy, modify, merge, publish, distribute, sublicense, and/or sell\n",
"copies of the Software, and to permit persons to whom the Software is\n",
"furnished to do so, subject to the following conditions:\n",
"\n",
"The above copyright notice and this permission notice shall be included in\n",
"all copies or substantial portions of the Software.\n",
"\n",
"THE SOFTWARE IS PROVIDED \"AS IS\", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR\n",
"IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,\n",
"FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE\n",
"AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER\n",
"LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,\n",
"OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN\n",
"THE SOFTWARE.\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false,
"jupyter": {
"outputs_hidden": false
}
},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as pt\n",
"\n",
"from cmath import exp, pi"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false,
"jupyter": {
"outputs_hidden": false
}
},
"outputs": [],
"source": [
"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)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false,
"jupyter": {
"outputs_hidden": false
}
},
"outputs": [],
"source": [
"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)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false,
"jupyter": {
"outputs_hidden": false
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEACAYAAABWLgY0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xu8VXP+x/HXRyGXxolSQyoi3YZSmibGHJJC5DJM+TE1\n/DAzv2HCuIzpNzFyH5coP4Qmg4kxxnXcSntCVNIh6YqT0oWQcUsX398f33N05OxzVu3Ld+213s/H\nYz/OXvus9vr0afXZa3/Wd32XOecQEZFk2iJ0ACIiUjgq8iIiCaYiLyKSYCryIiIJpiIvIpJgKvIi\nIgmWlyJvZnea2Qozez3L739iZqvM7NWqx7B8bFdEROrWME/vMxa4Gbi7jnUmO+eOztP2REQkgrwc\nyTvnXgA+rmc1y8e2REQkumL25Hua2Uwze8LMOhZxuyIiqZWvdk19ZgCtnXNfmNnhwMNAuyJtW0Qk\ntYpS5J1zn9V4/qSZ3WJmOzrnPtp4XTPTZDoiIpvIOVdrSzyf7RojS9/dzJrXeN4DsNoKfDXnXM6P\n4cOH5+V90vBQrpQn5am0c1WXvBzJm9l9QDmwk5m9CwwHtvL12t0O/NTMfgWsBb4EfpaP7dalsrKy\n0JtIDOUqGuUpGuUpumLkKi9F3jl3Uj2/Hw2Mzse2REQkusRe8TpkyJDQIZQM5Soa5Ska5Sm6YuTK\n6uvnFJuZubjFJCISZ2aGK8KJ11jJZDKhQygZylU0ylM0ylN0xchVYou8iIioXSMiUvJS2a4REZEE\nF3n1BaNTrqJRnqJRnqJTT15ERHKinryISIlTT15EJKUSW+TVF4xOuYpGeYpGeYpOPXkREcmJevIi\nIiVOPXkRkZRKbJFXXzA65Soa5Ska5Sk69eRFRCQn6smLiJQ49eRFRFIqsUVefcHolKtolKdolKfo\n1JMXEZGcqCcvIlLi1JMXEUmpxBZ59QWjU66iUZ6iUZ6iU09eRERyop68iEiJU09eRCSlElvk1ReM\nTrmKRnmKRnmKTj15ERHJSV568mZ2J9AfWOGc2yfLOjcBhwOfA0OccxVZ1lNPXkRkE9TVk2+Yp22M\nBW4G7s4SwOFAW+fcXmb2Q+BWoGeeti2SM+fgo49gyRL/87PPvv0AaNjw248mTaBZM2ja1D+aNIEt\n9N1YYiYvRd4594KZta5jlQFUfQA456aa2Q5m1tw5tyIf269NJpOhvLy8UG+fKGnK1YcfwqxZGx5v\nvQWLF/vivvXW0LKlL9jbb7/hsd12YAaVlRlatChn3TpYswZWrYIPPoCVK/3js8+gdWvYay/Yc0//\nc++9oVs32Gmn0H/z4knT/pSrYuQqX0fy9dkVWFxj+b2q1wpW5EXWr4fXXoPJk/1j6lRfiDt3hh/8\nALp0gRNOgN1288V9++3rfr9MBur6//jVV/DOO7BwISxYAHPmwEMPwYwZsPPO0KMH7L8/9OoF3btD\ngwb5/NuK1C5v4+SrjuQfq60nb2aPA1c456ZULU8AznfOzaxlXfXkZbMtWQKPPgqPPw5TpsCuu8JB\nB/nHj37kj7St1s5l4axfD/PmwfTpMG2a/8BZuhR694a+feGww/wHjcjmKkZPvj5LgJq7cUtgabaV\nhwwZQps2bQAoKyujS5cu33ylqR5ypGUtVy8vWgRLlpTz8MMwb16Gnj3h9NPLGTcOZs8OH1/1cseO\n0Lp1hhNOgHbtynnmGfjrXzOcey60bVvOoEHQqlWG5s3jEa+W47tc/byyspJ6Oefy8gDaALOy/O4I\n4Imq5z2Bl+t4H5cPkyZNysv7pEEp5mrVKuduvdW5Hj2c22UX584+27nnnnNuzZrCbbNQeVq3zrmJ\nE507/XTndtrJuQMOcG7UKOc+/rggmyu4UtyfQslXrqrqZq01NS9jAczsPmAK0M7M3jWzX5jZmWZ2\nRlXV/hfwjpktBG4Dfp2P7Ur6vPQSnHyyb7tMmACXXAKLFsHIkXDwwbDllqEj3HQNGsAhh8Dtt/s2\nzkUX+ZbO7rvDmWfC66+HjlBKmeaukdj7+mt47DG49lpfBM8+2xf6pk1DR1ZYy5bBmDFw223Qti2c\ncw4MGKBhmvJddfXkVeQlttatg3HjfHFv3BjOPx+OO86PUU+TtWvh4Yfhqqt8Tv73f30eVOylWion\nKKt5gkLqFrdcOQf/+Ad06gT33gu33upHpZx4YtgCHypPW27ph3q+8gpcfjlcfTXsuy/8/e8+V3ET\nt/0pzoqRq8QWeSlNkyZBz54wYgTcdBNMnOjHphd72GMcmUH//v4D7+qr4cor4cc/9uPwRbJRu0Zi\nYelSGDrUjyW/4gr42c/UjqjP+vXwl7/AsGFw+OE+by1ahI5KQkhlu0ZKw/r1MGqUbz+0awdvvgmD\nBqnAR9GgAZx2Gsyd66dN6NwZ/u///IlqkWqJ/a+kvmB0oXL1xhu+NfP3v/shgyNGwDbbBAklkrju\nUzvs4E9OT57sj+z79PHTK4QS1zzFkXrykkjOwc03+3Htv/ylnxOmQ4fQUZW+jh3hxRehXz8/T84t\nt+ioXtSTlyJbsQKGDPHT+d57r5+tUfJv7lwYPNhPjDZuHOy4Y+iIpJDUk5dYmDjRz/zYrRu88IIK\nfCG1b+9zvNdePt/Tp4eOSEJJbJFXXzC6QufKOT/twH/9lz96HzGiNKcfKLV9asst4frr4brr4Mgj\nYfTo4oyrL7U8hVSMXKXs2kEpttWrfd995kx4+WWomlxUiui442CffeD44/08OKNHp++q4TRTT14K\n5v334eijoVUrGDvW32FJwvn0U3/VsHPwwAPwve+FjkjyRT15KbpFi/zVmH36wP33q8DHQePGfqK3\n3Xf3/zZLloSOSIohsUVefcHo8p2rN9/0ReTXv4bLLkvOlARJ2KcaNvRDK08+2d+GcMGC/G8jCXkq\nFvXkpeRMnw5HHeUvzjnllNDRSG3M/IyeTZr4eewnTPA3HJdkUk9e8ubVV/2FOHfc4XvxEn9jx/q5\nbyZM0AVppSwO93iVhJs9G444wk8LrAJfOn7xC9/C6d3bX8egQp886slLzrlauBD69oU//9kP10uq\npO5Tp5zipy3u1w/eey/390tqngpBPXmJvWXL/AiaP/7Rn8yT0jR4sP+3PPxwP9FZWVnoiCRf1JOX\nzfbll/6GHv37+1vSSWlzzt8/94034KmnYOutQ0ckUeker5J3zm2Y9/3ee5MzTDLt1q+HgQN9gf/r\nX/XvWipSeTGU+oLRbU6u/vQnqKyEO+9MTyFIwz7VoIGftXL2bD8d9OZIQ57yRT15iaVHHoG77oKp\nU+N9kw/ZPNtuCw895G/o0rWrv7BNSpfaNbJJFi+G7t3h4YfhRz8KHY0U0lNP+dsLTp8Ou+wSOhqp\nSyrbNZJ/69b5Pvw556jAp0G/fn4G0ZNO0h2mSllii7z6gtFFzdWll/qv8hdcUNh44iqN+9TFF/sP\n95Ejo/+ZNOZpc6knL7ExZYqfrqCiwo+okXSoPhHbs6e/4K1jx9ARyaZST17q9dVX/gTcpZfCCSeE\njkZCGDPGT1nx8suleVevpFNPXnJyxRXQrh389KehI5FQ/vu/oUULuOaa0JHIpspLkTezfmY218zm\nm9mFtfx+sJm9b2avVj1Ozcd266K+YHR15eqNN/z846NHp2c8fDZp3qfM/D5www3+hjB1SXOeNlUx\ncpVzkTezLYBRQF+gEzDIzNrXsup459x+VY+7ct2uFJ5zfnTFZZfBrruGjkZCa9MGfvtbP7pKSkfO\nPXkz6wkMd84dXrV8EeCcc1fXWGcw0N05d1aE91NPPiYefBBGjIAZM/wJOJHVq6FzZ39U37dv6Gik\nWqF78rsCi2ssL6l6bWPHmVmFmT1gZi3zsF0poDVr4MIL4brrVOBlg0aN4Kab/ERm69aFjkaiyMcQ\nyto+PTY+FH8UuM85t9bMzgTGAb2zveGQIUNo06YNAGVlZXTp0oXy8nJgQw+rvuXq16Kun+bliooK\nhg4d+q3fV1SU0749NGiQIZOJV7yhljfet0LHE2p5222hZctyxo2Dtm2/+/va9qc4xR+n5RtvvHGz\n61smk6GyspJ6OedyegA9gadqLF8EXFjH+lsAq+r4vcuHSZMm5eV90mDjXH38sXPNmjk3e3aYeOJK\n+9QGL77oXKtWzq1e/d3fKU/R5StXVXWz1pqaj558A2Ae/sh8GTANGOScm1NjnRbOueVVz48FznfO\n9cryfi7XmCQ3l18O8+f7i2BEsjnySH/Lx//5n9CRSMHnkzezfsBI/FH6nc65q8zsUmC6c+5xM7sC\nOBpYC3wE/Mo5Nz/Le6nIB/TFF7D77pDJ6H6fUrcZM+Coo/ztH7fdNnQ06Vbwi6Gcc0855/Z2zu3l\nnLuq6rXhzrnHq55f7Jzr7Jzr6pzrna3A51PN3pXUrWau7rgDDjxQBb422qe+rVs3/7jnnm+/rjxF\nV4xc6YpX+caaNXDttfD734eORErFOef4ycv05Tu+NHeNfOO++/ydniZODB2JlArnYN99/VDbPn1C\nR5NemrtGIrn9dvjVr0JHIaXEDIYOhRtvDB2JZJPYIq++YHSZTIa5c2HuXBgwIHQ08aV9qnYnnQSv\nvAJvveWXlafo1JOXohkzBn7xC00jK5uuUSMYOPC7J2AlHtSTF9as8ROQTZ0Ke+wROhopRa+84gv9\nggWarTQE9eSlThMmQPv2KvCy+bp1898CX3opdCSyscQWefUFo7v55oxuCBKB9qnszODnP4e771ae\nNoV68lJwa9fCiy/CcceFjkRK3cCB8M9/wtdfh45EalJPPuWefhouuURfsyU/OnWCsWOhR4/QkaSL\nevKS1cMP6yhe8ufII+GJJ0JHITUltsirLxjNxInQpEkmdBglQftU/Y48EsaPz4QOo2SoJy8FtXgx\nfPyxRtVI/vTqBUuXwvLloSORaurJp9i4cf6r9QMPhI5EkuToo/1IG43YKh715KVWEyfCIYeEjkKS\nplcvmDIldBRSLbFFXv3T+r34IvzkJ8pVVMpTNI0aZVTkI1JPXgpm1Sp4/31o1y50JJI07dvDrFmw\nenXoSATUk0+tTAYuvlhfq6UwuneHm27yrRspPPXk5TtmzoT99gsdhSTVD34As2eHjkIgwUVe/dO6\nvfoqdO3qnytX0ShP0WQyGTp0gDlzQkcSf+rJS8HMmwcdO4aOQpJKRT4+1JNPqWbN/MmxFi1CRyJJ\ntHAhHHooVFaGjiQd6urJq8in0Gefwc47w+ef6wYPUhjr18N22/lRXI0ahY4m+VJ54lX90+wqK6F1\n6w0FXrmKRnmKJpPJ0KCBP5BYsSJ0NPGmnrwURGUltGkTOgpJuhYtYNmy0FGI2jUpdNdd8Pzzft5v\nkUIZMMDfHP6YY0JHknypbNdIdqtWQVlZ6Cgk6XQkHw+JLfLqn2a3cZFXrqJRnqKpzlOzZrByZdhY\n4q5kevJm1s/M5prZfDO7sJbfb2Vm481sgZm9ZGat8rFd2TwffwxNmoSOQpJum23gyy9DRyE5F3kz\n2wIYBfQFOgGDzKz9RqudBnzknNsLuBG4Jtft1qe8vLzQmyhZn3wCO+ywYVm5ikZ5iqY6T40aaZKy\n+hRjn8rHkXwPYIFzbpFzbi0wHhiw0ToDgHFVzx8Eeudhu7KZ1q2DLbcMHYUknYp8POSjyO8KLK6x\nvKTqtVrXcc6tB1aZ2Y552HZW6p9Gp1xFozxFU52nlSs1y2l9irFPNczDe9Q2bGfjMZAbr2O1rPON\nIUOG0KZqIHdZWRldunT55mtNdVLqW64Wdf00LfsLVDYsV1RUxCo+LZf2cvX+9M478NprGTKZeMUX\np+WKiorN+vPVzysjzBuR8zh5M+sJXOKc61e1fBHgnHNX11jnyap1pppZA2CZc27nLO+ncfIFdtJJ\n0L+//ylSKLfd5mc7ve220JEkX6HHyU8H9jSz1ma2FTAQeHSjdR4DBlc9PwF4Lg/blRzoc1QK7csv\nNW9NHORc5Kt67L8BngFmA+Odc3PM7FIz61+12p1AUzNbAAwFLsp1u/Wp+bVGvq1xY/j00w3LylU0\nylM01XlavVpFvj7F2Kfy0ZPHOfcUsPdGrw2v8fwr4MR8bEtyV1bmL4gSKSQV+XjQ3DUpdOWVfqz8\nVVeFjkSS7JxzoGVLOO+80JEkn+aukW9p0sRf9SpSSMuWwfe/HzoKSWyRV/80uyZN4MMPNywrV9Eo\nT9FU52n5ct15rD7F2KcSW+Qlu1at4N13Q0chSbd8uY7k40A9+RRavhz22Qfefz90JJJkO+zgb1Cj\nyfAKTz15+Zbmzf19Xj/7LHQkklTV7UDdtyC8xBZ59U+zM/P3eF20yC8rV9EoT9FkMhnmzIEOHXSj\n+PqoJy8F07YtLFgQOgpJquoiL+GpJ59Sf/iDn274kktCRyJJdO65/qTr+eeHjiQd1JOX7+ja1U8e\nJVIIb76pI/m4SGyRV/+0bvvtBzNn+ufKVTTKUzSTJmWYPt3vY1I39eSlYHbf3U9S9sEHoSORpFm8\nGL73Pdhll9CRCKgnn2p9+sBZZ8HRR4eORJJk7FiYMAHuvTd0JOmhnrzU6uCD4TnN7C95NmUK9OoV\nOgqpltgir/5p/Xr39kVeuYpGeYrm6aczKvIRqScvBdWtm5/DRjNSSr68/ba/knrffUNHItXUk0+5\no4/293odODB0JJIEo0bBjBm+Ly/Fo568ZHXEEfDoxnfkFdlMTzwBRx4ZOgqpKbFFXv3TaI49Fh59\nNMPq1aEjiT/tU3X7/HN48UVo1CgTOpSSoZ68FFzz5rDnnvDss6EjkVL35JPwwx/C9tuHjkRqUk9e\nGDUKpk+HceNCRyKlbMAAOO44GDw4dCTpU1dPXkVeWLoUOneG996DbbYJHY2Uog8+gL328le7Nm4c\nOpr0SeWJV/VPo5s/P0PPnvDgg6EjiTftU9ndf78/4dq4sfK0KdSTl6I54wy47bbQUUipuvtuOOWU\n0FFIbdSuEQDWrfN3i3rmGejUKXQ0UkqmTYMTT4SFC6Fhw9DRpFMq2zWyaRo2hNNOg9tvDx2JlJqR\nI+Hss1Xg4yqxRV59weiqc3X66XDPPZrmIBvtU9/13nt+6ORpp214TXmKTj15KarddvPTHIweHToS\nKRWjR8PJJ8MOO4SORLLJqSdvZk2A+4HWQCVwonPuk1rWWw+8BhiwyDl3TB3vqZ58QHPnwkEHwTvv\nwHbbhY5G4mzVKj9scsoU/1PCKWRP/iJggnNub+A54PdZ1vvcObefc65rXQVewmvf3hf5MWNCRyJx\nd8MNftikCny85VrkBwDV10mOA7IV8Fo/YQpJfcHoNs7VxRfDn/8MX3wRJp640j61wYcf+iul//jH\n7/5OeYquFHryOzvnVgA455YDzbKst7WZTTOzKWY2IMdtSoHttx8ccABcf33oSCSurr0WTjgB9tgj\ndCRSn3oHPZnZs0Dzmi8BDhi2Cdtp5Zxbbma7A8+Z2evOuXeyrTxkyBDatGkDQFlZGV26dKG8vBzY\n8Mmn5fwuV6tevvLKcvbfHzp2zLDjjuHji8NyeXl5rOIJtbxyJYwZU85rr0Xfn+IUf5yWq1/bnP+v\nmUyGyspK6pPridc5QLlzboWZtQAmOec61PNnxgKPOeceyvJ7nXiNifPO83f50ZWwUtNJJ/kj+BEj\nQkci1Qp54vVRYEjV88HAI7VsvMzMtqp63hToBbyZ43brtfERhWSXLVfDhsE//wmzZhU3nrjSPgWT\nJvk54y++OPs6ylN0xchVrkX+aqCPmc0DDgWuAjCzbmZWfe1kB+AVM5sJTASudM7NzXG7UgRNmsBl\nl/mLpNavDx2NhLZ2LfzmN35Uzbbbho5GotLcNVKnr7+Ggw+G44/3l65Lel13nb+5zJNPghV9vJzU\nRfPJS07mz4devfwNmlu3Dh2NhDBnjr9+4qWX/J3EJF5SOUGZ+oLR1Zerdu38Sdgzz4Q0f/6mdZ9a\nuxZ+/nN/ojVKgU9rnjZHKfTkJSV+9zs/cdnIkaEjkWK74gpo2tTfc0BKj9o1Etnbb0PPnvCvf0H3\n7qGjkWKYNg2OOgpmzoRddgkdjWSTynaN5N8ee/hZBwcOhP/8J3Q0UmgrV/qbgdx6qwp8KUtskVdf\nMLpNydUJJ0Dv3v6re9q+cKVpn1q3zn+YDxoExx67aX82TXnKlXryEks33ghvvQVXXhk6EimUYcP8\nMEld1Vr61JOXzbJ0Kfzwh77gH3986Ggknx54AC64AF55xZ9wlfjTOHkpiBkzoF8/eOop6NYtdDSS\nD//+t2/JPfMMdOkSOhqJKpUnXtUXjG5zc9Wtm5+8bMAAfyeppEv6PjVrlj/ROn58bgU+6XnKp2Lk\nSvdXl5wcdxwsWwaHHgqTJ8Ouu4aOSDbH4sX+Lk8jR8Ihh4SORvJJ7RrJi2uugbFj/df9nXcOHY1s\niqVLfWE/4ww499zQ0cjmqKtdoyN5yYsLLvBzzx92GDz3HOy4Y+iIJIolS3yBP/VUFfikUk9e8par\nSy+FPn3gJz/xLZykSdo+9e67/t/qjDPgoovy975Jy1MhaZy8lBQz37YZOBAOPNCPpZd4WrgQysvh\nrLP8vESSXOrJS0Hcequ/4ciTT8I++4SORmqaMsVf2/CnP/kbwkjpU09eiu6Xv/R9+UMPhbvv9uPp\nJbwHHvB3d9K/SXoktl2jvmB0hcrViSf6e8Seeipce23pz3VTyvuUc3D11f6+AM8+W9gCX8p5KjaN\nk5eSd8ABMHUqHHMMVFTAHXfANtuEjipdPvnEf9C++66/s1PLlqEjkmJST16K4osvfP/3zTfhb3+D\n9u1DR5QOM2f6aQr69fP3aN1669ARSSGkcloDiZdtt4V77vG9+h//2E+HoM/ywnHOn/w+7DC4/HIY\nNUoFPq0SW+TVF4yuWLky8/eJff55X+SPOQY++KAom86LUtmn3n0X+vb1rbHnn4ef/ay42y+VPMWB\nxslLIrVvDy+/7H/us48/wtdRfe6cgzFj/MRxBx+8IceSburJS1DTpvmj+512gltugXbtQkdUmt58\nE84+259kHTsWOncOHZEUk3ryEls9esD06dC/P/TqBcOH+zlwJJqPPoLf/tZPT9C/vx89owIvNSW2\nyKsvGF3oXDVsCEOH+pEg8+f7o/lbboG1a4OG9R2h81TT2rX+puodOvjnc+b4HDaMwaDoOOUp7tST\nl1TZbTc/vPLxx+GRR3wBGz8evv46dGTxsWYN3H67/yB85BF/YdMtt+g2fZKdevISW889BxdfDKtW\n+Ss1TzkFGjUKHVUYX30Ff/kLXHGFP5k6fLhvb4lAAXvyZvZTM3vDzNab2X51rNfPzOaa2XwzuzCX\nbUp6HHKI7zHfeis8/DC0aQMjRsCHH4aOrHgWL4Zhw6B1a5+D+++Hp59WgZfocm3XzAKOBf6dbQUz\n2wIYBfQFOgGDzKzgA7vUF4wuzrky81PiPvEETJwIb78NbdvCoEG+VbF+ffFiKVaevv7a/12PP97f\na/U//4FMxs/o2bNnUULISZz3p7iJfU/eOTfPObcAqPVrQpUewALn3CLn3FpgPDAgl+1KOnXqBHfd\n5Qv9gQf6G13svjv88Y8we3Zpj7V3zo8yOu88aNXK/zz0UKishJtu0nh32Xx56cmb2STgPOfcq7X8\n7nigr3PujKrlk4Eezrmzs7yXevIS2Wuv+V71Qw/Bllv6q2gHDPDtjAYNQkdXt6++8nO7P/00PPig\n/9YyaJC/6UrHjqGjk1KS03zyZvYs0LzmS4AD/uCceyzK9mt5TVVc8mLffeGGG+D66/0sl4884u92\ntGQJHHTQhse++4Yv+uvW+W8ckyf7wj55sh9BdNhhfhRRt26+0IvkU71F3jnXJ8dtLAFa1VhuCSyt\n6w8MGTKENm3aAFBWVkaXLl0oLy8HNvSw6luufi3q+mlerqioYOjQobGJZ3OXu3aFTz7JUF4Oe+5Z\nzuTJMH58hhtugE8+KadHDygry7DHHnDiieV06ABTp0Z//433rbrW79atnIUL4R//yDB3LixfXk5F\nBey4Y4bOnWHw4HLGjYNZs/z63buHz1++lpOyPxVj+cYbb9zs+pbJZKisrKQ++WzX/M45N6OW3zUA\n5gG9gWXANGCQc25OlvfKS7smk8l8kxipWxpy9f77vuc9axa8/rr/uXAhtGjh51ffbbcNP3faCbbf\n3j8aN/Y/AV56KUPXruWsW+cvQProI1i50j8++ABWrPD3tV240J8sbdsW9t4buneH/ff3R+o77BA2\nD8WQhv0pX/KVq7raNTkVeTM7BrgZaAqsAiqcc4eb2feBMc65/lXr9QNG4k/03umcu6qO91RPXopi\nzRo/RHHxYt/eqX7+8cd+aoXqx6ef+jZKw4bffjRpAs2a+QuRmjaFnXeGPfaAvfaCXXaBLXSpoRRJ\nwYp8IajIi4hsmlROUFazdyV1U66iUZ6iUZ6iK0auElvkRURE7RoRkZKXynaNiIgkuMirLxidchWN\n8hSN8hSdevIiIpIT9eRFREqcevIiIimV2CKvvmB0ylU0ylM0ylN06smLiEhO1JMXESlx6smLiKRU\nYou8+oLRKVfRKE/RKE/RqScvIiI5UU9eRKTEqScvIpJSiS3y6gtGp1xFozxFozxFp568iIjkRD15\nEZESp568iEhKJbbIqy8YnXIVjfIUjfIUnXryIiKSE/XkRURKnHryIiIpldgir75gdMpVNMpTNMpT\ndOrJi4hITtSTFxEpcerJi4ikVE5F3sx+amZvmNl6M9uvjvUqzew1M5tpZtNy2WZU6gtGp1xFozxF\nozxFVwo9+VnAscC/61nva6DcOdfVOdcjx21GUlFRUYzNJIJyFY3yFI3yFF0xctUwlz/snJsHYGa1\n9oJqMIrcGlq1alUxN1fSlKtolKdolKfoipGrYhVeBzxtZtPN7PQibVNEJPXqPZI3s2eB5jVfwhft\nPzjnHou4nV7OueVm1gx41szmOOde2PRwo6usrCzk2yeKchWN8hSN8hRdMXKVlyGUZjYJOM8592qE\ndYcDnzrnrs/ye42fFBHZRNmGUObUk99IrRsws22BLZxzn5nZdsBhwKXZ3iRboCIisulyHUJ5jJkt\nBnoCj5u/c3A6AAACLElEQVTZk1Wvf9/MHq9arTnwgpnNBF4GHnPOPZPLdkVEJJrYXfEqIiL5k+gr\nXs3sTzUuwnrKzFqEjimOzOwaM5tjZhVm9g8z+17omOIq6gWAaWVm/cxsrpnNN7MLQ8cTV2Z2p5mt\nMLPXC72tRBd54Brn3L7Oua7AE8Dw0AHF1DNAJ+dcF2AB8PvA8cRZ1AsAU8fMtgBGAX2BTsAgM2sf\nNqrYGovPU8Elusg75z6rsbgd/spb2YhzboJzrjo3LwMtQ8YTZ865ec65BWQZaJByPYAFzrlFzrm1\nwHhgQOCYYqlqCPnHxdhWPkfXxJKZjQB+DqwCDg4cTik4Ff+fU2RT7QosrrG8BF/4JaCSL/L1Xazl\nnBsGDKvqD54FXFL8KMOLclGbmf0BWOucuy9AiLGRpwsA06i2bzca2RFYyRd551yfiKv+Dd+Xv6Rw\n0cRXfXkys8HAEcAhxYkovjZhn5JvWwK0qrHcElgaKBapkuievJntWWNxADAnVCxxZmb9gAuAo51z\nX4WOp4SoL/9t04E9zay1mW0FDAQeDRxTnBlF2IcSPU7ezB4E2uFPuC4CfumcWxY2qvgxswXAVsCH\nVS+97Jz7dcCQYsvMjgFuBpriz/NUOOcODxtVfFQdMIzEH0De6Zy7KnBIsWRm9wHlwE7ACmC4c25s\nQbaV5CIvIpJ2iW7XiIiknYq8iEiCqciLiCSYiryISIKpyIuIJJiKvIhIgqnIi4gkmIq8iEiC/T/C\n1aMREJkloQAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"def fw_euler_step(y, t, h, f):\n",
" return y + h * f(t, y)\n",
"\n",
"plot_stability_region(-1, fw_euler_step)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false,
"jupyter": {
"outputs_hidden": false
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEACAYAAACwB81wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcFNW1wPHfGXZQGcWABJXxQRSIy7gRVGImigooIS5x\nj46JazSaGH0u8UWMJi5PnxrjEp9LXEIMT9EIboAyKlERFxaRTXQAw6YiioDIDOf9cXukwR7omaqu\n23X7fD+f/tDVXU6d44Uz1adv3RJVxRhjTPjKfAdgjDEmGVbwjTGmRFjBN8aYEmEF3xhjSoQVfGOM\nKRFW8I0xpkRELvgisr2IvCAi74rINBE5v5H9/iQic0RksohURj2uMcaYpmkZw8+oAy5U1ckisgXw\npoiMUdWZDTuIyCCgh6p+R0S+B9wF9Ivh2MYYY/IU+QxfVRer6uTM8y+AGUC3jXYbCjyY2Wci0FFE\nukQ9tjHGmPzF2sMXkQqgEpi40VvdgAVZ2//mm78UjDHGFFBsBT/TznkUuCBzpr/B2zn+E1vTwRhj\nEhRHDx8RaYkr9g+p6j9z7PIhsEPW9vbAwkZ+lv0iMMaYJlLVXCfWG4jrDP8+4F1VvbWR958ETgEQ\nkX7AclVd0tgPU9UgH1deeaX3GCw/y8/yC++Rr8hn+CJyAHASME1E3sa1ai4HurvarXer6tMiMlhE\n3gNWAqdFPW4a1dbW+g6hoCy/dLP8whe54Kvqv4AWeex3XtRjGWOMaT670jZB1dXVvkMoKMsv3Sy/\n8ElT+j9JEBEttpiMMaaYiQia4Je2Jg81NTW+Qygoyy/dLL/wWcE3xpgSYS0dY4xJOWvpGGOM2YAV\n/ASF3kO0/NLN8gufFXxjjCkR1sM3ZhNWr4Zly6B1a2jb1j1atgTZbLfUmOTk28OPZfE0Y9Jq3TqY\nORNefRWmTYPFizd8rF4N22wDdXXw5ZfusW4dtGkD7dvDTjtB797Qq9f6R8+e7heEMcXGzvATVFNT\nQ1VVle8wCiYN+a1aBS+9BK+95or8xInQqRPstx9UVkK3brDddusf5eXrz+Yb8qurgzVrYOVKmDvX\n/cLIfsybB3vvDYMGweDB7ueWpaB5mobxiyLk/OwM35iMujoYNw6GD4dRo2D33eGAA+Dcc+Ghh6Bz\n56b9vJYt3aNDB/ff7rffhu+vXg0vvwxPPw0nnACffba++B9+uPtkYIwPdoZvgqTqzt7/9jcYMQIq\nKuCkk+C446BLwjfXnDsXnnnG/bJ580047TT4xS9cO8iYOOR7hm8F3wRFFcaMgauugqVL4ZRT4MQT\nXV+9GMydC3feCX/9K+y/P5x3HgwYkI6WjyleduFVEQp9HrDP/FTdWfR++8Gvfw2//CXMmgW/+118\nxT6O/Hr0gBtvhPnz4Uc/gosvhl13hdGjXQ4+2d/P8FnBN6k3diz06+eK569/7WbbnHACtNjsXRr8\nad8eTj8dJk+Gm25ysQ8cCNOn+47MhMxaOia1PvnEFfgJE+D66+Hoo9PbGlm71rV6rrkGjj3WtaQ6\ndfIdlUkLa+mYoD36KOy2m5sjP20a/OQn6S32AK1awfnnw4wZbhpo795uBpExcUrxP5H0Cb2HmER+\nixe7M/krrnBF/5Zb3PTIJCSRX6dOcNttrk31xz9CdTV88UXBDwvY389SEEvBF5F7RWSJiExt5P0f\niMhyEXkr87gijuOa0lJT4y5i2mUX1/vef3/fERXOHnvAG2+4s/199oEpU3xHZEIQSw9fRPoDXwAP\nquruOd7/AfAbVf1RHj/LevjmG+6+G/7rv9y8+gEDfEeTrIcfdt9VXHUVnHOOreNjvinRK21VdYKI\ndN9cTHEcy5SWujq48EI3t37CBPjOd3xHlLyTT4a+fd2XuVOnwu23F/cMJFO8kuzh9xORt0XkKRHp\nk+Bxi0boPcS48/v0U7ccwezZbu0b38Xe5/jtvLNbrmH2bHfF8FdfxX8M+/sZvqQK/ptAd1XdE/gz\n8ERCxzUp9cknUFXlZquMHu0WMSt1W27p1uf58kt30dbKlb4jMmkT2zz8TEtnVK4efo59PwD2VtVl\nOd7TU089lYqKCgDKy8uprKz8epW7ht/Sth3u9hdfwLBhVRx8MAwcWINIccXne7u+Hh5+uIpZs+Cy\ny2rYcsviis+2C7/d8Ly2thaABx54INm1dESkAlfwd8vxXhdVXZJ53hcYoaoVjfwc+9K2hK1YAYce\n6nrWt9xiX1A2Zt06uOgiGD8eXnwRttrKd0TGp0QvvBKR4cArwM4iMl9EThORs0TkzMwux4jIOyLy\nNnALcFwcx02b7N/OIYqa38qVbvngPfYozmJfTONXVuaWZNh3X/dl7tq10X9mMeVXCKHnl4+4Zumc\nuJn3bwduj+NYJkx1dXDkkW5xsTvuKL5iX4xE3P+rIUPccst3323/38ym2Vo6pihcdBG88w489ZRN\nOWyqFSvgBz+AY46Byy/3HY3xwe54ZVJjxAgYOdJdWWrFvum23NLNZNpvP3ejlxM3+XnblDJbSydB\nofcQm5Pf9OnuVoOPPeYWQitmxTx+3/62+3R0wQXuvrrNUcz5xSH0/PJhBd9489lnrm9/442w556+\no0m/XXeFq692Z/iFuDDLpJ/18I03J5wAW2/tvng08VCFoUPdBWvXX+87GpMU6+GbojZqlOvZT825\nvqppLhG45x63qujAgfDDH/qOyBQTa+kkKPQeYr75ff6569vffTe0a1fYmOKUlvHr3BnuvdfdwP3T\nT/P/79KSX3OFnl8+rOCbxF1+ubua1s4+C2fQIPcYNsx3JKaYWA/fJOpf/3K3I5w+3fXvTeF89BH0\n6eOWXuhTkuvTlg67p60pOvX1cPbZcOutVuyT8K1vwW9/626eYudQBqzgJyr0HuLm8hs+HDp2dFeE\nplEax+/cc2H+fDdHf3PSmF9ThJ5fPqzgm0R89ZXrJ//hD7beS5JatYKbb3Z3DbO5+cZ6+CYRd93l\nlk8YM8Z3JKXpsMPcdyenn+47ElMI+fbwreCbglu92t2ecORIt869Sd6LL8IZZ8CMGbZeUYjsS9si\nFHoPsbH8/vIX2Gef9Bf7NI/fgQdCp07wxCZuLprm/PIRen75sIJvCqq+Hv70J7jsMt+RlDYRuPRS\nuO46m7FTyqylYwrqqafcl7Wvv25f1vq2bp1bYO222+Dgg31HY+JkLR1TFG6/3U0NtGLvX1kZXHwx\n/M//+I7E+GIFP0Gh9xA3zu+992DSJDgukDsYhzB+xx4Lr7wCixZ9870Q8tuU0PPLR1w3Mb9XRJaI\nSKNrH4rIn0RkjohMFpHKOI5ritudd8LPfpauBdJC16EDHH00PPSQ70iMD7H08EWkP/AF8KCq7p7j\n/UHAeap6uIh8D7hVVfs18rOshx+Aujro1g0mTHBTMk3xmDABzjzTrWdkrbYwJNrDV9UJwKYWYh0K\nPJjZdyLQUUS6xHFsU5xeegm2396KfTE64ABYu9a120xpSaqH3w1YkLX978xrJSX0HmJ2fiNGuH5x\nSEIZPxE49VR44IENXw8lv8aEnl8+krrjVa6PGo32baqrq6moqACgvLycyspKqqqqgPWDZtvFu11f\nDyNHVvHaa8URj21/c/uYY6oYMACOOaYGEf/x2HbTthue19bW0hSxzcMXke7AqEZ6+HcB41X1H5nt\nmcAPVHVJjn2th59y48a5C62sZVC8VGHnnd0nMbuBfPr5mIcv5D6TB3gSOCUTWD9gea5ib8Lw+OPp\nXQK5VIjAkCHu3sKmdMQ1LXM48Aqws4jMF5HTROQsETkTQFWfBj4QkfeAvwC/iOO4aZP9cSxEDfmN\nG+dWZwxNaOM3ZAiMHr1+O7T8NhZ6fvmIpYevqifmsc95cRzLFLf582HZMtj9G409U2z693cXxy1a\nBF27+o7GJMHW0jGxuv9+eO45eOQR35GYfBx1lFsn/4QTfEdiorC1dIwX48bBgAG+ozD56t/fXYhl\nSoMV/ASF3kMcP76GF14IdyXGEMcvu+CHmF+20PPLhxV8E5uPPnJL8GYuoTApsOee8P77sHy570hM\nEqyHb2IzciTcd9+GMz9M8TvoILjoIhg82Hckprmsh28SN2kS7Luv7yhMU/XrZxfJlQor+AkKvYc4\nZkxN0AU/1PHbfXeYOjXc/BqEnl8+rOCbWKjCrFmw996+IzFNtcceMGWK7yhMEqyHb2KxaJErHEuX\n+o7ENFVdHWy1lRu7LbbwHY1pDuvhm0TNmAG9evmOwjRHy5bQpw+8847vSEyhWcFPUMg9xJkzoWPH\nGt9hFFTI49enDzz+eI3vMAoq5PHLlxV8E4sZM2DHHX1HYZqrRw9YuNB3FKbQrOAnqOEmBiGaPRuO\nOKLKdxgFFfL49egBdXVVvsMoqJDHL19W8E0samvtCts069ED5s71HYUpNCv4CQq1h6jqlkWura3x\nHUpBhTp+4Ar+zJk1vsMoqJDHL19W8E1kH38MHTpAu3a+IzHN9a1vwZdfwsqVviMxhWTz8E1kb74J\nZ5wBb73lOxITRY8e7l4GPXv6jsQ0lc3DN4n58EPYfnvfUZiouna1mTqhs4KfoFB7iEuXQpcu4ebX\nIPT8WrasCbrghz5++YjrJuYDRWSmiMwWkUtyvH+qiCwVkbcyj5/FcVxTHJYscT1gk26dOrklMky4\nIvfwRaQMmA0cDCwEJgHHq+rMrH1OBfZW1fPz+HnWw0+Z88+H//gP+NWvfEdiorj6alizBq65xnck\npqmS7OH3Beao6jxVXQs8AgzNFVMMxzJFaOlSO8MPQadO8MknvqMwhRRHwe8GLMja/jDz2saOEpHJ\nIjJCREryK75Qe4iffuqKRaj5NQg9v8WLa4Iu+KGPXz5axvAzcp25b9yTeRIYrqprReQs4AFcCyin\n6upqKjKXbZaXl1NZWfn1ZdENg2bbxbM9fz507FjFmjXFEY9tN297q63gvfdqqKkpjnhsu/Hthue1\ntbU0RRw9/H7AMFUdmNm+FFBVvb6R/cuAZapa3sj71sNPmV694PHHoXdv35GYKCZNgnPOgTfe8B2J\naaoke/iTgJ4i0l1EWgPH487os4PZLmtzKPBuDMc1RWL5cijP+evbpMmWW8KKFb6jMIUUueCraj1w\nHjAGmA48oqozROQqETkis9v5IvKOiLyd2bc66nHTKPvjWEg+/9wVi1DzaxB6ftOm1QRd8EMfv3zE\n0cNHVZ8FdtnotSuznl8OXB7HsUxxWbfOrcHSvr3vSExU7dvbGX7obC0dE8nKlW5K5qpVviMxUdXX\nQ6tW7k+xSdSpYmvpmESsXOlWyjTp16KFK/hr1viOxBSKFfwEhdhDXLVq/bLIIeaXrRTya9vWtehC\nFPr45cMKvolkzRpo29Z3FCYu7drB6tW+ozCFYj18E8nUqXDSSTBtmu9ITBx22gmef96tjWTSw3r4\nJhFr1kCbNr6jMHFp1QrWrvUdhSkUK/gJCrGH+NVX0Lq1ex5iftlKIb/WrcMt+KGPXz6s4JtI6urc\nWaEJQ6tW7pe4CZP18E0k48bBtde6vq9Jv7594bbb4Hvf8x2JaQrr4ZtE1NVBy1iu1zbFoGVLd+GV\nCZMV/ASF2EPMLvgh5petFPJr0cKNaYhCH798WME3kaxb567QNGFo0cLO8ENmBT9BDTcxCEl9PZRl\n/haFmF+2Usgv5IIf+vjlwwq+iaS+3s7wQ1JW5j61mTBZwU9QiD3EdevWn+GHmF+2UsivRYtwC37o\n45cPK/gmEtX1Bd+kX1lZuC0dYwU/USH2ENetW792eoj5ZSuF/EJeBz/08cuHFXwTiardLCM0dt1j\nuKzgJyjEHmJ2SyfE/LKVQn4i4Rb80McvH7EUfBEZKCIzRWS2iFyS4/3WIvKIiMwRkVdFZMc4jmv8\nC7U4lKqQC76JoeCLSBnwZ+Aw4LvACSLSa6Pdfg4sU9XvALcAN0Q9bhqF2kO0Hn4YrIcfvjjO8PsC\nc1R1nqquBR4Bhm60z1DggczzR4GDYziuMaYA7Aw/XHEU/G7AgqztDzOv5dxHVeuB5SKyTQzHTpXQ\ne4iWX7rV1NQwaRIsXuw7ksIIffzyEcc6h7k+BG58jrDxPpJjn69VV1dTUVEBQHl5OZWVlV9/HGsY\nNNsuju0ZM2pYsgSgOOKx7WjbixfX8OKLcNZZxRGPbefebnheW1tLU0ReD19E+gHDVHVgZvtSQFX1\n+qx9nsnsM1FEWgCLVLVzIz/P1sNPkYcfhmefdX+a9PvRj+DnP4ehGzdlTVFLcj38SUBPEekuIq2B\n44EnN9pnFHBq5vlPgBdiOK4xpgBC/uK21EUu+Jme/HnAGGA68IiqzhCRq0TkiMxu9wLbisgc4FfA\npVGPm0bZH8dC0vCBLNT8GpRCfiF/uA59/PIRy72KVPVZYJeNXrsy6/ka4Ng4jmWKi50NhsfGNFx2\npW2CGr54CU3DWWGo+TUohfxCPsMPffzyYQXfRFJWZvO2Q2JrI4XNCn6CQuwhiqxfPz3E/LKVQn4h\nF/zQxy8fVvBNJHaGH5bsG9qY8NjQJijEHmJZ1g0zQswvWynkF/JN6UMfv3xYwTeRlNk9UINiZ/hh\ns6FNUIg9xOx7oIaYX7ZSyK++3nr4IbOCbyLJbumY9Kuvh5axXJ1jipEV/ASF2ENs0QLq6tzzEPPL\nVgr51ddbDz9kVvBNJK1a2Rl+SOrq7Aw/ZFbwExRiD7FlS1i71j0PMb9spZBfyAU/9PHLhxV8E0nL\nlutbOib91q6F1q19R2EKJfJ6+HGz9fDTZeJE+OUv4fXXfUdi4tCrFzz+OPTu7TsS0xRJrodvSljr\n1rBmje8oTFzsDD9sVvATFGIPsU2b9QU/xPyylUJ+a9aEW/BDH798WME3kWQXfJN+q1dDu3a+ozCF\nYj18E8nChbD33rBoke9ITBzat4ePPoIOHXxHYprCevgmEe3bw8qVvqMwcVCFL7+Etm19R2IKxQp+\ngkLsIXboAKtWuWIRYn7ZQs9vzJgaWrUK90rb0McvH5EKvohsLSJjRGSWiDwnIh0b2a9eRN4SkbdF\n5IkoxzTFpVUrt57OV1/5jsREtWoVbLWV7yhMIUXq4YvI9cAnqnqDiFwCbK2ql+bY73NVzeuvkvXw\n02ebbWDOHOjUyXckJoq5c+GQQ+D9931HYpoqqR7+UOCBzPMHgB83Fk/E45gi1rEjfPaZ7yhMVCtW\nwJZb+o7CFFLUgt9ZVZcAqOpi4FuN7NdGRF4XkVdEZGjEY6ZWqD3E8nJYvjzc/BqEnt9LL9UEXfBD\nH798bHaZJBEZC3TJfglQ4IomHGdHVV0sIjsBL4jIVFX9oLGdq6urqaioAKC8vJzKysqvlzZtGDTb\nLp5tVVi+vIqysuKIx7abt/3556BaQ01NccRj241vNzyvra2lKaL28GcAVaq6RES2A8ar6iZX4RCR\n+4FRqjqykfeth58yRx4JJ58MRx/tOxITxb33woQJcP/9viMxTZVUD/9JoDrz/FTgnzkCKReR1pnn\n2wL7A+9GPK4pIttuCx9/7DsKE9WyZfbFe+iiFvzrgUNEZBYwALgOQET2FpG7M/v0Bt4QkbeB54Fr\nVXVmxOOmUvbHsZB07gxLl4abX4PQ83v77ZqgC37o45ePSLc6UNVluEK/8etvAmdmnr8K7B7lOKa4\ndenipmWadFu2zI2lCZddaZughi9eQtNwhh9qfg1Cz0+1im9/23cUhRP6+OXDCr6JrGtXt4iaSbeF\nCwm64Bsr+IkKtYe4446wYEG4+TUIPb9582qCLvihj18+rOCbyLp1c8sj19f7jsQ01+rVbqXMkL+0\nNbYevolJt27w2muwww6+IzHNMX26u45iZknOn0s/Ww/fJGrHHWHePN9RmOaaOxd69PAdhSk0K/gJ\nCrmH2LMnjB5d4zuMggp5/ObOhTZtanyHUVAhj1++rOCbWPTubWf4aTZ3rs3QKQVW8BMU8jzgXr1g\n1aoq32EUVMjjN3s2DB5c5TuMggp5/PJlBd/EolcvmDHDdxSmuaZOhd3tevjgWcFPUMg9xJ49oba2\nhi+/9B1J4YQ6fkuWuFtUzplT4zuUggp1/JrCCr6JRevWbkrmtGm+IzFN1XB2L3ZfuuDZPHwTm9NP\nh732gl/8wnckpiluugnmz4dbb/UdiWkum4dvErfvvjBpku8oTFO9+SZUVvqOwiTBCn6CQu8hitQE\nXfBDHb8JE6B//3DzaxB6fvmwgm9is9NO8MEHsGKF70hMvubPhzVr3JfuJnzWwzexqqqCSy6BQYN8\nR2LyMXw4PPaYe5j0sh6+8WLAABg3zncUJl8N7RxTGqzgJyj0HmJNTU3QBT/E8Rs/Hg480D0PMb9s\noeeXj0gFX0SOEZF3RKReRPbaxH4DRWSmiMwWkUuiHNMUt332cWvqLFniOxKzOXPnwvLlsOeeviMx\nSYnUwxeRXYB1wF+Ai1T1rRz7lAGzgYOBhcAk4HhVzbnytvXw02/oUDjuODjxRN+RmE259VZ3odw9\n9/iOxESVSA9fVWep6hxgUwfqC8xR1XmquhZ4BBga5bimuA0eDE8+6TsKszmjRsERR/iOwiQpiR5+\nN2BB1vaHmddKTug9xIb8jjoKnnkGVq3yG0/cQhq/zz6DiRPdl+wNQsovl9Dzy0fLze0gImOBLtkv\nAQr8VlVH5XGMXGf/m+zZVFdXU1FRAUB5eTmVlZVfL23aMGi2XdzbfftW8fTTsO22xRGPbW+4vXhx\nFf37wxtvFEc8tt207YbntbW1NEUs8/BFZDzwm0Z6+P2AYao6MLN9KaCqen0jP8t6+AH43/+FsWNh\nxAjfkZhcDj8cjj8efvpT35GYOPiYh9/YwSYBPUWku4i0Bo4HrMMbuKOOgueeg5UrfUdiNrZoEbzy\nihsjU1qiTsv8sYgsAPoBo0XkmczrXUVkNICq1gPnAWOA6cAjqlqSt8rI/jgWouz8OnWCAw6ARx/1\nF0/cQhm/hx6Co4+GDh02fD2U/BoTen75iDpL5wlV3UFV26lqV1UdlHl9kaoekbXfs6q6i6p+R1Wv\nixq0SYdzzoHbb/cdhcmmCn/9K1RX+47E+GBr6ZiCqa93i3L94x/Qt6/vaAy4mTknnQRz5tgNT0Ji\na+kY71q0sLP8YnPzzXD22VbsS5UV/ASF3kPMld/Pf+4uwvr44+TjiVvax2/uXHj+eTjrrNzvpz2/\nzQk9v3xYwTcF1akTHHkk3Hmn70jMjTe6Yr/llr4jMb5YD98U3OzZsP/+rm+89da+oylNixdDnz4w\ncyZ07uw7GhM36+GborHzzm5BtRtv9B1J6br1VreYnRX70mYFP0Gh9xA3ld/vfgd33QVLlyYXT9zS\nOn4LF7orny+6aNP7pTW/fIWeXz6s4JtEdO/upgNee63vSErPZZfBGWdAZnkqU8Ksh28S09BHfust\nKz5JmTjRLaEwc6Z9WRsy6+GborPddq6tcM457opPU1jr1sEFF8Af/2jF3jhW8BMUeg8xn/wuvhj+\n/W/4+98LH0/c0jZ+w4e7op/viphpy6+pQs8vH1bwTaJatXK31LvwwjAuxipWS5fCf/6nm51TZv/K\nTYb18I0XDQX/wQd9RxIeVRgyBHbbzb4kLxX59vCt4BsvVq6EXXd16+wMHuw7mrDccQfcd59b8751\na9/RmCTYl7ZFKPQeYlPy69DBnd2fdhp88EHhYopTGsZvxgy48kr429+aXuzTkF8UoeeXDyv4xpvv\nfx8uv9xNG1y92nc06bdmjbvW4Q9/gF128R2NKUbW0jFeqcLJJ7svc++/35btbS5Vd1OTlSvh//7P\n/j+WGmvpmFQQgbvvdhdj3XWX72jSa9gwd3HVgw9asTeNi3pP22NE5B0RqReRvTaxX62ITBGRt0Xk\n9SjHTLPQe4jNza9DBxg50hWtMWNiDSlWxTp+998PDz8Mo0ZB+/bN/znFml9cQs8vH1HP8KcBRwIv\nbma/dUCVqu6pqnazO/MNPXu6on/SSWD/LvM3dixceik8/bSthGk2L5YevoiMB36jqm818v4HwD6q\n+kkeP8t6+CVs/Hg49lh44gk44ADf0RS3SZPg8MPhscfcF+CmdBVbD1+B50RkkoickdAxTQr98Ieu\nPXHkkfB6yTb/Nu+FF1yxv/deK/Ymf5st+CIyVkSmZj2mZf4c0oTj7K+q+wCDgXNFpH+zI06x0HuI\nceV32GHuwqEhQ9xZbLEolvF74gk4/ng3G2dIU/4Vbkax5FcooeeXj5ab20FVD4l6EFVdnPnzIxF5\nHOgLTGhs/+rqaioy6+eWl5dTWVlJVVUVsH7QbDvs7SOOqOKee2DAgBouuAB+//viis/X9iWX1HDP\nPTB2bBV77eU/Htv2s93wvLa2lqaIs4d/kaq+meO99kCZqn4hIh2AMcBVqppzPob18E22KVPc7RFP\nOcXN4ikr0YnEqvDf/+2WTRgzxt020pgGifTwReTHIrIA6AeMFpFnMq93FZHRmd26ABNE5G3gNWBU\nY8XemI3tsYfr5b/wAvzkJ+7ColLzySful96jj8LLL1uxN80XqeCr6hOquoOqtlPVrqo6KPP6IlU9\nIvP8A1WtzEzJ3E1Vr4sj8DTK/jgWokLl17kzPP+8u4lH//7uAiMffIzfyy/Dnnu6pRImTIAddijc\nsezvZ/hK9AOySZs2bdwFRmec4Yr+tdfC2rW+oyqc+nq45ho3RfWuu1w7p7WtfGkisrV0TOrU1sKZ\nZ7r19O+7DyorfUcUr8mT4fzzoUULN0W1WzffEZliV2zz8I2JTUUFPPecK4qHHgpXXBHGaptLlrhP\nMAMHwoknwrhxVuxNvKzgJyj0HmKS+Ym41SGnTIHZs6FHD7j5Zli1qnDHLFR+a9bADTfAd78LHTu6\n7yjOPtud4SfJ/n6Gzwq+SbWuXWHECLeWzMsvu8J/003pmM2zYoW741efPvCvf8Grr8KNN0J5ue/I\nTKish2+CMmUKXH21m9FywQVu/n6xtUVmzYI//9ndlergg11rypZHMFFYD9+UpD32cPPVx46F995z\nN/I+6CC35szy5f7i+uILtxrooYfCgQe61s3UqW55BCv2JilW8BMUeg+xmPLbbTdX5BcuhHPPdS2f\n7t3d7RQJWVmYAAAE00lEQVQfftjdR7epHySbkp8qvPuuay8dfLBrPd1xB/z0pzBvnptyuf32TTt+\noRXT+BVC6PnlY7Nr6RiTZm3bwtFHu8fy5e4se+RIuPhiV5T79YP99nN/7rWXu7irqerr3S+QmTPd\n49133ZXBqjBokGstHXQQbLFF/PkZ0xTWwzclSRXmz4fXXnOPV191LZYWLaBLF9huu/WPTp3cRV5f\nfrnhY+VKeP991zrq0gV69YLevd2f3/++e263GzRJyLeHbwXfmAxVN3Nm8eINH8uWuatc27Zd/2jT\nBtq1g512csseRLm1oDFRWcEvQjU1NV8vcxoiyy/dLL/0slk6xhhjNmBn+MYYk3J2hm+MMWYDVvAT\nFPo8YMsv3Sy/8FnBN8aYEmE9fGOMSTnr4RtjjNlA1JuY3yAiM0Rksog8JiJbNbLfQBGZKSKzReSS\nKMdMs9B7iJZfull+4Yt6hj8G+K6qVgJzgMs23kFEyoA/A4cB3wVOEJFeEY+bSpMnT/YdQkFZfulm\n+YUvUsFX1XGqui6z+RqQa/2/vsAcVZ2nqmuBR4ChUY6bVst9rs+bAMsv3Sy/8MXZw/8Z8EyO17sB\nC7K2P8y8ZowxJkGbXR5ZRMYCXbJfAhT4raqOyuzzW2Ctqg7P9SNyvFaS03Bqa2t9h1BQll+6WX7h\nizwtU0ROBc4EDlLVNTne7wcMU9WBme1LAVXV6xv5eSX5y8AYY6LIZ1pmpBugiMhA4D+BA3MV+4xJ\nQE8R6Q4sAo4HTmjsZ+YTtDHGmKaL2sO/DdgCGCsib4nIHQAi0lVERgOoaj1wHm5Gz3TgEVWdEfG4\nxhhjmqjorrQ1xhhTGEV7pa2IXCQi60RkG9+xxElEfi8iU0TkbRF5VkS28x1TnPK9GC+tROQYEXlH\nROpFZC/f8cQh9AsjReReEVkiIlN9xxI3EdleRF4QkXdFZJqInL+p/Yuy4IvI9sAAYJ7vWArgBlXd\nQ1X3BJ4CrvQdUMw2ezFeyk0DjgRe9B1IHErkwsj7cfmFqA64UFX7APsB525q/Iqy4AM3Axf7DqIQ\nVPWLrM0OwLrG9k2jPC/GSy1VnaWqc8g93TiNgr8wUlUnAJ/6jqMQVHWxqk7OPP8CmMEmrnOKNEun\nEERkCLBAVaeJhPJvakMicg1wCrAc+KHncArpZ7gCYopXrgsj+3qKxUQgIhVAJTCxsX28FPxNXMx1\nBXA5cMhG76XK5i5WU9UrgCsy/dJfAsOSj7L5YrgYr6jlk19A7MLIAIjIFsCjwAUbdRE24KXgq+oh\nuV4XkV2BCmCKuNP77YE3RaSvqi5NMMRIGssvh7/j+vjDChdN/DaXX+ZivMHAQclEFK8mjF8IPgR2\nzNreHljoKRbTDCLSElfsH1LVf25q36Jq6ajqO8DXs1ZE5ANgL1UNpv8mIj1V9b3M5lBczy0YeV6M\nF4rUffrMoUkXRqaYEMZ45XIf8K6q3rq5HYv1S9sGSniDdJ2ITBWRybiZSBf4DihmOS/GC4WI/FhE\nFgD9gNEikmvBwNQohQsjRWQ48Aqws4jMF5HTfMcUFxE5ADgJOCgz1futzElX7v3twitjjCkNxX6G\nb4wxJiZW8I0xpkRYwTfGmBJhBd8YY0qEFXxjjCkRVvCNMaZEWME3xpgSYQXfGGNKxP8DxGqSfXtB\nnI4AAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"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)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false,
"jupyter": {
"outputs_hidden": false
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAW8AAAEACAYAAAB8nvebAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYFdWZx/HvAQRFxR5UkBGhRzGjAoqKGgFNQzAorvFR\nE+NEOpHEjEkMD7jjRBITE02IuBE1EVFcyBgziriBiYUKKgK2IAFEDYsLiwooKkt3n/nj0KSD9/Z2\n69apOvf3eZ5+2upu6r6+nPtS/dapc4y1FhERyZZWvgMQEZHmU/EWEckgFW8RkQxS8RYRySAVbxGR\nDFLxFhHJoDaFnsAY0w54Dmi77Xx/ttb+rNDziohIfiaOed7GmPbW2s+MMa2BmcDF1trZBZ9YRERy\niqVtYq39bNt/tsNdfevJHxGRIoqleBtjWhljXgVWAdOtta/EcV4REcktrivvWmvt4UBX4BhjzCFx\nnFdERHIr+IZlfdbaj40xEXAi8Pf63zPGqJUiItIC1lqz49fimG2yF7DVWrvBGLMLMBj4dZ4ACn25\nTBszZgxjxozxHUbqhZwna2HuXFiwANasgbVrv/gZ4LjjYOBA97H//mC+8NZ1Qs5VnLKcJ5PnLz+O\nK+8uwD3GmFa4NsyfrLVPxHDe4Cxbtsx3CJkQYp42boQHHoDbb4cNG1xx7tQJOneG3r1h773d8d57\nw5YtMGMGPPss/M//wE47waBBrpCfeKL7M3VCzFUxhJingou3tXYBcEQMsYgEZ/58V7AnT4avfAV+\n/WsYPBhaNXK36cADYfhwd6W+ZIkr5FOnwiWXuHN897v5r8alNMQyz7tJL2SMLfW2SRRFVFRU+A4j\n9bKep9paePBBuO02WLECvvc9uOAC6Nq18HMvWADnn+/O9Yc/wOLF2c5VUrI8powxOXveKt4iMVq9\nGoYNg48+gquuglNOgTaxTgtwbZVrr3XF+5Zb4Oyz4z2/pEu+4q21TRIURZHvEDIhq3l6+mk4/HDo\n2xdmzoQzzoi/cAO0beuK96OPwsiREeedB+vWxf86IcnqmGqIirdIgbZsgUsvda2R+++HX/zC3WQs\ntmOOcVffHTtCv37w2WeN/xkJh9omIgV4800491zo0gUmTIC99vITx3/9F5SVwa23+nl9KR61TURi\ndv/9cOyx7gbio4/6K9zgivaUKTBtmr8YJFkq3gkKse9WDGnPk7WuNTJmDDzzDPz4x/6m7dXlqqwM\n7r7btW4++shPLGmW9jHVEireIs1gLVxxBfzpT/D883DYYb4j+qevfhXOPBN++EPfkUgS1PMWaaLa\nWrj4Ynj5ZXjqKdhzT98RfdHnn8MRR8A118A3v+k7GomD5nmLFKC62j1s8+ab8Pjj0KGD74jymzMH\nhg6FV1+Ffff1HY0USjcsUyDEvlsxpC1PW7bAt74F777rrrjTVLhz5apvXxfvHXckH09apW1MxUHF\nW6QBmza5PvLmzW42x667+o6oaYYMgRde8B2FFJPaJiJ5bNwIp5/uVvu7995kHryJy/r1sN9+buZJ\nluKWL1LbRKQZ1q+Hr30N/uM/4L77slcAy8pc7FVVviORYlHxTlCIfbdi8J2ntWvd+tlHHQV33gmt\nW3sNp0EN5WrAALVO6vgeU8Wg4i1Sz3vvuXW3hw6FceMaX3c7zVS8w6aet8g2y5a5jRIuuACuvNJ3\nNIVbvtwtXvX++9q4IcvU8xZpwLx5cPzx8JOfhFG4Abp1c736t97yHYkUg4p3gkLsuxVD0nl68EE3\ntW7sWLdOSZY0lCtj1DqpE+J7T8VbSlZNDVx+OYwe7RaYCnFHmgED3MYQEh71vKUkrVvnnkLcvBn+\n93/9LudaTE88ATfd5Hb5kWxSz1tkm7lz3Y28L33JFbVQCzfAypXuYR0Jj4p3gkLsuxVDsfK0datb\ng/ukk9znm27K3sM3O2osV8uXQ/fuycSSZiG+94qwPapI+ixc6Ha86dSptFbbW7HCPSkq4VHPW4JW\nXQ2/+x385jdw3XUwfHhpzXk+7ji303xFhe9IpKXy9bx15S3Bev55t6tMp04we7Zb66PUqG0SLvW8\nExRi360YCs3T6tUwbJibTXL11TB9eriFu6FcVVfDqlXQtWty8aRViO89FW8JxtatcPPN0KsXdO4M\nixbBOeeUVpukvnffdXnI+k1Zya3gnrcxpitwL7APUAP8wVp7c46fU89biubJJ2HkSDctbtw4OOQQ\n3xH599xz7lF/PaSTbcXseVcDI621VcaY3YC5xphp1trFMZxbpEGLFsGoUW79jrFj4eSTS/dKe0fq\nd4et4LaJtXaVtbZq239vBBYBJTIRq3lC7LsVQ1Py9NFHbhGp4493U+EWLIBTTim9wt1QrlasUPGu\nE+J7L9aetzGmHOgDvBzneUXqVFfD+PFw8MGux71oEYwYAW3b+o4sfXTlHbbY5nlva5lEwLXW2kdz\nfF89bynI3/7mrrb33tv1tQ891HdE6TZkiPuH7aSTfEcihSjqPG9jTBvgz8CkXIW7TmVlJeXl5QCU\nlZXRp08fKrY9PVD3a42Odbzj8apVcN55EQsXwvjxFXz96zBjRkQUpSO+tB4vWgTduqUnHh037TiK\nIiZOnAiwvV7mEsuVtzHmXuADa+3IBn6m5K+8oyja/pcl+dXlqaYGbr/drUMyfLibs73rrr6jS5d8\nY8pal6s1a2C33ZKPK22y/N4r2pW3MaY/cB6wwBjzKmCBq6y1TxV6bild8+e77ch22QWiCHr29B1R\ntixcCF26qHCHTGubSKpUV8MNN8CNN8Kvfw3f/W7pzSCJwy23wGuvwR//6DsSKZTWNpHUW7zYPdbe\noYNbc7tbN98RZVcUwZln+o5CikmPxyeo7qaE/Ctr4fe/d1t2DRsGV14ZqXA3Ua4xVVsLM2ZoJcH6\nQnzv6cpbvNq4ES68EF5/HV58EQ480F01SsstWAAdO5bOmuWlSj1v8WbRIjjrLLcl2W23uZuTUrhx\n41xu77jDdyQSB+1hKanyyCPu0fZRo2DCBBXuOD37LAwc6DsKKTYV7wSF2Hdridtuc5skPPmkm02y\nI+Wp6XbMVU2NW01Q/e5/FeKYUs9bElNbC1dd5a66X3gh3A0SfKqqcvO799nHdyRSbOp5SyK2boXK\nSli2DKZMgT339B1RmMaOhbffdr/dSBjU8xZvqqvhvPNgwwZ45hkV7mJSv7t0qHgnKMS+W2NqatwV\n98cfw8MPN+3GZCnmqaXq56q62rWj1O/+ohDHlHreUjS1tW4O93vvwdSp0K6d74jCNm+eeyp1r718\nRyJJUM9bimbkSJg9G556SgskJeH6692mwzd/YQdZyTL1vCVRd93lrrYfe0yFOymzZrm581IaVLwT\nFGLfLZeZM92u5VOmwL/9W/P/fKnkKQ71c/X669C7t79Y0izEMaXiLbFasQLOPhvuuQcOOsh3NKXj\n00/h/ffhgAN8RyJJUc9bYrN1K/TvD+ecA5dc4jua0jJnjtu84rXXfEcicVPPW4ruuuvcanajRvmO\npPQsXKjdhkqNineCQuy71ZkzB8aPdzcqC935JuQ8xa0uVyreDQtxTKl4S8E+/xy+/W246SatIe2L\ninfpUc9bCnbFFW7NksmTfUdSusrLYfp0t5mFhCVfz1vFWwry1ltw9NFumlqXLr6jKU0bN0KnTvDJ\nJ9C6te9oJG66YZkCIfbdLrvM3aCMs3CHmKdiiaKIpUuhRw8V7oaEOKa0tom0WBS59TTuv993JKVt\n1Sr91lOK1DaRFrEWjjzS9bvPOcd3NKVt4kS3FOw99/iORIpBbROJ1ZNPulUDzz7bdySyahV07uw7\nCkmaineCQuq7XX89XH554XO6cwkpT8UWRRGrV6t4NybEMaXiLc02axasXKmr7rRYtUp7VpYi9byl\n2U4/HYYMgYsu8h2JAAwa5DZ2HjzYdyRSDPl63pptIs2yYoXbauvBB31HInXWrWvZ0ruSbbG0TYwx\ndxljVhtj5sdxvlCF0HebNAm+8Q1o3754rxFCnpISRRGffAIdOviOJN1CHFNx9bzvBobEdC5JKWvd\ntLTKSt+RSH2ffAK77+47CklabD1vY0x34DFr7aF5vq+ed8bNnAnf+55bBKkYs0ykZdq3h7VrYddd\nfUcixaB53lKwSZPg/PNVuNOkuho2by5uG0vSKdEblpWVlZSXlwNQVlZGnz59qKioAP7Zkwr5uKqq\nihEjRqQmnuYcP/tsxEMPwUsvFf/16vcn0/L/n9bjTz+F9u0rmDEjHfGk9XjcuHGZqTdRFDFx4kSA\n7fUyF7VNEhRF0fa/rKyZNw/OPReWLCn+a2U5T0l75JGI4cMr+OAD35GkW5bHVNGXhDXGlOOKd879\nq1W8s+3nP4cNG2DsWN+RSH3vvAPHHAPvvus7EimWova8jTEPALOALxljVhhjvhPHeSU9pk6Fk0/2\nHYXsaPNmaNfOdxTiQyzF21r7LWvtv1tr21lru1lr747jvKGp38vNkvXrYdEiGDAgmdfLap58eP75\nSMW7CUIcU5ptIo166SXo2xfatvUdieyouhp22sl3FOKD1jaRRv30p1BTA7/8pe9IZEdz57q59/Pm\n+Y5EikXzvKXFZs6Efv18RyG5VFdDG61QVJJUvBOUxb5bTQ3Mng3HHpvca2YxT7688kqkvSubIMQx\npeItDfrHP2DPPaFjR9+RSC61tbryLlUq3gnK4kMCr78OPXsm+5pZzJMvhx1WoSvvJghxTKl4S4MW\nLoRevXxHIfnU1EArvYtLkv7aE5TFvtvrrydfvLOYJ19efTVS8W6CEMeU/tqlQUuWwEEH+Y5C8rFW\nqzyWKs3zlgZ17AhvvAF77eU7Esll2jT47W/dZwmT5nlLs338MWzZ4mabSDrV1urKu1SpeCcoa323\n5cuhe/fki0PW8uTT/PmRincThDimVLwlr7riLemlTmTpUvFOUNbmmq5aBV26JP+6WcuTT717V+jK\nuwlCHFMq3pLXmjXQqZPvKKQxKt6lScU7QVnru/kq3lnLk0/z50e+Q8iEEMeUirfkpSvv9FPPu3Sp\neCcoa323Dz/0M00wa3ny6dBD1fNuihDHlIq35PXxx7DHHr6jEJFcVLwTlLW+24YN0KFD8q+btTz5\npJ5304Q4plS8JS9deYukl9Y2kbw6dICVK1XA0+zxx2H8ePdZwqS1TaTZPv8c2rf3HYWI5KLinaAs\n9d2qq/1tsZWlPPmmnnfThDimVLwlp88/h1120dN7ImmlnrfktHYtHHwwfPCB70ikIep5h089b2mW\nrVuhbVvfUYhIPireCcpS36262k+/G7KVJ9/U826aEMdULMXbGHOiMWaxMeYNY8zlcZxT/PJZvEWk\ncQUXb2NMK+BWYAjQEzjXGKMta3PI0voKPot3lvLk26GHVvgOIRNCHFNxXHkfDSy11i631m4FJgOn\nx3Be8ai2FlqpqSaSWnG8PfcFVtY7fmfb12QHIfbdikF5aroZMyJeftl3FOkX4piK4xfjXDOBc84J\nrKyspLy8HICysjL69Omz/deZuuSGfFxVVZWqeBo6nj074rPPANIRj45zH69Z45buTUs8aT2uqqpK\nVTwNHUdRxMSJEwG218tcCp7nbYz5MjDGWnvituMrAGutvX6Hn9M87wxZvBjOOMN9lvTSPO/wFXOe\n9ytAD2NMd2NMW+CbwJQYziue6d9akfQquHhba2uAHwHTgIXAZGvtokLPG6K6X42yoE0bqKnx89pZ\nypNvmufdNCGOqVgmg1lrnwL+M45zSTq0aeOmC4pIOmkyWILqbk5kgc/inaU8+aZ53k0T4phS8Zac\n2rRx65uISDqpeCcoS323nXeGzZv9vHaW8uSbet5NE+KYUvGWnHbZxa3pLSLppPW8JSdroXVr1zpp\n3dp3NJKP5nmHT+t5S7MY41onmzb5jkREclHxTlDW+m677Qaffpr862YtTz6p5900IY4pFW/Ja489\nYMMG31GISC4q3gnK2lzTDh3g44+Tf92s5cknzfNumhDHlIq35NWhg668RdJKxTtBWeu7lZXBunXJ\nv27W8uTT/PmRFhBrghDHlIq35NWpE6xd6zsKaYjJtZq+lAQV7wRlre/WqZNb7D9pWcuTT+p5N02I\nY0rFW/LyVbyledQ2KU0q3gnKWt+tc2dYvTr5181annxasCDyHUImhDimVLwlr65dYeXKxn9O/NKV\nd2nS2iaS13vvwRFHwKpVviORfJ56Cm68EZ5+2nckUixa20SabZ99YP16rS6YZsboyrtUqXgnKGt9\nt1atYL/9YMWKZF83a3nySfO8mybEMaXiLQ064ABYutR3FJJPq1a68i5VKt4JyuJc0549YeHCZF8z\ni3ny5fDDK1S8myDEMaXiLQ3q1Qtef913FJJPq1ZQW+s7CvFBxTtBWey79eyZfPHOYp58ee21iJoa\n31GkX4hjSsVbGnTIIbBkiXaST6tWrVDxLlGa5y2N6tUL7rkHjjzSdySyo5deghEj3GcJk+Z5S4v1\n6wczZ/qOQnJp0waqq31HIT6oeCcoq323/v1h1qzkXi+refKhqipS8W6CEMeUirc0ql+/ZIu3NF3r\n1rofUaoK6nkbY84CxgAHA0dZa+c18LPqeWeUtbDvvvDcc9Cjh+9opL433oCTT9aDVCErVs97AfB1\nYEaB55EUMwaGDoXHH/cdieyoXTvYvNl3FOJDQcXbWrvEWrsU0GZMTZDlvtspp8DUqcm8VpbzlLS5\ncyMV7yYIcUyp5y1NMniwm472ySe+I5H6dtpJV96lqk1jP2CMmQ50rv8lwAKjrbWPNefFKisrKS8v\nB6CsrIw+ffpsX3Og7l/G0I/rpCWeph7PmRNx0EHw9NMVnHVWcV+voqLC+/9vVo5POKGCzz9PTzxp\nPa77Wlriaeg4iiImTpwIsL1e5hLLQzrGmGeBUbphGbYJE2DKFHjkEd+RSB1roW1b+PRT91nCk8RD\nOup7N6LuX9esOvtsiKLib0qc9TwlacaMiN13VzurMSGOqYKKtzHmDGPMSuDLwFRjzJPxhCVptPvu\ncNpp8MADviOR+lS8S5PWNpFm+dvfYORIqKryHYnU6dkTJk+G3r19RyLFoLVNJBYVFe4qT2udpMce\ne8CGDb6jkKSpeCcohL5bq1ZwySVw/fXFe40Q8pSUKIro1Kn49yGyLsQxpeItzVZZCbNnJ789muTW\nuTOsXu07Ckmaet7SItdd5zZpuOce35HINde4zz/7md84pDjU85ZYXXSRe1z+rbd8RyK68i5NKt4J\nCqnvVlYGo0bBZZfFf+6Q8lRsURSxzz4q3o0JcUypeEuLjRwJ8+a5B3fEn86dYdUq31FI0tTzloI8\n9JDrf8+Z4zYGkOS9847bX1RX32FSz1uK4qyzYLfd4K67fEdSuvbd160s+MEHviORJKl4JyjEvpsx\nMH48jB4Ny5bFc84Q81QsURRhDBxyiKZuNiTEMaXiLQXr3RsuvdTN/66t9R1NaerZU8W71KjnLbGo\nqXGPzn/96+5GpiRr3Dh480249VbfkUjc1POWomrd2j2w86tfwfz5vqMpPbryLj0q3gkKse9W3/77\nuyvAM84o7OZZ6HmKU12uVLwbFuKYUvGWWJ13ntu04eyzYetW39GUji5dXL7XrvUdiSRFPW+JXU0N\nnH46dOvmZqJIMgYMgF/8wt17kHCo5y2Jad3a7bYTRa6NIsno2xdeesl3FJIUFe8Ehdh3y6dDB3ji\nCbjxRrjzzub92VLKU6Hq56qiAp591lsoqRbimGrjOwAJV3k5/PWvMHAgtGsHw4b5jihsX/kKnH++\n633vtJPvaKTY1POWolu8GAYNgt/9Dr75Td/RhO2II9xc7379fEcicVHPW7w56CCYNs09vHP77b6j\nCdvAgWqdlAoV7wSF2Hdrql694Pnn3dX36NHQ0C9hpZyn5toxVyreuYU4plS8JTEHHOB2nf/rX13/\ne8sW3xGF57jj4OWX3SqDEjb1vCVxn30G554LGzbA5Mmwzz6+IwrLUUfB2LFw/PG+I5E4qOctqdG+\nPfzlL6649O3r2ikSH7VOSoOKd4JC7Lu1VOvW8POfwx/+4B6l/+1v/9kHV56aLleuBg7U1nQ7CnFM\nqXiLVyedBLNnu+3UTjtNezHGYcAAty3dpk2+I5FiUvFOUIUWncipWzfXOjnsMPexenWF75AyI9eY\n2n13N7vnxReTjyetQnzvFVS8jTE3GGMWGWOqjDEPG2M6xBWYlJa2bd2iSo89BmPGwDe+oT0ZC6FH\n5cNX6JX3NKCntbYPsBS4svCQwhVi3y1uRx8N48ZFdO3qtlebNKnhOeGlLt+Y0k3LfxXie6+g4m2t\nfcZaW7dr4UtA18JDklLXrp2b6vboo25hq0GD3CP20nT9+7u+d3W170ikWGKb522MmQJMttY+kOf7\nmuctzVZd7dYEv/Za+P734aqrYNddfUeVDV27wqxZ7p6CZFeL53kbY6YbY+bX+1iw7fOp9X5mNLA1\nX+EWaak2beDii+G11+Dtt+Hgg+HBB9VKaYru3WH5ct9RSLE0uiSstfaEhr5vjBkGDAUGNXauyspK\nysvLASgrK6NPnz7b7wLX9aRCPq6qqmLEiBGpiSetx/+6RrX7/htvRFx4IVx0UQU/+Qlcd13Ej34E\nF17oP16fx3Vfy/X9du1g+fIKjjsuPfH6Oh43blxm6k0URUycOBFge73MpaC2iTHmRGAscLy19sNG\nfrbk2yZRFG3/y5L8GstTTQ1MnAhXXw1Dh8Ivf1m6j9g3lKsrrnDTBkePTjamNMryey9f26TQ4r0U\naAvUFe6XrLUX5fnZki/eEq8NG1zhnjABLrkERoyAnXf2HVV6/P73UFUFd9zhOxIpRFHWNrHWHmit\n7W6tPWLbR87CLVIMe+wBN9zg9m2cPdv1wx96SP3wOup5h01PWCaofp9S8mtunnr0cAtdTZjgHvT5\n6lfh738vTmxp01CuunVT8a4T4ntPxVuCMXAgzJ0LZ57p9nO89FL45BPfUfnTvTusWKHfREKl9bwl\nSKtXw+WXu40fxo51KxeaL3QNw9exI7zxBuy1l+9IpKW0nreUlM6d3YyUBx90S8+edRasWeM7quSp\n7x0uFe8Ehdh3K4Y48zRggGul9OjhViz8v/+L7dSp0Fiu1Pd2QnzvqXhL8Nq1g+uvh4cfhssug29/\nG9at8x1VMnTlHS71vKWkfPopXHmlW/Tq4YfdNmwhGzsWVq6EceN8RyItpZ63CG5Rq5tvhptucrv4\nTJrkO6LiqptxIuFR8U5QiH23YkgiT2ec4da7HjMGRo3K7tKpjeXKWrecQKkL8b2n4i0lq1cveOUV\nmD/frZHy0Ue+I4rfrFlw7LG+o5BiUM9bSl51tbuROWUKPP00HHCA74jic9RRbkOLAQN8RyItVZSF\nqZoZgIq3pNr48XDddTB9ulsnJes2bnTz3T/8UAt2ZZluWKZAiH23YvCVp4sucqsUDhrkVuPLgoZy\nNXs29Omjwg1hvvca3YxBpJQMG+ZmpAwZ4qYTfvnLviNquRdeULskZGqbiOTwxBOukD/0EGR0DX++\n9jX48Y/h1FMb/1lJL/W8RZopiuCcc9waKUOH+o6meaqr3aJU//gH7Lmn72ikEOp5p0CIfbdiSEue\nKirgscfgO99xT2OmUb5cLVjgdo9X4XbSMqbipJ63SAOOOcZNHzzpJPdo/fnn+46oadTvDp/aJiJN\nsHgxnHCC28z3Bz/wHU3DNm50s0xuvRVOPNF3NFIo9bxFCvT22zB4sJtSeMklvqPJ78ILYcsWuPtu\n35FIHNTzToEQ+27FkNY87b8/PPcc/PGPbk2UNFyL7JirqVNdm+emm/zEk1ZpHVOFUM9bpBm6dnUF\n/IQT3P6YN9wArVv7jspZuxa+/323e1CHDr6jkWJT20SkBT76yG2tZi3cdx/su6/feKx18ey/P/zm\nN35jkXipbSISo44d3RoogwfDkUe6Ra18mjTJbTR87bV+45DkqHgnKMS+WzFkJU+tW7vZJ3/5C1x8\nsXuacdOmZGOIooh589wN1Pvu0zom+WRlTDWHirdIgfr1cwtZrV7t5oUvWpTM61ZXw733uumA48e7\nDZaldKjnLRITa+Guu9wemddeCxdcADvtVJzXWrLEPTC0xx4wYYK7kSphUs9bpMiMgeHD3WyUyZOh\nvBx++lO3AXBcamvdHpz9+0NlpZsWqMJdmgoq3saYnxtjXjPGvGqMecoYs09cgYUoxL5bMWQ9Twcf\n7Ba1mjYN1q1zTzuefjo8+WTL95P87DN45hk3RXHyZHjxRfjv/4YZM6I4Qw9W1sdULoVeed9grT3M\nWns48DhwTQwxBasqKyv8exZKnnr2hFtucbu3n3oqXH019OgBv/oVzJwJS5fChg25H/bZtMn9A3DN\nNXD88dCpk3sw6LTT4Pnn4cAD3c+FkqtiCzFPBT2kY63dWO9wV6C2sHDCtn79et8hZEJoedp1V9dO\nGT4c5syBO+90UwvXrHEP1mzeDHvv7T46dXKPts+ZA4ccAgMHuqLfv787z45Cy1WxhJingp+wNMb8\nAjgfWA8MLDgikYD17es+6tu0yRXxumIObgaLnpKUhjRavI0x04HO9b8EWGC0tfYxa+3VwNXGmMuB\nHwNjihFoCJYtW+Y7hEwotTztvDPst5/7aK5Sy1VLhZin2KYKGmO6AY9ba3vn+b7mCYqItECuqYIF\ntU2MMT2stW9uOzwdyPt4Qq4XFxGRlinoytsY82fgS7gblcuBH1hr348pNhERySOxJyxFRCQ+esLS\nE2PMJcaYWmNMR9+xpJEx5gZjzCJjTJUx5mFjjOZe1GOMOdEYs9gY88a2yQKSgzGmqzHmb8aYvxtj\nFhhjLvYdU1xUvD0wxnQFBuNaTZLbNKCntbYPsBS40nM8qWGMaQXcCgwBegLnGmMO8htValUDI621\nhwDHAj8MJVcq3n7cCFzqO4g0s9Y+Y62te+jrJUArePzT0cBSa+1ya+1WYDJuwoDswFq7ylpbte2/\nN+ImVXjeOiMeKt4JM8acCqy01i7wHUuGfBd40ncQKbIvUH+5q3cIpCAVkzGmHOgDvOw3knhoD8si\naODBpqvNJ8zmAAABIElEQVSBq4ATdvheSWrsAbBtPzMa2GqtfcBDiGmVa8xo5kEDjDG7AX8GfrLD\nsh6ZpeJdBNbaE3J93RjTCygHXjPGGFwrYK4x5mhr7ZoEQ0yFfHmqY4wZBgwFBiUTUWa8A3Srd9wV\neM9TLKlnjGmDK9yTrLWP+o4nLpoq6JEx5h/AEdbadb5jSRtjzInAWOB4a+2HvuNJE2NMa2AJ8FXg\nfWA2cK61NqE9fLLFGHMv8IG1dqTvWOKknrdflhJumzTiFmA3YLoxZp4xZrzvgNLCWlsD/Ag3I2ch\nMFmFOzdjTH/gPGDQtn0H5m27MMg8XXmLiGSQrrxFRDJIxVtEJINUvEVEMkjFW0Qkg1S8RUQySMVb\nRCSDVLxFRDJIxVtEJIP+H0uU6Ejf+wSoAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"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)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false,
"jupyter": {
"outputs_hidden": false
}
},
"outputs": [],
"source": []
}
],
"metadata": {
"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.8.3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}