ketch · March 9, 2017 15:31
diff --git a/BDF2_simple.ipynb b/BDF2_simple.ipynb
 {
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simple implementation of BDF2 for a linear, scalar ODE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the BDF2 formula is\n",
    "\n",
    "$$y^{n+1} = \\frac{4}{3} y^n - \\frac{1}{3} y^{n-1} + \\frac{2}{3} h f(y^{n+1})$$\n",
    "\n",
    "For a scalar ODE with $f(y) = \\alpha y$, this means\n",
    "\n",
    "$$y^{n+1} = \\frac{4}{3} y^n - \\frac{1}{3} y^{n-1} + \\frac{2}{3} h \\alpha y^{n+1},$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "so\n",
    "\n",
    "$$y^{n+1} = \\left(1-\\frac{2}{3}h\\alpha\\right)^{-1}\\left(\\frac{4}{3} y^n - \\frac{1}{3} y^{n-1}\\right).$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def bdf2(h, T, alpha, y0):\n",
    "    # Variable meanings:\n",
    "    # y_nm1, y_n, y_np1: y_{n-1}, y_n, y_{n+1}\n",
    "    # h: Time step\n",
    "    # T: Final time\n",
    "    # y0: Initial value\n",
    "    \n",
    "    t = 0\n",
    "    times = [0]\n",
    "    y_n = y0\n",
    "    \n",
    "    # Euler step to start\n",
    "    y_np1 = y_n + h*alpha*y_n\n",
    "    \n",
    "    t += h\n",
    "    times.append(t)\n",
    "    y_nm1, y_n = y_n, y_np1\n",
    "    \n",
    "    while t<T:\n",
    "        if T-t < h:  # Don't step past the end time\n",
    "            h = T-t\n",
    "        \n",
    "        y_np1 = ( (4./3)*y_n - (1./3)*y_nm1 ) / (1.-(2./3)*h*alpha)\n",
    "        t += h\n",
    "        times.append(t)\n",
    "        y_nm1, y_n = y_n, y_np1\n",
    "\n",
    "    exact_solution = np.exp(T*alpha)*y0\n",
    "    error = np.abs(y_n - exact_solution)\n",
    "    return y_n, error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "y0 = 1.\n",
    "T = 1.\n",
    "\n",
    "alpha = 1.j\n",
    "y_values = []\n",
    "err_values = []\n",
    "h_values = (0.5, 0.25, 0.125, 1./16, 1./32,)\n",
    "\n",
    "for h in h_values:\n",
    "    y, error = bdf2(h, T, alpha, y0)\n",
    "    y_values.append(y)\n",
    "    err_values.append(error)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x110ee2cd0>"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEACAYAAAC6d6FnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlYVeXax/HvAhQcUXEWEBSVQRxRU7OcZzK1jmmZQ2+d\n6pSZlTlmZqllg1YezdQ8lZWncsIJ55xzVmZUHAAHBGSSQYbn/WOZx3SDyLAH9v25rq5Lts9a+4Fk\n/fZa93rupSmlEEIIYX1sTD0BIYQQpiEBIIQQVkoCQAghrJQEgBBCWCkJACGEsFISAEIIYaUkAIQQ\nwkpJAAghhJWSABBCCCslASCEEFbKztQTKEjNmjWVm5ubqachhBAW5dixY/FKqVoPGmfWAeDm5sbR\no0dNPQ0hhLAomqZdLMw4uQQkhBBWyiwDQNM0f03TliQnJ5t6KkIIUWaZZQAopQKUUi85OjqaeipC\nCFFmmXUNwJDs7GxiYmLIzMw09VTKBAcHB5ydnSlXrpyppyKEMDKLC4CYmBiqVKmCm5sbmqaZejoW\nTSlFQkICMTExuLu7m3o6QggjM8tLQAXJzMzEyclJDv4lQNM0nJyc5GxKCCtlcQEAyMG/BMnPUggz\nk50Jp1aBER7Xa3GXgMyBra0tvr6+d75+5plnmDRpkglnJIQoE6L+gA1vQuI5cPIA57al+nYSAEVQ\noUIFTp48WeCY3NxcbG1t73ydk5ODnd2Df9yFHSeEKENuJsDWqXDqZ6juDiPXlvrBHyQASpSbmxvD\nhg1j27ZtTJw4kcWLF9OqVSv27dvH8OHDGTp0KGPHjiU+Pp5atWrx3Xff4erqyujRo3FwcODEiRN0\n7tyZQYMG8cYbbwD6JZo9e/ZQpUoVE393QogSpxSc/Am2ToOsFOjyNjz2NpSrYJS3t+gAmBkQQujl\nlBLdp3f9qszw9ylwTEZGBq1atbrz9eTJkxk2bBgATk5OHD9+HIDFixdz69atO+0s/P39GTVqFKNG\njWL58uWMGzeOtWvXAvrdTQcOHMDW1hZ/f38WLlxI586dSUtLw8HBoUS/RyGEGYg/o1/uubAXXDqA\n/wKo7WXUKVh0AJhKQZeA/goCQ18fPHiQ1atXAzBy5EgmTpx45++efvrpO5eMOnfuzIQJE3j22WcZ\nMmQIzs7OJf0tCCFMJScL9s2HvZ/qn/QHzoc2o8DG+PfkWHQAPOiTuilUqlSpwK8Ls92kSZMYMGAA\nmzZtonPnzgQGBuLp6Vmi8xRCmMCF/bBhPMRHQvOh0GcOVKljsulY5G2glqpTp0788ssvAKxcuZIu\nXboYHHfu3Dl8fX159913adeuHeHh4cacphCipKUnwrrXYEV//Qzg2d/hqeUmPfiDhZ8BmMq9NYC+\nffsyd+7cB2731VdfMWbMGObNm3enCGzI/Pnz2bVrFzY2Nvj4+NCvX78Sm7sQwoiUgtP/hcApkHED\nOr8Bj0+C8hVNPTMANGWExQYPS9M0f8Dfw8PjxTNnzvzt78LCwvDyMm6hpKyTn6kQpSDhHGycAFG7\noYGfXuSt29wob61p2jGllN+DxpnlJSDpBiqEsFg5t2DPp7CoE8Qeh/6fwgtbjXbwfxhyCUgIIUrK\npUMQMB6uh4H3IOj7MVStZ+pZ5UsCQAghiivjBmx/H46tAEcXGL4KmvU19aweSAJACCGKSikIWQ2b\nJ0F6PHR8DbpOBvvKxdptdGI6LjVKv1AsASCEEEVx4wJsfAvObof6reG536Bey2LtMjoxnblbwtkS\nfJUtb3ShSZ3SbQEjASCEEA8jNxsOLoTdc8HGVr/O3/5F/c9FlJKZzcJdZ/lu3wVsbOC1bh40qF76\n/YDM8i4gcxcdHU23bt3w9vbGx8eHBQsWPNT2Xbt2vdMfqCSMHj2a3377rcT2J4TIR8xRWNIVts8A\njx7wr8PwyMtFPvjn5Obx46GLdJu3m2/+iGJgy3rsersrb/ZqSsXypf/5XM4AisDOzo7PPvuMNm3a\nkJqaStu2benVqxfe3t6l/t4l0S5aWk4L8ZAyk2HHLDiyFKrUg2ErwWtgsXa5OyKOjzaGcSYujfbu\nNfhugBctnKuV0IQLR44CRVCvXj3q1dNv7apSpQpeXl7Exsby6quv0qFDB3bt2kVSUhLLli2jS5cu\nZGRkMGbMGE6dOoWnpycZGRkG95uZmckrr7zC0aNHsbOz4/PPP6dbt26sWLGC1atXk5aWRm5uLrt3\n7+b1119n27ZtuLi4UL58+Tv7OHbsGBMmTCAtLY2aNWuyYsUK6tWrR9euXf/Wmvqtt94yys9KCIum\nFISug83vws046PAydJ8K9kW/Nh95LZWPNobxR+R1GjpVZPFzbenjU8ckT+ez7ADYPAmuBpXsPuv6\nQr8Ht3X4y4ULFzhx4gQdOnQA9E/Xhw8fZtOmTcycOZPt27ezaNEiKlasSFhYGKdPn6ZNmzYG97Vw\n4UI0TSMoKIjw8HB69+5NZGQkAMePH+f06dPUqFGD1atXExERQWhoKNeuXcPb25uxY8eSnZ3N66+/\nzrp166hVqxarVq1i6tSpLF++HOBvramFEA+QFA2b3obILVC3BQz/GRoY/t0tjPi0LL7YFsnPhy9R\nyd6OaQO8GNmxIfZ2Ra8dFJdlB4CJpaWlMXToUObPn0/VqlUBGDJkCABt27blwoULAOzZs4dx48YB\n0KJFC1q0aGFwf/v27eP1118HwNPTk4YNG94JgF69elGjRo07+xs+fDi2trbUr1+f7t27AxAREUFw\ncDC9evUC9KeS/XWmAve3qhZCGJCbA38uhl2zAQW9P9I/+dsW7XCZmZ3Ld/svsHDXWTKyc3m+oxvj\nejShRqXyD964lFl2ADzEJ/WSlp2dzdChQ+/07P+Lvb09oD83OCcnp8B9rFmzhpkzZwKwdOnSAscW\npq20UgofHx8OHjxY5H0IYdVij0PAG3D1NDTtC/3nQTXXIu1KKcXGoCvM3RxOzI0MenjWZnJ/Lzxq\nF2+NQEmSu4CKQCnFCy+8gJeXFxMmTHjg+Mcee4yffvoJgODgYE6fPg3A4MGDOXnyJCdPnsTPz48u\nXbqwcuVKACIjI7l06RLNmjUzuL9Vq1aRm5vLlStX2LVrFwDNmjXj+vXrdwIgOzubkJCQEvmehSjT\nslL1S8pLe0BaHPzjexj+S5EP/icu3eCpxQd57acTVLa348cXOrBsdDuzOviDpZ8BmMj+/fv54Ycf\n8PX1vdMWevbs2fmOf+WVVxgzZgxeXl54eXnRtq3hhz2/+uqrvPLKK/j6+mJnZ8eKFSvunFHcbfDg\nwezcuRNvb29cXV3p2LEjAOXLl+e3335j3LhxJCcnk5OTw/jx4/HxMb8H5whhNsI3wqZ3IOUytHsB\nerwHDkVrRBmblMEnW8JZd/IyNSvbM3eIL0/7uWBrY/wCb2GYZTvov/j5+al7i5bSurjkyc9UWKXk\nWNg8EcI3QG0fvV2zS7si7SotK4dFu8+ydO95AF7s0oiXuzamsr1pPmMXth20nAEIIaxLXi4c/hZ2\nztL/3HMmdPwX2JZ76F3l5il+PRrNp1sjiU/LYlCr+kzs60mDaqW/irckSAAIIazHlVN6kffyCfDo\nCQM+g+puRdrVvjPxfLgxlPCrqbRtWJ1vn29La9fqJTvfUiYBIIQo+7LSYPccOLQIKjrB0GX6Q9mL\nsPjqbFwaczaFsSM8DufqFfh6RGsG+NYzyUKu4rLIAFBKWeQP2xyZcw1IiBIRGah37UyOhrZjoOcM\nqPDwn9Rv3LzF/O2R/PjnJSqUs2VSP09Gd3LDoZzpFnIVl9ECQNO0RsBUwFEp9VRR9+Pg4EBCQgJO\nTk4SAsWklCIhIQEHBwdTT0WIkpdyBba8q7dyqOUJYwPB9ZGH3k1WTi4/HLzIlzvOkJaVw/D2rrzZ\nqyk1K99/h56lKVQAaJq2HBgIxCmlmt/1el9gAWALLFVK5bsySykVBbygaVqx2lY6OzsTExPD9evX\ni7MbcZuDgwPOzs6mnoYQJScvF44uhx0fQE4WdJ8OncaB3cOtvFVKERhylTmbw7mYkM7jTWsxdYAX\nTUu5R78xFfYMYAXwNfD9Xy9ommYLLAR6ATHAEU3T1qOHwZx7th+rlIor9myBcuXK4e7uXhK7EkKU\nNVeD9SJv7FFo1BUGfA5OjR96N0ExyczaGMrh84k0qV2ZFWPa0bVZ7RKfrqkVKgCUUns0TXO75+X2\nwNnbn+zRNO0XYJBSag762YIQQhjHrXT442M4+DU4VIMh34Lv0w9d5L2SnMG8wAhWH4/FqVJ5Pnyy\nOc+0c8HOtmw2TShODaABEH3X1zFAh/wGa5rmBHwEtNY0bfLtoDA07iXgJQBX16ItwxZCWJEz22Hj\nBEi6CK2fg16zoGKNh9pF+q0cFv8RxZI958jLg5cfb8yr3RpT1eHh1wZYEqMVgZVSCcDLhRi3BFgC\n+krg0p6XEMJCpV6DwMkQ/DvUbAqjN4Lbow+1i7w8xe/HY5gXGEFcahYDWtRjUl9PozyQ3RwUJwBi\nAZe7vna+/ZoQQpSevDw4/h/9sYzZGdB1Cjw6Huwe7q6cg+cS+HBjKCGXU2jpUo1Fz7WhbcOHO3Ow\ndMUJgCNAE03T3NEP/M8AI0pkVkIIYUhcGASMh+hD4NYFBn4BNZs81C7Ox99kzqYwtoZeo76jAwue\naYV/i/rYmGnDttJU2NtAfwa6AjU1TYsBZiillmma9hoQiH7nz3KlVIn0HtY0zR/w9/DwKIndCSEs\nXXYG7JkH+xfoj2Mc9G9oNeKhirzJ6dks2HGG7w9ewN7Ohnf6NOOFR90teiFXcVlcN1AhhJU5t0sv\n8iZGQcsR0HsWVKpZ6M2zc/P48dBFFuw4Q3JGNsP8XJjQuym1q5TdBZDSDVQIYdluxkPgFDi9Cmo0\nhufXQ6PHC725UortYXHM2RRGVPxNOns4MbW/N971q5bipC2LBIAQwrwoBSd+hG3T9SZuj02ELm9B\nucJ/Yg+5nMxHG8M4cC6BRrUqsXy0H92a1Zb2MfcwywCQGoAQVup6JGwYDxf3g2tHGDgfansWevO4\nlEw+3RrBr8dicKxQjplP+DCigyvlyuhCruIyywBQSgUAAX5+fi+aei5CCCPIzoR9X8C+z6FcBfD/\nElqPBJvCHbgzbuXy7d4oFv9xjuzcPP7vUXde69YEx4pleyFXcZllAAghrMj5vfqn/oSzevuGPrOh\ncuH67uTlKdadiuWTLRFcSc6kr09dJvf3pKFTpVKedNkgASCEMI30RNg6DU6u1J/K9dxq8OhR6M2P\nXEjkww2hnIpJxreBI/OHtaJDI6fSm28ZZJYBIDUAIcowpeDUL7B1KmQmw6MT4PGJ+qWfQriUkM7c\nLWFsCrpK3aoOfP6PljzZqoFVLuQqLlkHIIQwnoRz+uWe83vAuT34z4c6PoXaNDkjm4W7zrJi/wVs\nbTRefrwxLz7mTsXyZvk51qRkHYAQwnzk3NJX8e6ZB3YOep/+tmMKVeTNyc3j58OX+GL7GW6k32Jo\nG2fe6dOMOlXL7kIuY5EAEEKUrosH9Ye0xEeAz2DoOxeq1H3gZkopdkdc56NNYZyNS6ODew2mD/Sm\neQNHI0zaOkgACCFKR3qi3rHz+Pfg6AojfoWmvQu1afjVFD7aGMbeM/G4OVVkyci29PKuIwu5SpgE\ngBCiZCkFQb/pvfrTE/Xn8XadBOUffGvm9dQsPt8Wyaojl6hsb8f0gd6MfKQh5e1kIVdpMMsAkLuA\nhLBQiVGw8S04txMatIWRa6Cu7wM3y8zOZfn+8/x71zkys3N5vqMbb/RoQvVKD/cgd/FwzDIAZCWw\nEBYmNxsOfKU/l9emHPSbB+1eAJuCWy0rpQg4fYWPN4cTm5RBT686TO7vSeNalY00cetmlgEghLAg\n0Yf1Im9cKHg9Af0+hqr1H7jZ8Us3mLUhlBOXkvCqV5V5T7Wgk0fh2zyL4pMAEEIUTUYS7JgJR7+D\nqg1g+C/QrN8DN4tOTOeTwAgCTl2mVhV7PhnagqFtnbGVhVxGJwEghHg4SkHIGtgyCW5eh0dehW5T\nwL7gyzapmdn8e/c5lu07j40G47p78M/HG1PJXg5DpiI/eSFE4d24CJvehjNboV5LGLEK6rcucJPc\nPMWqI9F8vi2C+LRbDG7dgHf6NKN+tcK1fhClxywDQO4CEsLM5GbDoUWwew6gQZ850P4lsC34EHIy\nOompa4IIuZxCO7fqLBvVjpYu1YwzZ/FAZhkAcheQEGYk5phe5L0WBM36Q/954Ohc4CbJ6dl8HBjO\nz4cvUauyPV8Nb83AFvVkIZeZMcsAEEKYgcwU2DkLDn+rt24Y9iN4DoQCDuJKKX4/HsucTWHcSL/F\nmE7uvNmrCVUc5MEs5kgCQAjxd0pBWABsngipV/VLPd2ngUPBD1OPvJbKtDXBHL6QSBvXanz/Qnt8\n6kvfHnMmASCE+J+kaNj0DkRuhjq+MGwlOLctcJObWTl8ueMMy/adp7KDHXOH+PIPPxfpz28BJACE\nEJCbA4eXwM4PAQW9Zum3dxZQ5FVKERhyjQ8CQricnMk//JyZ1M+LGtK+wWJIAAhh7S6f0Iu8V05B\nk97Q/1Oo3rDATaIT05mxPoSd4XE0q1OFX4e3pp1bDSNNWJQUCQAhrFVWKuyaDX8uhkq14OkV4P1k\ngUXerJxcvt0TxVc7z2JrozG1vxejO7tRzla6dVoiswwAWQcgRCkL36Rf60+JBb+x0OM9qFDw/fn7\nz8YzfV0wUddv0t+3LtMHelPPURZzWTKzDABZByBEKUm5rN/dExYAtb3h6e/ApX2Bm8SlZPLhxjDW\nn7qMa42KfDemHd2a1TbShEVpMssAEEKUsLxcOLIUdsyCvGzoMQM6vQ62+d+fn5un+OHgBT7bGklW\nTh7jejTh1a6NcShXcItnYTkkAIQo666chg3jIfYYNO6uP5C9hnuBm9zdwqFLk5p8MKg57jUf/EQv\nYVkkAIQoq27d1Hv3HPw3VKwBQ5dB86EFFnmT07P5JDCcn263cPh6RGsG+EoLh7JKAkCIsihyq/5o\nxuRL0GYU9HxfD4F8KKVYfTyW2dLCwapIAAhRlqRehc3vQuhaqNkMxmyGhp0K3CTyWirT1gZz+Hwi\nraWFg1WRABCiLMjLg2PLYftMyMmCbtOg8xtgl/+q3PRbOSzYcYZle89Tyd6OOUN8GSYtHKyKBIAQ\nlu5aiL6SN+YIuD8GA+eDU+N8hyul2Bp6jZnr/9fC4d2+njhVtjfipIU5MMsAkIVgQhTCrXTY8wkc\n+AocHGHwN9BiWIFF3ujEdN5fH8IOaeEgAE0pZeo55MvPz08dPXrU1NMQwvyc3Q4bJkDSRWj1HPT6\nACo55Ts8KyeXpXvP89XOM9hoGm/2bCotHMowTdOOKaX8HjTOLM8AhBD5SIuDLZMh+Ddw8oBRG8C9\nS4GbHDgbz7TbLRz6NddbOMjzeAVIAAhhGfLy4MT3sO09yM6AxydBlwlgl/91+7jUTD7aGMa6k9LC\nQRgmASCEuYsL11fyXjoIDR+FgV9Arab5Ds/NU/x46CKfBkZICwdRIAkAIcxVdibs/RT2zQf7yjBo\nIbR6tsAi76noJKauDSI4NoVHPWrywSAfGtWqbMRJC0siASCEOYraDRvehMQoaPEM9PkIKtXMd3hy\nejbztoaz8k9p4SAKTwJACHNyMx4Cp8LpX6BGIxi5Fhp3y3f4vS0cRndyY0KvptLCQRSKBIAQ5kAp\nOLkStk6DrDTo8jY89jaUy/9uHWnhIIpLAkAIU4s/AwHj4eI+cHkE/OdDba98h6ffyuHLHWdZujdK\nWjiIYpEAEMJUcrJg7+ew73P9k77/Amj9PNjkvzhra8hVZgaEEpuUwdNtnZnUT1o4iKKTABDCFM7v\n1Yu8CWeg+VPQdw5Uzv8e/ejEdGYGhLA97HYLh5c7SgsHUWwSAEIYU3oibJ0OJ3+Eag3h2d+hSc98\nh9/KyePbvVF3WjhM7e8lLRxEiZEAEMIYlILTqyBwCmQmw6NvwmMToXzFfDc5cDae6euCOXf9Jn19\n6vKev7RwECXLLANAuoGKMiXhnH655/wf4NxOv9Zfxyff4XGpmczeGMbav1o4jG5HN09p4SBKnlkG\ngFIqAAjw8/N70dRzEaLIcm7B/gWwZ57es2fAZ9B2bL5F3tw8xco/LzIvMIKs7DzGdffg1W4e0sJB\nlBqzDAAhLN7Fg/pDWuIjwPtJ6DsXqtbLd/ip6CSmrQ0mKDZZWjgIo5EAEKIkZdyAbTPg+H/A0RVG\n/Bea9sl3+L0tHL4a3pqBLaSFgzAOCQAhSoJSEPw7bJmk3+nT8TXoNgXKV8pnuGLNCb2FQ+JNaeEg\nTEMCQIjiSjwPG9+Cczugfht47neo1zLf4Wdut3D483wirVyqsWJMe5o3kBYOwvgkAIQoqtxs/Xm8\nf3wMNnbQ7xNo939gY7hoe28Lh9mDfXmmnbRwEKYjASBEUUQf1ou8caHgOVA/+Ds2yHf4ttBrvL8+\nhNikDJ5q68xkaeEgzIAEgBAPIzMZts+Eo8uhan145ifwHJDv8HtbOPz3nx1p7y4tHIR5kAAQojCU\ngtC1sPlduHkdOrwM3aeCfRWDw+9t4TClvydjOrtLCwdhViQAhHiQGxdh09twZivUbQHDf4EGbfId\nfuBcPNPXSgsHYf4kAITIT24OHPo37J4DaNBnNrT/J9ga/rW5u4WDS40K0sJBmD0JACEMiTkGG96A\nq0HQtB/0nwfVXAwOlRYOwlJJAAhxt8wU2PkhHF4ClevAP74Hrycgn5W5p2OSmLpGb+HQ2cOJDwY1\np7G0cBAWQgJAiL+EBcCmiZB6Bdq/CN2ngYPhBVrJGdl8GhjBj39epGZle74c3hp/aeEgLIwEgBDJ\nMfqBP2Ij1GkOw34AZz+DQ5VSrDt5mQ83hpJ48xajOroxoXdTqkoLB2GBJACE9crLhT+/0S/5qDzo\n9QE88irYGj6Yn4+/yfS1wew7G09LaeEgygAJAGGdLp/UV/JeOQkePfVe/dXdDA7Nysnlmz+i+HrX\nWextbZg1yIcRHRpiKy0chIWTABDWJSsNds2GPxdBxZrw1HfgMzjfIu+hqASmrgni3PWbDGhRjxkD\nvald1cHIkxaidEgACOsRsRk2vg0pMdB2DPR8HypUMzg08eYt5mwK49djMbjUqMCKMe3o2kzu6Rdl\ni9ECQNO0J4EBQFVgmVJqq7HeW1i5lMuweaJ+l08tLxi7FVw7GByqlOL347F8tDGU1MwcXunamHHd\nm1ChvNzTL8qeQgWApmnLgYFAnFKq+V2v9wUWALbAUqXU3Pz2oZRaC6zVNK068CkgASBKV16u3rRt\n+0zIy4Ye70HH18GuvMHh566nMXVNEIeiEmnbsDqzB/vSrK7hXj9ClAWFPQNYAXwNfP/XC5qm2QIL\ngV5ADHBE07T16GEw557txyql4m7/edrt7YQoPVeD9CJv7DFo1A0Gfg41Ghkcmpmdy6Ld51i0+xwO\n5WykT7+wGoUKAKXUHk3T3O55uT1wVikVBaBp2i/AIKXUHPSzhb/R9BUyc4HNSqnjxZm0EPm6dRN2\nz4WDC6FCdRiyFHyfyrfIe+BsPFPXBnM+/iaDWtVn2gBvalWRPv3COhSnBtAAiL7r6xjA8IVV3etA\nT8BR0zQPpdRiQ4M0TXsJeAnA1dW1GNMTVufMNtg4AZIuQZvnoedMqGi4935CWhYfbQxj9YlYGjpV\n5IcX2tOlSS0jT1gI0zJaEVgp9SXwZSHGLQGWAPj5+anSnpcoA1KvwpbJELIaajaF0ZvArbPBoXl5\nil+PRTN7Uzjpt3J4vbsH/5LGbcJKFScAYoG72yM6335NCOPIy4Nj3+lF3pxM6DYVOr8BdoYv4Zy5\nlsrUNcEcvpBIe7cazB7SHI/aUuQV1qs4AXAEaKJpmjv6gf8ZYESJzEqIB7kWChvGQ/Sf4NYFBs6H\nmh4Gh2Zm5/L1zrN8s+cclezt+GRoC55q6yxFXmH1Cnsb6M9AV6CmpmkxwAyl1DJN014DAtHv/Fmu\nlAopiUlpmuYP+Ht4GP6FFlYsOwP++AQOfAn2VeHJRdByeL5F3j2R15m+LpiLCekMadOAqf295GHs\nQtymKWW+l9n9/PzU0aNHTT0NYS7O7YQNb8KNC9DqWeg1Cyo5GRwal5rJhxvCWH/qMo1qVuLDJ5vT\nyaOmcecrhIlomnZMKWW4pe1dpBWEMH9p1yFwMgT9CjUaw6gAcH/M4NC8PMXPRy7x8eZwMrPzGN+z\nCS8/3liKvEIYIAEgzFdeHpz8EbZO1+/vf/xdeHQClDPcjC38agpTVgdx/FISHRs58eFgeTqXEAUx\nywCQGoDgegQEjIdLB6BhZxj4BdRqZnBoxq1cFuw4w9K9UVRxsOOzp1sypE0DeTqXEA9glgGglAoA\nAvz8/F409VyEkWVnwt7PYN8XUL4SPPG1fr3fxsbg8F0RcUxfG0zMjQz+4efM5H5eVK9kuNePEOLv\nzDIAhJWK+kMv8iaeA99/QJ/ZUNnw6txrKZl8EBDKxqArNK5ViVUvPUKHRoYLwkIIwyQAhOndTICt\n0+DUT1DdHUaugcbdDQ7NzVOs/PMi87ZEkJWbx1u9mvLS442wt5MirxAPyywDQGoAVkIpOPmTfvDP\nSoEub8Fj70C5CgaHh1xOZsqaYE5FJ/GoR00+fLI5bjUrGXnSQpQdZhkAUgOwAvFn9ZW8F/aCSwfw\nXwC1vQwOvZmVw/ztkSzff4HqFcux4JlWPNGyvhR5hSgmswwAUYblZMG++bD3U7CroLdwaDMq3yLv\n9tBrzFgfQmxSBsPbu/BuX0+qVZQirxAlQQJAGM+F/fqn/vhIaD4U+syBKnUMDr2SnMHM9aFsCblK\n0zqV+e3ljvi5GW7tLIQoGgkAUfrSE2Hbe3DiB6jmCs/+Bk16GRyam6f4/uAFPg2MIFcpJvZtxv89\n2ojydobPEIQQRScBIEqPUnr7hi2TIeOG3qr58UlQvqLB4UExyUxZE0RQbDKPN63FrEHNcXUyPFYI\nUXxmGQByF1AZkHBOfzpX1G5o4AfPr4O6zQ0OTcvK4bOtEfznwAWcKtvz9YjWDPCtJ0VeIUqZWQaA\n3AVkwXJe/USyAAAPrUlEQVRu6a2a98wD2/LQ/1PwGws2hu/TDwy5yvvrQ7iaksmzHVx5p48njhXK\nGXnSQlgnswwAYaEuHdL791wPA+9B0PdjqFrP4NDYpAxmrAthe9g1POtWYeGzbWjjWt3IExbCukkA\niOLLuKE/lvHYd+DoAsNXQbO+Bofm5Oax4sAFPt8WiVIwpb8nYzq7U85WirxCGJsEgCg6pfQHsW+e\nBOnx0PE16DoZ7A23YD4ZncSU1UGEXkmhu2dtPhjkg3N1KfIKYSoSAKJoblyAjW/B2e1QvzU8+yvU\nb2VwaEpmNp8FRvD9oYvUrmLPomfb0Ld5XSnyCmFiZhkAcheQGcvNhoMLYfdcvbDb92No/6LBIq9S\nis3BepH3eloWozq68VbvplRxkCKvEObALANA7gIyUzFHIeANuBYMzQZA/0/A0dng0OjEdN5bF8yu\niOv41K/Kt8/70dKlmpEnLIQoiFkGgDAzmcmwYxYcWQpV6sGwleA10ODQ7Nw8lu07z/ztkdhoGtMG\neDG6kxt2UuQVwuxIAIj8KQVh62Hzu5B6FTr8E7pPA/sqBocfu3iDqWuCCL+aSm/vOrz/hA/1qxlu\n7SyEMD0JAGFYUjRsehsit0BdX3hmJTRoa3Bocno2HweG8/PhS9St6sA3I9vSx6eukScshHhYEgDi\n73Jz4M/FsGs2oKD3R9DhZbC9/5+KUoqA01f4ICCUxJtZjO3szpu9mlLZXv5ZCWEJ5DdV/M/lE3qR\n98opaNIHBnyqd+804GLCTaatDWbvmXhaODuyYkw7mjdwNPKEhRDFIQEgICsVdn4Eh7+BSrXh6f/o\nrRwM3Kd/KyePb/dG8eWOM5SzteF9f29GdnTD1kbu6RfC0phlAMg6ACMK3wib3oGUy9DuBejxHjgY\n/iR/5EIiU1YHcSYujX7N6zLD34e6jg5GnrAQoqSYZQDIOgAjSI6FzRMhfAPU9tE/9bu0Mzg0Kf0W\nczeH88uRaBpUq8CyUX708DL8JC8hhOUwywAQpSgvFw5/Cztn6X/uORM6/gts71+dq5Ri7clYPtwQ\nRlJGNi891ojxPZtQsbz8sxGiLJDfZGty5ZRe5L18Ahr3gAGfQQ13g0PPx99k2tog9p9NoJVLNX4Y\n7It3/apGnrAQojRJAFiDrDTYPQcOLYKKTjB0mf5QdgNF3qycXL75I4qvd53F3s6GWU82Z0R7Vyny\nClEGSQCUdZGBetfO5GhoOxp6vg8VDD945VBUAlPWBBF1/SYDW9TjvYHe1K4qRV4hyioJgLIq5Qps\neRdC10EtTxgbCK6PGByaePMWszeF8duxGFxqVGDFmHZ0bVbbyBMWQhibBEBZk5cLR5fDjg8gJwu6\nT4dO48Cu/H1DlVL8fjyWjzaGkpqZw6tdG/N69yZUKG/4+b1CiLJFAqAsuRqsF3ljj4L74zDwC3Bq\nbHDouetpTF0TxKGoRNo2rM7swb40q2u4yZsQomySACgLbqXDHx/Dwa/1RVyDl0CLf+Rb5F28O4qF\nu85iX86G2YN9eaadCzZS5BXC6phlAMhK4IdwZjtsnABJF6H1c9BrFlSsYXDo3UVe/5b1mT7Qi9pV\npMgrhLUyywCQlcCFkHoNAidD8O/g1ARGbwS3Rw0OvXG7yPurFHmFEHcxywAQBcjLg+P/ge0zIDsD\nuk6GR98EO/v7hiqlWHMilg83hpGSkc0rXRszToq8QojbJAAsSVwYBIyH6EPg1kUv8tZsYnDo+fib\nTF0TxIFzCbR2rcacIb541pWVvEKI/5EAsATZGbBnHuxfoD+OcdC/odWIQq3k/fD2Sl4p8goh7iUB\nYO7O7dKLvIlR0HI49P4QKtU0OPTw+USmrAnibFwaA1rUY4as5BVCFEACwFzdjIfAKXB6FdRoDM+v\nh0aPGxyalH6LOZvCWXU0GufqFfhuTDu6SZFXCPEAEgDmRik48SNsm643cXtsInR5C8rd/0leKcW6\nk5eZtSGUpIxs/vl4I8b3aCpFXiFEoUgAmJPrkbBhPFzcD64dYeB8qO1pcOiFeP2ZvPvOxtPKpRo/\nDvHFq54UeYUQhScBYA6yM2HfF7DvcyhXAfy/hNYjwcbmvqG3cvJYsuccX+48i72ttGsWQhSdBICp\nnd8DG96EhLPg+zT0mQ2VDV+/P3Ihkcmrbxd5fevxnr83daTIK4QoIgkAU0lPhK3T4ORKqO4Gz60G\njx4Gh977TN7lo/3o7inP5BVCFI8EgLEpBad+ga1TITMZHp0Aj70D5SsaGKpYf0ov8t5Il2fyCiFK\nlhxJjCnhnF7kPb8HnNuD/3yo42Nw6MUEvci790w8LV2q8Z+xzfGp72jkCQshyjKzDIAy1w0055a+\ninfPPLBzgAGfQ9sx+RZ5v90bxZc7zlDO1oYPBvnwbIeGUuQVQpQ4swyAMtUN9OIBvX9PfAT4DIa+\nc6FKXYNDj17QV/JGXkujv29d3hvoQ11HKfIKIUqHWQZAmZCeqHfsPP49OLrCiF+haW+DQ5PTs5m7\nJZyfD1+iQbUKLBvlRw8vKfIKIUqXBEBJUwqCftN79acn6s/j7ToJylcyMFQRcPoKHwSEciP9Fi92\ncWd8z6ZUspf/LUKI0idHmpKUGAUb34JzO6F+G/3WznotDA69lJDOtHXB7Im8TktnR1aMaUfzBlLk\nFUIYjwRAScjNhgNfwh+fgE056DcP2r0ANvf35MnO1Yu8C7brRd73/b0Z2dFNirxCCKOTACiu6MMQ\n8AbEhYKXP/T7BKrWNzj02MVEpqwOJuJaKn196jLjCW/qOVYw8oSFEEInAVBUGUmwYyYc/U4/4D/z\nM3j2Nzg0OSObT7aEs/LPS9R3dODb5/3o5S1FXiGEaUkAPCylIGQNbJkEN6/DI69Ct8n6k7ruG6rY\ncPoKMwNCSbyZxQuPujOhlxR5hRDmQY5ED+PGRdj0NpzZCvVawohVUL+1waHRielMWxvMH5HX8W0g\nRV4hhPmRACiM3Gw49G/YNQc0G+gzB9q/BLb3//iyc/NYtu8887dHYqtpzPD35nkp8gohzJAEwIPE\nHNOLvNeCoGk/6D8PqrkYHHr80g2mrA4i/GoqfXzq8P4TPlLkFUKYLQmA/GSmwM5ZcPhbvXXDsB/B\ncyBo93+ST8n8X5G3blUHloxsS28fw+0ehBDCXEgA3EspCAuAzRMh9Sq0fxG6TweH+x+3qJRiU9BV\n3g8IISEtizGd3JnQuymVpcgrhLAAcqS6W1I0bHoHIjdDHV8YthKc2xocGp2YznvrgtkVcZ3mDaqy\nfFQ7fJ2lyCuEsBwSAAC5OXD4G9j5EaCg1yz99s58irzL953ni+2R2Gga0wd6M6pjQ+xs72/tLIQQ\n5kwC4PIJvch75RR49IIBn0H1hgaHnrh0g8m3i7w9verwwSAf6leTIq8QwjJZbwBkpcKu2fDnYqhU\nC55eAd5P5lvk/TQwgh8OXaROFQe+GdmWPlLkFUJYOOsMgPBN+rX+lFjwGws93oMK1e4bppRic/BV\n3l8fwvW0LEZ1dOOt3k2p4lDOBJMWQoiSZbQA0DTNC3gDqAnsUEotMtZ735FyWT/wh2+A2t7w9Hfg\n0t7g0Jgb6by3LoSd4XH41K/Kt8/70dLl/pAQQghLVagA0DRtOTAQiFNKNb/r9b7AAsAWWKqUmpvf\nPpRSYcDLmqbZAN8DxguAvFw4shR2zIK8bOgxAzq9Drb3f5LPyc3ju/0X+HxbJJoG0wZ4MbqTmxR5\nhRBlTmHPAFYAX6MfuAHQNM0WWAj0AmKAI5qmrUcPgzn3bD9WKRWnadoTwCvAD8Wcd+FdOa0XeS8f\nh8bd9SJvjUYGh56MTmLK6iBCr6TQ06s2Mwc1p4EUeYUQZVShAkAptUfTNLd7Xm4PnFVKRQFomvYL\nMEgpNQf9bMHQftYD6zVN2wj8VNRJF8qtm7B7Dhz8N1SsAUOWgu9TBou8qZnZfLY1kv8cvEDtKvYs\nfq4NfXzqohkYK4QQZUVxagANgOi7vo4BOuQ3WNO0rsAQwB7YVMC4l4CXAFxdXYs2s8it+qMZky9B\nm+eh50w9BO6hlCIw5Coz1ocQlypFXiGEdTFaEVgptRvYXYhxS4AlAH5+fqpIb3bkWyhXAcZshoad\nDA6JTcpgxrpgtofF4V2vKt+M9KOVFHmFEFakOAEQC9zdFtP59mum9+RisK8Mdvb3/VVObh4rDuhF\nXqVgan8vxnSWIq8QwvoUJwCOAE00TXNHP/A/A4wokVkVVyUngy+fjkli8uogQi6n0N2zNh8M8sG5\nekUjT04IIcxDYW8D/RnoCtTUNC0GmKGUWqZp2mtAIPqdP8uVUiElMSlN0/wBfw8Pj5LY3Z0i7/cH\nL1Czsj2Lnm1D3+ZS5BVCWDdNqaJdZjcGPz8/dfTo0WLtY8vtlbzXUjMZ+UhD3u7TjKpS5BVClGGa\nph1TSvk9aFyZbQVxOSmDGetD2BZ6Dc+6VVj0XBtau1Y39bSEEMJslMkAWL7vPJ9ujUApmNLfkzGd\n3SknRV4hhPgbswyA4tYAwq+m0MG9Bh8Mao5LDSnyCiGEIWWyBpCVk0t5Wxsp8gohrJJV1wDs7WxN\nPQUhhDB7cmFcCCGslASAEEJYKbMMAE3T/DVNW5KcnGzqqQghRJlllgGglApQSr3k6Oho6qkIIUSZ\nZZYBIIQQovRJAAghhJWSABBCCCtllusA/loJDKRomnamkJs5AqaqGhvjvUv6PUpif0XdR1G2e9ht\nagLxD/ke1syUvz9FYer5lvb7F3f/DQs1SilVJv4DlpTl9y7p9yiJ/RV1H0XZ7mG3AY6a6t+DJf5n\nyt8fS5xvab+/sb6/snQJKKCMv3dJv0dJ7K+o+yjKdqb8/2sNLO3na+r5lvb7G+X7M+teQEIUlaZp\nR1UheqEIYc3K0hmAEHdbYuoJCGHu5AxACCGslJwBCCGElZIAEEIIKyUBIIQQVkoCQFgVTdMaaZq2\nTNO030w9FyFMTQJAWAxN05ZrmhanaVrwPa/31TQtQtO0s5qmTSpoH0qpKKXUC6U7UyEsg1m2ghAi\nHyuAr4Hv/3pB0zRbYCHQC4gBjmiath6wBebcs/1YpVSccaYqhPmTABAWQym1R9M0t3tebg+cVUpF\nAWia9gswSCk1Bxho3BkKYVnkEpCwdA2A6Lu+jrn9mkGapjlpmrYYaK1p2uTSnpwQ5kzOAIRVUUol\nAC+beh5CmAM5AxCWLhZwuetr59uvCSEeQAJAWLojQBNN09w1TSsPPAOsN/GchLAIEgDCYmia9jNw\nEGimaVqMpmkvKKVygNeAQCAM+K9SKsSU8xTCUkgzOCGEsFJyBiCEEFZKAkAIIayUBIAQQlgpCQAh\nhLBSEgBCCGGlJACEEMJKSQAIIYSVkgAQQggrJQEghBBW6v8BrUmKQrRTWC8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x110ff9a90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.loglog(h_values,err_values)\n",
    "plt.loglog(h_values,np.array(h_values)**2)\n",
    "plt.legend(['Errors','2nd-order'])"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
 }
	{
	"cells": [
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": [
	"# Simple implementation of BDF2 for a linear, scalar ODE"
	]
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": [
	"Note that the BDF2 formula is\n",
	"\n",
	"$$y^{n+1} = \\frac{4}{3} y^n - \\frac{1}{3} y^{n-1} + \\frac{2}{3} h f(y^{n+1})$$\n",
	"\n",
	"For a scalar ODE with $f(y) = \\alpha y$, this means\n",
	"\n",
	"$$y^{n+1} = \\frac{4}{3} y^n - \\frac{1}{3} y^{n-1} + \\frac{2}{3} h \\alpha y^{n+1},$$"
	]
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": [
	"so\n",
	"\n",
	"$$y^{n+1} = \\left(1-\\frac{2}{3}h\\alpha\\right)^{-1}\\left(\\frac{4}{3} y^n - \\frac{1}{3} y^{n-1}\\right).$$"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 8,
	"metadata": {
	"collapsed": true,
	"deletable": true,
	"editable": true
	},
	"outputs": [],
	"source": [
	"def bdf2(h, T, alpha, y0):\n",
	" # Variable meanings:\n",
	" # y_nm1, y_n, y_np1: y_{n-1}, y_n, y_{n+1}\n",
	" # h: Time step\n",
	" # T: Final time\n",
	" # y0: Initial value\n",
	" \n",
	" t = 0\n",
	" times = [0]\n",
	" y_n = y0\n",
	" \n",
	" # Euler step to start\n",
	" y_np1 = y_n + halphay_n\n",
	" \n",
	" t += h\n",
	" times.append(t)\n",
	" y_nm1, y_n = y_n, y_np1\n",
	" \n",
	" while t<T:\n",
	" if T-t < h: # Don't step past the end time\n",
	" h = T-t\n",
	" \n",
	" y_np1 = ( (4./3)y_n - (1./3)y_nm1 ) / (1.-(2./3)halpha)\n",
	" t += h\n",
	" times.append(t)\n",
	" y_nm1, y_n = y_n, y_np1\n",
	"\n",
	" exact_solution = np.exp(Talpha)y0\n",
	" error = np.abs(y_n - exact_solution)\n",
	" return y_n, error"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 30,
	"metadata": {
	"collapsed": false
	},
	"outputs": [],
	"source": [
	"y0 = 1.\n",
	"T = 1.\n",
	"\n",
	"alpha = 1.j\n",
	"y_values = []\n",
	"err_values = []\n",
	"h_values = (0.5, 0.25, 0.125, 1./16, 1./32,)\n",
	"\n",
	"for h in h_values:\n",
	" y, error = bdf2(h, T, alpha, y0)\n",
	" y_values.append(y)\n",
	" err_values.append(error)"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 34,
	"metadata": {
	"collapsed": false
	},
	"outputs": [
	{
	"data": {
	"text/plain": [
	"<matplotlib.legend.Legend at 0x110ee2cd0>"
	]
	},
	"execution_count": 34,
	"metadata": {},
	"output_type": "execute_result"
	},
	{
	"data": {
	"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEACAYAAAC6d6FnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlYVeXax/HvAhQcUXEWEBSVQRxRU7OcZzK1jmmZQ2+d\n6pSZlTlmZqllg1YezdQ8lZWncsIJ55xzVmZUHAAHBGSSQYbn/WOZx3SDyLAH9v25rq5Lts9a+4Fk\n/fZa93rupSmlEEIIYX1sTD0BIYQQpiEBIIQQVkoCQAghrJQEgBBCWCkJACGEsFISAEIIYaUkAIQQ\nwkpJAAghhJWSABBCCCslASCEEFbKztQTKEjNmjWVm5ubqachhBAW5dixY/FKqVoPGmfWAeDm5sbR\no0dNPQ0hhLAomqZdLMw4uQQkhBBWyiwDQNM0f03TliQnJ5t6KkIIUWaZZQAopQKUUi85OjqaeipC\nCFFmmXUNwJDs7GxiYmLIzMw09VTKBAcHB5ydnSlXrpyppyKEMDKLC4CYmBiqVKmCm5sbmqaZejoW\nTSlFQkICMTExuLu7m3o6QggjM8tLQAXJzMzEyclJDv4lQNM0nJyc5GxKCCtlcQEAyMG/BMnPUggz\nk50Jp1aBER7Xa3GXgMyBra0tvr6+d75+5plnmDRpkglnJIQoE6L+gA1vQuI5cPIA57al+nYSAEVQ\noUIFTp48WeCY3NxcbG1t73ydk5ODnd2Df9yFHSeEKENuJsDWqXDqZ6juDiPXlvrBHyQASpSbmxvD\nhg1j27ZtTJw4kcWLF9OqVSv27dvH8OHDGTp0KGPHjiU+Pp5atWrx3Xff4erqyujRo3FwcODEiRN0\n7tyZQYMG8cYbbwD6JZo9e/ZQpUoVE393QogSpxSc/Am2ToOsFOjyNjz2NpSrYJS3t+gAmBkQQujl\nlBLdp3f9qszw9ylwTEZGBq1atbrz9eTJkxk2bBgATk5OHD9+HIDFixdz69atO+0s/P39GTVqFKNG\njWL58uWMGzeOtWvXAvrdTQcOHMDW1hZ/f38WLlxI586dSUtLw8HBoUS/RyGEGYg/o1/uubAXXDqA\n/wKo7WXUKVh0AJhKQZeA/goCQ18fPHiQ1atXAzBy5EgmTpx45++efvrpO5eMOnfuzIQJE3j22WcZ\nMmQIzs7OJf0tCCFMJScL9s2HvZ/qn/QHzoc2o8DG+PfkWHQAPOiTuilUqlSpwK8Ls92kSZMYMGAA\nmzZtonPnzgQGBuLp6Vmi8xRCmMCF/bBhPMRHQvOh0GcOVKljsulY5G2glqpTp0788ssvAKxcuZIu\nXboYHHfu3Dl8fX159913adeuHeHh4cacphCipKUnwrrXYEV//Qzg2d/hqeUmPfiDhZ8BmMq9NYC+\nffsyd+7cB2731VdfMWbMGObNm3enCGzI/Pnz2bVrFzY2Nvj4+NCvX78Sm7sQwoiUgtP/hcApkHED\nOr8Bj0+C8hVNPTMANGWExQYPS9M0f8Dfw8PjxTNnzvzt78LCwvDyMm6hpKyTn6kQpSDhHGycAFG7\noYGfXuSt29wob61p2jGllN+DxpnlJSDpBiqEsFg5t2DPp7CoE8Qeh/6fwgtbjXbwfxhyCUgIIUrK\npUMQMB6uh4H3IOj7MVStZ+pZ5UsCQAghiivjBmx/H46tAEcXGL4KmvU19aweSAJACCGKSikIWQ2b\nJ0F6PHR8DbpOBvvKxdptdGI6LjVKv1AsASCEEEVx4wJsfAvObof6reG536Bey2LtMjoxnblbwtkS\nfJUtb3ShSZ3SbQEjASCEEA8jNxsOLoTdc8HGVr/O3/5F/c9FlJKZzcJdZ/lu3wVsbOC1bh40qF76\n/YDM8i4gcxcdHU23bt3w9vbGx8eHBQsWPNT2Xbt2vdMfqCSMHj2a3377rcT2J4TIR8xRWNIVts8A\njx7wr8PwyMtFPvjn5Obx46GLdJu3m2/+iGJgy3rsersrb/ZqSsXypf/5XM4AisDOzo7PPvuMNm3a\nkJqaStu2benVqxfe3t6l/t4l0S5aWk4L8ZAyk2HHLDiyFKrUg2ErwWtgsXa5OyKOjzaGcSYujfbu\nNfhugBctnKuV0IQLR44CRVCvXj3q1dNv7apSpQpeXl7Exsby6quv0qFDB3bt2kVSUhLLli2jS5cu\nZGRkMGbMGE6dOoWnpycZGRkG95uZmckrr7zC0aNHsbOz4/PPP6dbt26sWLGC1atXk5aWRm5uLrt3\n7+b1119n27ZtuLi4UL58+Tv7OHbsGBMmTCAtLY2aNWuyYsUK6tWrR9euXf/Wmvqtt94yys9KCIum\nFISug83vws046PAydJ8K9kW/Nh95LZWPNobxR+R1GjpVZPFzbenjU8ckT+ez7ADYPAmuBpXsPuv6\nQr8Ht3X4y4ULFzhx4gQdOnQA9E/Xhw8fZtOmTcycOZPt27ezaNEiKlasSFhYGKdPn6ZNmzYG97Vw\n4UI0TSMoKIjw8HB69+5NZGQkAMePH+f06dPUqFGD1atXExERQWhoKNeuXcPb25uxY8eSnZ3N66+/\nzrp166hVqxarVq1i6tSpLF++HOBvramFEA+QFA2b3obILVC3BQz/GRoY/t0tjPi0LL7YFsnPhy9R\nyd6OaQO8GNmxIfZ2Ra8dFJdlB4CJpaWlMXToUObPn0/VqlUBGDJkCABt27blwoULAOzZs4dx48YB\n0KJFC1q0aGFwf/v27eP1118HwNPTk4YNG94JgF69elGjRo07+xs+fDi2trbUr1+f7t27AxAREUFw\ncDC9evUC9KeS/XWmAve3qhZCGJCbA38uhl2zAQW9P9I/+dsW7XCZmZ3Ld/svsHDXWTKyc3m+oxvj\nejShRqXyD964lFl2ADzEJ/WSlp2dzdChQ+/07P+Lvb09oD83OCcnp8B9rFmzhpkzZwKwdOnSAscW\npq20UgofHx8OHjxY5H0IYdVij0PAG3D1NDTtC/3nQTXXIu1KKcXGoCvM3RxOzI0MenjWZnJ/Lzxq\nF2+NQEmSu4CKQCnFCy+8gJeXFxMmTHjg+Mcee4yffvoJgODgYE6fPg3A4MGDOXnyJCdPnsTPz48u\nXbqwcuVKACIjI7l06RLNmjUzuL9Vq1aRm5vLlStX2LVrFwDNmjXj+vXrdwIgOzubkJCQEvmehSjT\nslL1S8pLe0BaHPzjexj+S5EP/icu3eCpxQd57acTVLa348cXOrBsdDuzOviDpZ8BmMj+/fv54Ycf\n8PX1vdMWevbs2fmOf+WVVxgzZgxeXl54eXnRtq3hhz2/+uqrvPLKK/j6+mJnZ8eKFSvunFHcbfDg\nwezcuRNvb29cXV3p2LEjAOXLl+e3335j3LhxJCcnk5OTw/jx4/HxMb8H5whhNsI3wqZ3IOUytHsB\nerwHDkVrRBmblMEnW8JZd/IyNSvbM3eIL0/7uWBrY/wCb2GYZTvov/j5+al7i5bSurjkyc9UWKXk\nWNg8EcI3QG0fvV2zS7si7SotK4dFu8+ydO95AF7s0oiXuzamsr1pPmMXth20nAEIIaxLXi4c/hZ2\nztL/3HMmdPwX2JZ76F3l5il+PRrNp1sjiU/LYlCr+kzs60mDaqW/irckSAAIIazHlVN6kffyCfDo\nCQM+g+puRdrVvjPxfLgxlPCrqbRtWJ1vn29La9fqJTvfUiYBIIQo+7LSYPccOLQIKjrB0GX6Q9mL\nsPjqbFwaczaFsSM8DufqFfh6RGsG+NYzyUKu4rLIAFBKWeQP2xyZcw1IiBIRGah37UyOhrZjoOcM\nqPDwn9Rv3LzF/O2R/PjnJSqUs2VSP09Gd3LDoZzpFnIVl9ECQNO0RsBUwFEp9VRR9+Pg4EBCQgJO\nTk4SAsWklCIhIQEHBwdTT0WIkpdyBba8q7dyqOUJYwPB9ZGH3k1WTi4/HLzIlzvOkJaVw/D2rrzZ\nqyk1K99/h56lKVQAaJq2HBgIxCmlmt/1el9gAWALLFVK5bsySykVBbygaVqx2lY6OzsTExPD9evX\ni7MbcZuDgwPOzs6mnoYQJScvF44uhx0fQE4WdJ8OncaB3cOtvFVKERhylTmbw7mYkM7jTWsxdYAX\nTUu5R78xFfYMYAXwNfD9Xy9ommYLLAR6ATHAEU3T1qOHwZx7th+rlIor9myBcuXK4e7uXhK7EkKU\nNVeD9SJv7FFo1BUGfA5OjR96N0ExyczaGMrh84k0qV2ZFWPa0bVZ7RKfrqkVKgCUUns0TXO75+X2\nwNnbn+zRNO0XYJBSag762YIQQhjHrXT442M4+DU4VIMh34Lv0w9d5L2SnMG8wAhWH4/FqVJ5Pnyy\nOc+0c8HOtmw2TShODaABEH3X1zFAh/wGa5rmBHwEtNY0bfLtoDA07iXgJQBX16ItwxZCWJEz22Hj\nBEi6CK2fg16zoGKNh9pF+q0cFv8RxZI958jLg5cfb8yr3RpT1eHh1wZYEqMVgZVSCcDLhRi3BFgC\n+krg0p6XEMJCpV6DwMkQ/DvUbAqjN4Lbow+1i7w8xe/HY5gXGEFcahYDWtRjUl9PozyQ3RwUJwBi\nAZe7vna+/ZoQQpSevDw4/h/9sYzZGdB1Cjw6Huwe7q6cg+cS+HBjKCGXU2jpUo1Fz7WhbcOHO3Ow\ndMUJgCNAE03T3NEP/M8AI0pkVkIIYUhcGASMh+hD4NYFBn4BNZs81C7Ox99kzqYwtoZeo76jAwue\naYV/i/rYmGnDttJU2NtAfwa6AjU1TYsBZiillmma9hoQiH7nz3KlVIn0HtY0zR/w9/DwKIndCSEs\nXXYG7JkH+xfoj2Mc9G9oNeKhirzJ6dks2HGG7w9ewN7Ohnf6NOOFR90teiFXcVlcN1AhhJU5t0sv\n8iZGQcsR0HsWVKpZ6M2zc/P48dBFFuw4Q3JGNsP8XJjQuym1q5TdBZDSDVQIYdluxkPgFDi9Cmo0\nhufXQ6PHC725UortYXHM2RRGVPxNOns4MbW/N971q5bipC2LBIAQwrwoBSd+hG3T9SZuj02ELm9B\nucJ/Yg+5nMxHG8M4cC6BRrUqsXy0H92a1Zb2MfcwywCQGoAQVup6JGwYDxf3g2tHGDgfansWevO4\nlEw+3RrBr8dicKxQjplP+DCigyvlyuhCruIyywBQSgUAAX5+fi+aei5CCCPIzoR9X8C+z6FcBfD/\nElqPBJvCHbgzbuXy7d4oFv9xjuzcPP7vUXde69YEx4pleyFXcZllAAghrMj5vfqn/oSzevuGPrOh\ncuH67uTlKdadiuWTLRFcSc6kr09dJvf3pKFTpVKedNkgASCEMI30RNg6DU6u1J/K9dxq8OhR6M2P\nXEjkww2hnIpJxreBI/OHtaJDI6fSm28ZZJYBIDUAIcowpeDUL7B1KmQmw6MT4PGJ+qWfQriUkM7c\nLWFsCrpK3aoOfP6PljzZqoFVLuQqLlkHIIQwnoRz+uWe83vAuT34z4c6PoXaNDkjm4W7zrJi/wVs\nbTRefrwxLz7mTsXyZvk51qRkHYAQwnzk3NJX8e6ZB3YOep/+tmMKVeTNyc3j58OX+GL7GW6k32Jo\nG2fe6dOMOlXL7kIuY5EAEEKUrosH9Ye0xEeAz2DoOxeq1H3gZkopdkdc56NNYZyNS6ODew2mD/Sm\neQNHI0zaOkgACCFKR3qi3rHz+Pfg6AojfoWmvQu1afjVFD7aGMbeM/G4OVVkyci29PKuIwu5SpgE\ngBCiZCkFQb/pvfrTE/Xn8XadBOUffGvm9dQsPt8Wyaojl6hsb8f0gd6MfKQh5e1kIVdpMMsAkLuA\nhLBQiVGw8S04txMatIWRa6Cu7wM3y8zOZfn+8/x71zkys3N5vqMbb/RoQvVKD/cgd/FwzDIAZCWw\nEBYmNxsOfKU/l9emHPSbB+1eAJuCWy0rpQg4fYWPN4cTm5RBT686TO7vSeNalY00cetmlgEghLAg\n0Yf1Im9cKHg9Af0+hqr1H7jZ8Us3mLUhlBOXkvCqV5V5T7Wgk0fh2zyL4pMAEEIUTUYS7JgJR7+D\nqg1g+C/QrN8DN4tOTOeTwAgCTl2mVhV7PhnagqFtnbGVhVxGJwEghHg4SkHIGtgyCW5eh0dehW5T\nwL7gyzapmdn8e/c5lu07j40G47p78M/HG1PJXg5DpiI/eSFE4d24CJvehjNboV5LGLEK6rcucJPc\nPMWqI9F8vi2C+LRbDG7dgHf6NKN+tcK1fhClxywDQO4CEsLM5GbDoUWwew6gQZ850P4lsC34EHIy\nOompa4IIuZxCO7fqLBvVjpYu1YwzZ/FAZhkAcheQEGYk5phe5L0WBM36Q/954Ohc4CbJ6dl8HBjO\nz4cvUauyPV8Nb83AFvVkIZeZMcsAEEKYgcwU2DkLDn+rt24Y9iN4DoQCDuJKKX4/HsucTWHcSL/F\nmE7uvNmrCVUc5MEs5kgCQAjxd0pBWABsngipV/VLPd2ngUPBD1OPvJbKtDXBHL6QSBvXanz/Qnt8\n6kvfHnMmASCE+J+kaNj0DkRuhjq+MGwlOLctcJObWTl8ueMMy/adp7KDHXOH+PIPPxfpz28BJACE\nEJCbA4eXwM4PAQW9Zum3dxZQ5FVKERhyjQ8CQricnMk//JyZ1M+LGtK+wWJIAAhh7S6f0Iu8V05B\nk97Q/1Oo3rDATaIT05mxPoSd4XE0q1OFX4e3pp1bDSNNWJQUCQAhrFVWKuyaDX8uhkq14OkV4P1k\ngUXerJxcvt0TxVc7z2JrozG1vxejO7tRzla6dVoiswwAWQcgRCkL36Rf60+JBb+x0OM9qFDw/fn7\nz8YzfV0wUddv0t+3LtMHelPPURZzWTKzDABZByBEKUm5rN/dExYAtb3h6e/ApX2Bm8SlZPLhxjDW\nn7qMa42KfDemHd2a1TbShEVpMssAEEKUsLxcOLIUdsyCvGzoMQM6vQ62+d+fn5un+OHgBT7bGklW\nTh7jejTh1a6NcShXcItnYTkkAIQo666chg3jIfYYNO6uP5C9hnuBm9zdwqFLk5p8MKg57jUf/EQv\nYVkkAIQoq27d1Hv3HPw3VKwBQ5dB86EFFnmT07P5JDCcn263cPh6RGsG+EoLh7JKAkCIsihyq/5o\nxuRL0GYU9HxfD4F8KKVYfTyW2dLCwapIAAhRlqRehc3vQuhaqNkMxmyGhp0K3CTyWirT1gZz+Hwi\nraWFg1WRABCiLMjLg2PLYftMyMmCbtOg8xtgl/+q3PRbOSzYcYZle89Tyd6OOUN8GSYtHKyKBIAQ\nlu5aiL6SN+YIuD8GA+eDU+N8hyul2Bp6jZnr/9fC4d2+njhVtjfipIU5MMsAkIVgQhTCrXTY8wkc\n+AocHGHwN9BiWIFF3ujEdN5fH8IOaeEgAE0pZeo55MvPz08dPXrU1NMQwvyc3Q4bJkDSRWj1HPT6\nACo55Ts8KyeXpXvP89XOM9hoGm/2bCotHMowTdOOKaX8HjTOLM8AhBD5SIuDLZMh+Ddw8oBRG8C9\nS4GbHDgbz7TbLRz6NddbOMjzeAVIAAhhGfLy4MT3sO09yM6AxydBlwlgl/91+7jUTD7aGMa6k9LC\nQRgmASCEuYsL11fyXjoIDR+FgV9Arab5Ds/NU/x46CKfBkZICwdRIAkAIcxVdibs/RT2zQf7yjBo\nIbR6tsAi76noJKauDSI4NoVHPWrywSAfGtWqbMRJC0siASCEOYraDRvehMQoaPEM9PkIKtXMd3hy\nejbztoaz8k9p4SAKTwJACHNyMx4Cp8LpX6BGIxi5Fhp3y3f4vS0cRndyY0KvptLCQRSKBIAQ5kAp\nOLkStk6DrDTo8jY89jaUy/9uHWnhIIpLAkAIU4s/AwHj4eI+cHkE/OdDba98h6ffyuHLHWdZujdK\nWjiIYpEAEMJUcrJg7+ew73P9k77/Amj9PNjkvzhra8hVZgaEEpuUwdNtnZnUT1o4iKKTABDCFM7v\n1Yu8CWeg+VPQdw5Uzv8e/ejEdGYGhLA97HYLh5c7SgsHUWwSAEIYU3oibJ0OJ3+Eag3h2d+hSc98\nh9/KyePbvVF3WjhM7e8lLRxEiZEAEMIYlILTqyBwCmQmw6NvwmMToXzFfDc5cDae6euCOXf9Jn19\n6vKev7RwECXLLANAuoGKMiXhnH655/wf4NxOv9Zfxyff4XGpmczeGMbav1o4jG5HN09p4SBKnlkG\ngFIqAAjw8/N70dRzEaLIcm7B/gWwZ57es2fAZ9B2bL5F3tw8xco/LzIvMIKs7DzGdffg1W4e0sJB\nlBqzDAAhLN7Fg/pDWuIjwPtJ6DsXqtbLd/ip6CSmrQ0mKDZZWjgIo5EAEKIkZdyAbTPg+H/A0RVG\n/Bea9sl3+L0tHL4a3pqBLaSFgzAOCQAhSoJSEPw7bJmk3+nT8TXoNgXKV8pnuGLNCb2FQ+JNaeEg\nTEMCQIjiSjwPG9+Cczugfht47neo1zLf4Wdut3D483wirVyqsWJMe5o3kBYOwvgkAIQoqtxs/Xm8\nf3wMNnbQ7xNo939gY7hoe28Lh9mDfXmmnbRwEKYjASBEUUQf1ou8caHgOVA/+Ds2yHf4ttBrvL8+\nhNikDJ5q68xkaeEgzIAEgBAPIzMZts+Eo8uhan145ifwHJDv8HtbOPz3nx1p7y4tHIR5kAAQojCU\ngtC1sPlduHkdOrwM3aeCfRWDw+9t4TClvydjOrtLCwdhViQAhHiQGxdh09twZivUbQHDf4EGbfId\nfuBcPNPXSgsHYf4kAITIT24OHPo37J4DaNBnNrT/J9ga/rW5u4WDS40K0sJBmD0JACEMiTkGG96A\nq0HQtB/0nwfVXAwOlRYOwlJJAAhxt8wU2PkhHF4ClevAP74Hrycgn5W5p2OSmLpGb+HQ2cOJDwY1\np7G0cBAWQgJAiL+EBcCmiZB6Bdq/CN2ngYPhBVrJGdl8GhjBj39epGZle74c3hp/aeEgLIwEgBDJ\nMfqBP2Ij1GkOw34AZz+DQ5VSrDt5mQ83hpJ48xajOroxoXdTqkoLB2GBJACE9crLhT+/0S/5qDzo\n9QE88irYGj6Yn4+/yfS1wew7G09LaeEgygAJAGGdLp/UV/JeOQkePfVe/dXdDA7Nysnlmz+i+HrX\nWextbZg1yIcRHRpiKy0chIWTABDWJSsNds2GPxdBxZrw1HfgMzjfIu+hqASmrgni3PWbDGhRjxkD\nvald1cHIkxaidEgACOsRsRk2vg0pMdB2DPR8HypUMzg08eYt5mwK49djMbjUqMCKMe3o2kzu6Rdl\ni9ECQNO0J4EBQFVgmVJqq7HeW1i5lMuweaJ+l08tLxi7FVw7GByqlOL347F8tDGU1MwcXunamHHd\nm1ChvNzTL8qeQgWApmnLgYFAnFKq+V2v9wUWALbAUqXU3Pz2oZRaC6zVNK068CkgASBKV16u3rRt\n+0zIy4Ye70HH18GuvMHh566nMXVNEIeiEmnbsDqzB/vSrK7hXj9ClAWFPQNYAXwNfP/XC5qm2QIL\ngV5ADHBE07T16GEw557txyql4m7/edrt7YQoPVeD9CJv7DFo1A0Gfg41Ghkcmpmdy6Ld51i0+xwO\n5WykT7+wGoUKAKXUHk3T3O55uT1wVikVBaBp2i/AIKXUHPSzhb/R9BUyc4HNSqnjxZm0EPm6dRN2\nz4WDC6FCdRiyFHyfyrfIe+BsPFPXBnM+/iaDWtVn2gBvalWRPv3COhSnBtAAiL7r6xjA8IVV3etA\nT8BR0zQPpdRiQ4M0TXsJeAnA1dW1GNMTVufMNtg4AZIuQZvnoedMqGi4935CWhYfbQxj9YlYGjpV\n5IcX2tOlSS0jT1gI0zJaEVgp9SXwZSHGLQGWAPj5+anSnpcoA1KvwpbJELIaajaF0ZvArbPBoXl5\nil+PRTN7Uzjpt3J4vbsH/5LGbcJKFScAYoG72yM6335NCOPIy4Nj3+lF3pxM6DYVOr8BdoYv4Zy5\nlsrUNcEcvpBIe7cazB7SHI/aUuQV1qs4AXAEaKJpmjv6gf8ZYESJzEqIB7kWChvGQ/Sf4NYFBs6H\nmh4Gh2Zm5/L1zrN8s+cclezt+GRoC55q6yxFXmH1Cnsb6M9AV6CmpmkxwAyl1DJN014DAtHv/Fmu\nlAopiUlpmuYP+Ht4GP6FFlYsOwP++AQOfAn2VeHJRdByeL5F3j2R15m+LpiLCekMadOAqf295GHs\nQtymKWW+l9n9/PzU0aNHTT0NYS7O7YQNb8KNC9DqWeg1Cyo5GRwal5rJhxvCWH/qMo1qVuLDJ5vT\nyaOmcecrhIlomnZMKWW4pe1dpBWEMH9p1yFwMgT9CjUaw6gAcH/M4NC8PMXPRy7x8eZwMrPzGN+z\nCS8/3liKvEIYIAEgzFdeHpz8EbZO1+/vf/xdeHQClDPcjC38agpTVgdx/FISHRs58eFgeTqXEAUx\nywCQGoDgegQEjIdLB6BhZxj4BdRqZnBoxq1cFuw4w9K9UVRxsOOzp1sypE0DeTqXEA9glgGglAoA\nAvz8/F409VyEkWVnwt7PYN8XUL4SPPG1fr3fxsbg8F0RcUxfG0zMjQz+4efM5H5eVK9kuNePEOLv\nzDIAhJWK+kMv8iaeA99/QJ/ZUNnw6txrKZl8EBDKxqArNK5ViVUvPUKHRoYLwkIIwyQAhOndTICt\n0+DUT1DdHUaugcbdDQ7NzVOs/PMi87ZEkJWbx1u9mvLS442wt5MirxAPyywDQGoAVkIpOPmTfvDP\nSoEub8Fj70C5CgaHh1xOZsqaYE5FJ/GoR00+fLI5bjUrGXnSQpQdZhkAUgOwAvFn9ZW8F/aCSwfw\nXwC1vQwOvZmVw/ztkSzff4HqFcux4JlWPNGyvhR5hSgmswwAUYblZMG++bD3U7CroLdwaDMq3yLv\n9tBrzFgfQmxSBsPbu/BuX0+qVZQirxAlQQJAGM+F/fqn/vhIaD4U+syBKnUMDr2SnMHM9aFsCblK\n0zqV+e3ljvi5GW7tLIQoGgkAUfrSE2Hbe3DiB6jmCs/+Bk16GRyam6f4/uAFPg2MIFcpJvZtxv89\n2ojydobPEIQQRScBIEqPUnr7hi2TIeOG3qr58UlQvqLB4UExyUxZE0RQbDKPN63FrEHNcXUyPFYI\nUXxmGQByF1AZkHBOfzpX1G5o4AfPr4O6zQ0OTcvK4bOtEfznwAWcKtvz9YjWDPCtJ0VeIUqZWQaA\n3AVkwXJe/USyAAAPrUlEQVRu6a2a98wD2/LQ/1PwGws2hu/TDwy5yvvrQ7iaksmzHVx5p48njhXK\nGXnSQlgnswwAYaEuHdL791wPA+9B0PdjqFrP4NDYpAxmrAthe9g1POtWYeGzbWjjWt3IExbCukkA\niOLLuKE/lvHYd+DoAsNXQbO+Bofm5Oax4sAFPt8WiVIwpb8nYzq7U85WirxCGJsEgCg6pfQHsW+e\nBOnx0PE16DoZ7A23YD4ZncSU1UGEXkmhu2dtPhjkg3N1KfIKYSoSAKJoblyAjW/B2e1QvzU8+yvU\nb2VwaEpmNp8FRvD9oYvUrmLPomfb0Ld5XSnyCmFiZhkAcheQGcvNhoMLYfdcvbDb92No/6LBIq9S\nis3BepH3eloWozq68VbvplRxkCKvEObALANA7gIyUzFHIeANuBYMzQZA/0/A0dng0OjEdN5bF8yu\niOv41K/Kt8/70dKlmpEnLIQoiFkGgDAzmcmwYxYcWQpV6sGwleA10ODQ7Nw8lu07z/ztkdhoGtMG\neDG6kxt2UuQVwuxIAIj8KQVh62Hzu5B6FTr8E7pPA/sqBocfu3iDqWuCCL+aSm/vOrz/hA/1qxlu\n7SyEMD0JAGFYUjRsehsit0BdX3hmJTRoa3Bocno2HweG8/PhS9St6sA3I9vSx6eukScshHhYEgDi\n73Jz4M/FsGs2oKD3R9DhZbC9/5+KUoqA01f4ICCUxJtZjO3szpu9mlLZXv5ZCWEJ5DdV/M/lE3qR\n98opaNIHBnyqd+804GLCTaatDWbvmXhaODuyYkw7mjdwNPKEhRDFIQEgICsVdn4Eh7+BSrXh6f/o\nrRwM3Kd/KyePb/dG8eWOM5SzteF9f29GdnTD1kbu6RfC0phlAMg6ACMK3wib3oGUy9DuBejxHjgY\n/iR/5EIiU1YHcSYujX7N6zLD34e6jg5GnrAQoqSYZQDIOgAjSI6FzRMhfAPU9tE/9bu0Mzg0Kf0W\nczeH88uRaBpUq8CyUX708DL8JC8hhOUwywAQpSgvFw5/Cztn6X/uORM6/gts71+dq5Ri7clYPtwQ\nRlJGNi891ojxPZtQsbz8sxGiLJDfZGty5ZRe5L18Ahr3gAGfQQ13g0PPx99k2tog9p9NoJVLNX4Y\n7It3/apGnrAQojRJAFiDrDTYPQcOLYKKTjB0mf5QdgNF3qycXL75I4qvd53F3s6GWU82Z0R7Vyny\nClEGSQCUdZGBetfO5GhoOxp6vg8VDD945VBUAlPWBBF1/SYDW9TjvYHe1K4qRV4hyioJgLIq5Qps\neRdC10EtTxgbCK6PGByaePMWszeF8duxGFxqVGDFmHZ0bVbbyBMWQhibBEBZk5cLR5fDjg8gJwu6\nT4dO48Cu/H1DlVL8fjyWjzaGkpqZw6tdG/N69yZUKG/4+b1CiLJFAqAsuRqsF3ljj4L74zDwC3Bq\nbHDouetpTF0TxKGoRNo2rM7swb40q2u4yZsQomySACgLbqXDHx/Dwa/1RVyDl0CLf+Rb5F28O4qF\nu85iX86G2YN9eaadCzZS5BXC6phlAMhK4IdwZjtsnABJF6H1c9BrFlSsYXDo3UVe/5b1mT7Qi9pV\npMgrhLUyywCQlcCFkHoNAidD8O/g1ARGbwS3Rw0OvXG7yPurFHmFEHcxywAQBcjLg+P/ge0zIDsD\nuk6GR98EO/v7hiqlWHMilg83hpGSkc0rXRszToq8QojbJAAsSVwYBIyH6EPg1kUv8tZsYnDo+fib\nTF0TxIFzCbR2rcacIb541pWVvEKI/5EAsATZGbBnHuxfoD+OcdC/odWIQq3k/fD2Sl4p8goh7iUB\nYO7O7dKLvIlR0HI49P4QKtU0OPTw+USmrAnibFwaA1rUY4as5BVCFEACwFzdjIfAKXB6FdRoDM+v\nh0aPGxyalH6LOZvCWXU0GufqFfhuTDu6SZFXCPEAEgDmRik48SNsm643cXtsInR5C8rd/0leKcW6\nk5eZtSGUpIxs/vl4I8b3aCpFXiFEoUgAmJPrkbBhPFzcD64dYeB8qO1pcOiFeP2ZvPvOxtPKpRo/\nDvHFq54UeYUQhScBYA6yM2HfF7DvcyhXAfy/hNYjwcbmvqG3cvJYsuccX+48i72ttGsWQhSdBICp\nnd8DG96EhLPg+zT0mQ2VDV+/P3Ihkcmrbxd5fevxnr83daTIK4QoIgkAU0lPhK3T4ORKqO4Gz60G\njx4Gh977TN7lo/3o7inP5BVCFI8EgLEpBad+ga1TITMZHp0Aj70D5SsaGKpYf0ov8t5Il2fyCiFK\nlhxJjCnhnF7kPb8HnNuD/3yo42Nw6MUEvci790w8LV2q8Z+xzfGp72jkCQshyjKzDIAy1w0055a+\ninfPPLBzgAGfQ9sx+RZ5v90bxZc7zlDO1oYPBvnwbIeGUuQVQpQ4swyAMtUN9OIBvX9PfAT4DIa+\nc6FKXYNDj17QV/JGXkujv29d3hvoQ11HKfIKIUqHWQZAmZCeqHfsPP49OLrCiF+haW+DQ5PTs5m7\nJZyfD1+iQbUKLBvlRw8vKfIKIUqXBEBJUwqCftN79acn6s/j7ToJylcyMFQRcPoKHwSEciP9Fi92\ncWd8z6ZUspf/LUKI0idHmpKUGAUb34JzO6F+G/3WznotDA69lJDOtHXB7Im8TktnR1aMaUfzBlLk\nFUIYjwRAScjNhgNfwh+fgE056DcP2r0ANvf35MnO1Yu8C7brRd73/b0Z2dFNirxCCKOTACiu6MMQ\n8AbEhYKXP/T7BKrWNzj02MVEpqwOJuJaKn196jLjCW/qOVYw8oSFEEInAVBUGUmwYyYc/U4/4D/z\nM3j2Nzg0OSObT7aEs/LPS9R3dODb5/3o5S1FXiGEaUkAPCylIGQNbJkEN6/DI69Ct8n6k7ruG6rY\ncPoKMwNCSbyZxQuPujOhlxR5hRDmQY5ED+PGRdj0NpzZCvVawohVUL+1waHRielMWxvMH5HX8W0g\nRV4hhPmRACiM3Gw49G/YNQc0G+gzB9q/BLb3//iyc/NYtu8887dHYqtpzPD35nkp8gohzJAEwIPE\nHNOLvNeCoGk/6D8PqrkYHHr80g2mrA4i/GoqfXzq8P4TPlLkFUKYLQmA/GSmwM5ZcPhbvXXDsB/B\ncyBo93+ST8n8X5G3blUHloxsS28fw+0ehBDCXEgA3EspCAuAzRMh9Sq0fxG6TweH+x+3qJRiU9BV\n3g8IISEtizGd3JnQuymVpcgrhLAAcqS6W1I0bHoHIjdDHV8YthKc2xocGp2YznvrgtkVcZ3mDaqy\nfFQ7fJ2lyCuEsBwSAAC5OXD4G9j5EaCg1yz99s58irzL953ni+2R2Gga0wd6M6pjQ+xs72/tLIQQ\n5kwC4PIJvch75RR49IIBn0H1hgaHnrh0g8m3i7w9verwwSAf6leTIq8QwjJZbwBkpcKu2fDnYqhU\nC55eAd5P5lvk/TQwgh8OXaROFQe+GdmWPlLkFUJYOOsMgPBN+rX+lFjwGws93oMK1e4bppRic/BV\n3l8fwvW0LEZ1dOOt3k2p4lDOBJMWQoiSZbQA0DTNC3gDqAnsUEotMtZ735FyWT/wh2+A2t7w9Hfg\n0t7g0Jgb6by3LoSd4XH41K/Kt8/70dLl/pAQQghLVagA0DRtOTAQiFNKNb/r9b7AAsAWWKqUmpvf\nPpRSYcDLmqbZAN8DxguAvFw4shR2zIK8bOgxAzq9Drb3f5LPyc3ju/0X+HxbJJoG0wZ4MbqTmxR5\nhRBlTmHPAFYAX6MfuAHQNM0WWAj0AmKAI5qmrUcPgzn3bD9WKRWnadoTwCvAD8Wcd+FdOa0XeS8f\nh8bd9SJvjUYGh56MTmLK6iBCr6TQ06s2Mwc1p4EUeYUQZVShAkAptUfTNLd7Xm4PnFVKRQFomvYL\nMEgpNQf9bMHQftYD6zVN2wj8VNRJF8qtm7B7Dhz8N1SsAUOWgu9TBou8qZnZfLY1kv8cvEDtKvYs\nfq4NfXzqohkYK4QQZUVxagANgOi7vo4BOuQ3WNO0rsAQwB7YVMC4l4CXAFxdXYs2s8it+qMZky9B\nm+eh50w9BO6hlCIw5Coz1ocQlypFXiGEdTFaEVgptRvYXYhxS4AlAH5+fqpIb3bkWyhXAcZshoad\nDA6JTcpgxrpgtofF4V2vKt+M9KOVFHmFEFakOAEQC9zdFtP59mum9+RisK8Mdvb3/VVObh4rDuhF\nXqVgan8vxnSWIq8QwvoUJwCOAE00TXNHP/A/A4wokVkVVyUngy+fjkli8uogQi6n0N2zNh8M8sG5\nekUjT04IIcxDYW8D/RnoCtTUNC0GmKGUWqZp2mtAIPqdP8uVUiElMSlN0/wBfw8Pj5LY3Z0i7/cH\nL1Czsj2Lnm1D3+ZS5BVCWDdNqaJdZjcGPz8/dfTo0WLtY8vtlbzXUjMZ+UhD3u7TjKpS5BVClGGa\nph1TSvk9aFyZbQVxOSmDGetD2BZ6Dc+6VVj0XBtau1Y39bSEEMJslMkAWL7vPJ9ujUApmNLfkzGd\n3SknRV4hhPgbswyA4tYAwq+m0MG9Bh8Mao5LDSnyCiGEIWWyBpCVk0t5Wxsp8gohrJJV1wDs7WxN\nPQUhhDB7cmFcCCGslASAEEJYKbMMAE3T/DVNW5KcnGzqqQghRJlllgGglApQSr3k6Oho6qkIIUSZ\nZZYBIIQQovRJAAghhJWSABBCCCtllusA/loJDKRomnamkJs5AqaqGhvjvUv6PUpif0XdR1G2e9ht\nagLxD/ke1syUvz9FYer5lvb7F3f/DQs1SilVJv4DlpTl9y7p9yiJ/RV1H0XZ7mG3AY6a6t+DJf5n\nyt8fS5xvab+/sb6/snQJKKCMv3dJv0dJ7K+o+yjKdqb8/2sNLO3na+r5lvb7G+X7M+teQEIUlaZp\nR1UheqEIYc3K0hmAEHdbYuoJCGHu5AxACCGslJwBCCGElZIAEEIIKyUBIIQQVkoCQFgVTdMaaZq2\nTNO030w9FyFMTQJAWAxN05ZrmhanaVrwPa/31TQtQtO0s5qmTSpoH0qpKKXUC6U7UyEsg1m2ghAi\nHyuAr4Hv/3pB0zRbYCHQC4gBjmiath6wBebcs/1YpVSccaYqhPmTABAWQym1R9M0t3tebg+cVUpF\nAWia9gswSCk1Bxho3BkKYVnkEpCwdA2A6Lu+jrn9mkGapjlpmrYYaK1p2uTSnpwQ5kzOAIRVUUol\nAC+beh5CmAM5AxCWLhZwuetr59uvCSEeQAJAWLojQBNN09w1TSsPPAOsN/GchLAIEgDCYmia9jNw\nEGimaVqMpmkvKKVygNeAQCAM+K9SKsSU8xTCUkgzOCGEsFJyBiCEEFZKAkAIIayUBIAQQlgpCQAh\nhLBSEgBCCGGlJACEEMJKSQAIIYSVkgAQQggrJQEghBBW6v8BrUmKQrRTWC8AAAAASUVORK5CYII=\n",
	"text/plain": [
	"<matplotlib.figure.Figure at 0x110ff9a90>"
	]
	},
	"metadata": {},
	"output_type": "display_data"
	}
	],
	"source": [
	"plt.loglog(h_values,err_values)\n",
	"plt.loglog(h_values,np.array(h_values)**2)\n",
	"plt.legend(['Errors','2nd-order'])"
	]
	}
	],
	"metadata": {
	"anaconda-cloud": {},
	"kernelspec": {
	"display_name": "Python 2",
	"language": "python",
	"name": "python2"
	}
	},
	"nbformat": 4,
	"nbformat_minor": 1
	}