{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "

# Worked Solutions 6: Implicit Methods

" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will consider the pair of (linear) first-order ordinary differential equations\n", "\\begin{eqnarray}\n", "\\frac{dx}{dt} &=& \\phantom{-}398 x + 798 y , \\\\\n", "\\frac{dy}{dt} &=& -399 x - 799 y , \\quad \\quad x(0) = 1, \\; \\, y(0) = 0 ,\n", "\\end{eqnarray}\n", "which define the functions $x(t)$ and $y(t)$. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "1(a). Exact solution
\n", "The equations can be written in matrix form, as\n", "$$\n", "\\frac{d\\mathbf{x}}{dt} = \\mathbf{A} \\mathbf{x}, \\quad \\text{where} \n", "\\quad \\mathbf{x} = \\begin{pmatrix}x(t) \\\\ y(t)\\end{pmatrix} , \\quad\n", "\\mathbf{A} = \\begin{pmatrix} 398 & 798 \\\\ -399 & -799 \\end{pmatrix}.\n", "$$\n", "\n", "The matrix $\\mathbf{A}$ has a pair of eigenvectors $\\mathbf{e}_i$ and eigenvalues $\\lambda_i$ such that $\\mathbf{A} \\mathbf{e}_i = \\lambda \\mathbf{e}_i$. In our case, $\\lambda_1 = -400, \\lambda_2 = -1$ and $\\mathbf{e}_1 = \\begin{pmatrix} -1 \\\\ 1 \\end{pmatrix}$, $\\mathbf{e}_2 = \\begin{pmatrix} 2 \\\\ -1 \\end{pmatrix}$. It is straightforward to show that the general solution of a 2D set of linear ODEs can be written in terms of the eigenvectors and eigenvalues as follows:\n", "$$\n", "\\mathbf{x}(t) = A_1 e^{\\lambda_1 t} \\mathbf{e}_1 + A_2 e^{\\lambda_2 t} \\mathbf{e}_2 ,\n", "$$\n", "where $A_1$, $A_2$ are constants. The initial condition $\\mathbf{x}(0) = \\begin{pmatrix} 1 \\\\ 0 \\end{pmatrix}$ implies that $A_1 = A_2 = 1$, which gives the exact solution $x(t) = 2 e^{-t} - e^{-400t}$, $y(t) = -e^{-t} + e^{-400t}$. \n", "\n", "(Alternatively, you may have substituted the exact solution back into the ODEs to verify that it is a valid solution.)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEKCAYAAAD0Luk/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYVNW97vHvDxFlEkdGmQdRJhVEQAwdkyg4RqMeYiQ3\nxhhjNDcmPsYcc86RnCTXmDPEGJNrSMxgEmOiiVESULhqe0RRGaRBGWSSeQggCsjQNL/7x6qyiqaL\n7qaqa++q/X6eZz817d61ej/d691rrb32NndHRESSqVnUBRARkegoBEREEkwhICKSYAoBEZEEUwiI\niCSYQkBEJMEUAiIFZmYHzKzXEf7saDNbVOgyieSiEJBYMLN3zOwDM3vfzHakHh9owu8bY2Zrmmjz\nDZ58Uzsw3H2Gu5/eNMUSOVTzqAsgkuLAJe7+QpG+z2hEZX0E224ozdaUSKklIHFSZ+VpZj81syey\nXt9nZtNTz483s8lmttnMtqaed85a9wQz+6WZrUt9/hczawVMATpntTo61vG9F5vZW6nP15jZ17M+\nu8nMlprZFjP7q5l1ylH2F8zs81mv/5eZvZR6/mLqd56f+o5rardQzKx/ahvvmtkCM7ss67NfmdmD\nZva31M/PNLOeDdnRImkKASkFdwADzeyzZnY+cAPw2dRnzYBfAl2BbsAHwE+yfvZ3QEvgdKA98EN3\n/wAYB6x397bufpy7b6zje38B3OTuxwEDgecBzOwC4P8AVwOdgNXAY434fRzA3cekXg9KleHx7M/N\nrDkwGXgGOAX438Dvzaxv1rb+CbgHOB5YDnyvEeUQUQhIrPzVzLaljnq3mdmNAO6+G5gA/BB4BLjN\n3TekPtvm7k+6+1533wXcC3wEIHV0fhFws7u/7+417v5SI8qzDxhgZm3d/T13n5d6/zrgYXevcvdq\n4J+BkWbW7Qh/71zdRyOB1u5+n7vvT3WV/Q34dNY6T7r7HHc/APweOPMIyyAJpRCQOLnC3U909xNS\njw+nP3D3WcAKQoWZPmLGzFqa2c9SA8vbgReB483MgFOBbe7+/hGW51PAJcCqVJfMuan3OwOrssq2\nC9gKdDnC78mlE1B78HpVre/JbsF8ALQpcBmkzCkEJE5yDqia2a1AC2A9cFfWR3cAfYFz3P14Uq2A\n1LbWACea2XF1bLLeAdnUEfYnCV0xT5EJn/VA96yytQZOAtbWsZldQKus14eMPRzGekI3V7ZuwLpG\nbEPksBQCEntm1g/4DvAZwljAN8xscOrjtsBu4H0zOxGYmP65VD//VOCnqQHk5qkxBYBNwEk5AgIz\nO9rMrjOz49y9BtgB1KQ+/gNwg5kNNrNjCOMDr7p7XaeczgOuSrVY+gA31vp8I5BrTsFrwAdm9o1U\n2SuAS1PfL1IQCgGJk8mps1zSy5/N7Cjgt8C97v6muy8D7gZ+a2ZHA/cTjrS3AK8QzvrJNgHYDywm\nVPxfBXD3JYTKdEVq/KGuI/QJwMpUN9MXCWMBuPtzwL8CfyEclfcExmf9XHYr44dANaGy/xVhoDrb\nROCRVBmuzv4gNd5wGXBx6vd7EJjg7kvr+B6RI2L53lTGzE4lDNZ1AA4AP3f3Qyb5pCb+jCM0jz+X\nNcgmIiIRKcRksf3A1919npm1AeaY2TR3X5xewczGAb3dvW9qcO0hYEQBvltERPKQd3eQu29MH9W7\n+05gEYeeJXEFobWAu78GtDOzDvl+t4iI5KegYwJm1oNwnvJrtT7qwsGnuq2j8KfTiYhIIxUsBFJd\nQU8AX021CEREJOYKcgG51PT2J4DfuvtTdayyjoPPdz6VHOc6m5nOeBARaSR3b8yFCz9UqJbAL4GF\n7v6jHJ8/TepaL2Y2Atju7ptybczd8X378Jdewv/7v/GXX8YPHAjvJ2i55557Ii9DHBbtB+0L7YvD\nL/nIuyVgZucRJvEsMLM3COcu302YUenuPsndp6SuyLiMcIroDfVu+K67YOpUOP98mDQJ9u6FL30J\nbrwRTjop32KLiAgFCAF3fxk4qgHr3daoDVdWwq9+BSNGgDvMng0/+Qn06QOf/CTceisMG3aEpRYR\nEYjrjOEPPoAlS+DM1AURzeCcc+DXv4alS6F/f7j6ajj3XHjkEdizJ9LiNpWKioqoixAL2g8Z2hcZ\n2heFkfeM4UIzM/cZM+D222HWrNwr1tTAlCmhdTB3Lnz+86G7qEePopVVRCQOzAyPeGC4sF57DYYP\nP/w6Rx0Fl10GzzwDM2aEMYOhQ+GKK2DaNDhwoDhlFREpYfEMgddfD109DdWvH/zwh7B6NVx6Kdx5\nZ+gy+tGPYPv2piuniEiJi2cINKQlUJfWreGmm2DePPjlL+HVV6FnT7j5Zpg/v/DlFBEpcfEcE2jX\nDrZtg2YFyKiNG+HnP4ef/SwEwq23wlVXQYsW+W9bRCQG8hkTiGcIjBoFL79c2A1XV8NTT4WB5CVL\nQovhi1+ELrqEkYiUtvIbGG7duvDbPProcFrpCy/A9OmwZQsMGgTXXBPmJMQsDEVEiiGeIdDUXTUD\nBoQWwTvvQEUFfPnLIRB++lN4/0jvSS4iUnriGQLHHFOc7znuuDBG8NZb8MAD8Pzz0L17uDTFq6+q\ndSAiZS/ZIZBmBhdcAE88AYsXh1NOJ0yAwYPhxz+Gd98tbnlERIpEIVBbhw7h4nVLloR5Bq+8Es4q\nmjABXnpJrQMRKSsKgVyaNQutgz/8AZYtg7PPDmcTnXEG/Nd/hYFlEZESF88QiNs5/CefDF/7Gixc\nGOYczJ8frmY6fjw895wuUSEiJSueIRCHlkBdzGD0aPjNb2DlyvD8a18LgfCd74TLVoiIlBCFwJE6\n4QS47TaoqoI//SnMTD7rLPjEJ+DRR2H37qhLKCJSL4VAvszCzW1+8hNYty6cXvqb34SZyF/6UrgO\nkgaTRSSm4hkCcRsTaKhjjw3jBM8+G1oIXbvC9dfDwIHwH/8RWgsiIjESzxAopZZALl27wre+BW+/\nDQ89BIsWwemnw+WXw5NPwr59UZdQREQh0OTM4Pzzw6Wt16wJVzC9/3449dQwqFxVFXUJRSTBFALF\n1KYNfO5z8OKLYRJa69ahZTB4cOguWrcu6hKKSMLEMwRKdUygMfr0ge9+N5xq+uMfhxnKgwaFs4se\neQR27Ii6hCKSAPEMgXJtCdSlWTMYMwZ+8YvQErjpJnj88TCm8JnPhHso798fdSlFpEwpBOKkZUu4\n9lqYPBmWLoURI+Df/i0Ewte/DnPn6nRTESkohUBcnXIKfOUr8Prr4aY3rVvDpz4VTjf9/vfDILOI\nSJ4UAqXgtNPCZSmWLw+nm65YAUOGwEc/CpMmwdatUZdQREpUPEMgCQPDR6JZs3C66aRJsH59aClM\nnw69esGll8Lvf68BZRFplHiGgFoC9Tv22DDn4PHHYe3aMFP50UfD/INrrw0T0vbsibqUIhJzBQkB\nM3vYzDaZ2fwcn48xs+1mNje1/MthN6gQaJy2bcPlKf7+99Bl9PGPh9tlduoU5iU8+6zOMBKROpkX\n4GwTMxsN7AQecffBdXw+BrjD3S9vwLbc334b+vbNu1yJt25daCn84Q/wzjtw9dWhxXDeeaFrSUTK\ngpnh7nYkP1uQmsDdZwD13Yi34QXUmEBhdOkCt98ermT6yivQuTPccgv06AF33hnOPNIppyKJVszD\nwZFmNs/M/m5mZxx2TXUHFV7v3uGCdm++CVOmhH08YUIIhDvugFdfVSCIJFBBuoMAzKw7MDlHd1Ab\n4IC7f2Bm44AfuXu/HNtx37Yt3LRFmpZ7CIXHHw/Lrl2hy+iaa+Dcc9VlJFIi8ukOKkoI1LHuSmCo\nu2+r4zO/5+674eijAaioqKCioqIgZZTDcIe33oInngiB8P77YXLaNdfAyJEKBJEYqayspLKy8sPX\n3/72t2MRAj0IITCojs86uPum1PPhwJ/cvUeO7bhXV0Pz5gUplxyhhQszLYTt2zOBMGqUAkEkZiJv\nCZjZo0AFcBKwCbgHaAG4u08ys1uBW4BqYDfwNXd/Lce2vFDBJAWyaFGmhbBlSwiEq64KE9cU1iKR\nizwECkkhEHNLloRAePJJWLUqzFS+8spwCeyWLaMunUgiKQQkGqtXw1//GgJh7ly48MIQCJdcAu3a\nRV06kcRQCEj0tmwJl8B+8slw57SRI0MgXHEFdOwYdelEyppCQOJl506YOjUEwtSpcMYZIRCuvDLM\nVxCRglIISHzt2wfPPx8C4amnoH370Dq47DIYNkxnGokUgEJASkNNTZiZPHkyPP00vPtuGFi+7LJw\n0btWraIuoUhJUghIaVq+PBMIs2fDRz4Cl18egqFz56hLJ1IyFAJS+rZvh2eeCYHwzDPhRjmXXRZC\n4cwzwY7o71skERQCUl6qq+Hll0MgTJ4cbo6T7ja64IJwQx0R+ZBCQMqXe5iglg6EqqoQBJdcAuPG\nhTupiSScQkCSY8uWcNrp1KnhjmlduoQwuPjicF2j1IUHRZJEISDJVFMTbowzZUoIhWXLwllG48aF\nRYPLkhAKARGAjRtD62DKFJg+Hbp1Cy2Eiy+GESN0sTspWwoBkdr27w9zEqZODaGwalVoJVx8MYwd\nq0tZSFlRCIjUZ/36cOrp1KmhldCrF1x0Ubjo3ahRuqWplDSFgEhjVFfDzJkhDKZNC/dLOP/8EAgX\nXgj9+2tegpQUhYBIPrZtg+eeC4EwbRocOJAJhI99DE4+OeoSihyWQkCkUNxh6dJMILz4IvTtmwmF\nUaOgRYuoSylyEIWASFPZty8MMKdDYfHicI2jCy8Md1NT15HEgEJApFi2bg2Xxp42LYwp7NsXZjBf\ncEHoOurePeoSSgIpBESi4A4rV4bxhOefD0ubNiEM0sHQvn3UpZQEUAiIxIE7vPVWJhRefBG6ds2E\nwpgxuveyNAmFgEgc7d8Pc+eGQHjuuTC2cMYZmVA47zxo2TLqUkoZUAiIlIK9e8P8hHQoVFWFW2yO\nGROWkSMVCnJEFAIipWjHjnDfhBdfDMv8+XDWWZlQGDUKWreOupRSAhQCIuVg1y545ZVMKLzxBgwa\nBBUVIRTOOw/ato26lBJDCgGRcrR7dxhHqKwMoTB7dhhTGDMmBMPo0RpoFkAhIJIMe/aE+yekWwqv\nvQb9+oXJa6NHh5aCro6aSAoBkSTatw9mzYIZM8Ly8stw0kkhENJLv36a0ZwACgERCRe+W7QoEwoz\nZsDOnaGFkA6Fs8/WtY/KUOQhYGYPA5cCm9x9cI51HgDGAbuAz7n7vBzrKQRECmXt2tBCSLcUli6F\noUMzoTBypMYVykAcQmA0sBN4pK4QMLNxwG3ufomZnQv8yN1H5NiWQkCkqbz/fhhsTrcUZs2C3r0z\nYwojR4brH6kLqaREHgKpQnQHJucIgYeAF9z9j6nXi4AKd99Ux7oKAZFiqa4Op6KmQ2HmzPD+yJHh\nvswjR4YJbZrEFmulEAKTgXvd/ZXU6/8HfMPd59axrkJAJCrusHp1CIOZM0Or4c034fTTQyCkw6Fn\nT7UWYiSfEGhe6MIUwsSJEz98XlFRQUVFRWRlEUkUs9Ad1L07jB8f3tu9O1wDaeZM+Mtf4M47w3WR\narcWNLu5aCorK6msrCzItqLqDloMjFF3kEgJcoc1aw5uLSxYAKedlmktnHsu9Omj1kKRxKU7qAch\nBAbV8dnFwK2pgeERwP0aGBYpI3v2ZFoLM2eGSW07d4YWwvDhYTnnHOjUKeqSlqXIQ8DMHgUqgJOA\nTcA9QAvA3X1Sap0HgbGEU0RvqGs8ILWeQkCkHGzcGM4+ev31zGPr1plQGD48nK563HFRl7TkRR4C\nhaQQEClT7rB8eSYQXn8d5s0L4w/ZrYXBg+GYY6IubUlRCIhIaaquDndjS4fCrFmwbBkMHBgC4Zxz\nQmuhf39oHsvzWGJBISAi5WPXrjC+kG4xzJ0L69eHFsLZZ4dQGDo0nLZ69NFRlzYWFAIiUt7eey9M\naps7F+bMCY+rV4cWw9ChmXAYMCCR10ZSCIhI8uzYEcYU0qEwZw6sXBnuuZAdDIMGlf0Yg0JARARC\nV1JV1cHBsGxZmMOQDoazzgrB0KZN1KUtGIWAiEguu3eH+zfPmROWqipYuBBOPRXOPDMsQ4aEx86d\nS3KCm0JARKQx9u+HJUtCd1L2AgeHwplnhlZEzAegFQIiIvlyhw0bQhhUVWWCYc2acCZSdqthyJBY\n3YdBISAi0lR27QrXRkqHQlVVeN2+fQiDwYPDGMOgQeF6SUcdVfQiKgRERIqppiYMOM+bFwIhvWzc\nGFoNgwaF01fT4dCpU5OONSgERETiYOfOMAM6OxgWLAj3f04HQnoZOBDati3I1yoERETiyh02bTo0\nGBYtCl1KtcOhX79GD0QrBERESk1NDaxYcWg4rF4NvXqF2c8DBoTJbwMGQN++OcNBISAiUi727Amn\nr771VmZZuDCcpZQjHKxFC4WAiEhZ2707Ew4LF2YCYu1abM8ehYCISCLt3o21aqUQEBFJqnzGBJoV\nujAiIlI6FAIiIgmmEBARSTCFgIhIgikEREQSTCEgIpJgCgERkQRTCIiIJJhCQEQkwRQCIiIJphAQ\nEUmwgoSAmY01s8Vm9raZ3VXH52PMbLuZzU0t/1KI7xURkfw0z3cDZtYMeBD4GLAemGVmT7n74lqr\n/o+7X57v94mISOEUoiUwHFjq7qvcvRp4DLiijvWa7i7LIiJyRAoRAl2ANVmv16beq22kmc0zs7+b\n2RkF+F4REclT3t1BDTQH6ObuH5jZOOCvQL9cK0+cOPHD5xUVFVRUVDR1+URESkZlZSWVlZUF2Vbe\nN5UxsxHARHcfm3r9TcDd/b7D/MxKYKi7b6vjM91URkSkEaK+qcwsoI+ZdTezFsB44OlaBeyQ9Xw4\nIXwOCQARESmuvLuD3L3GzG4DphFC5WF3X2RmN4ePfRJwtZndAlQDu4F/yvd7RUQkf7rHsIhIiYu6\nO0hEREqUQkBEJMEUAiIiCaYQEBFJMIWAiEiCKQRERBJMISAikmAKARGRBFMIiIgkmEJARCTBFAIi\nIgmmEBARSTCFgIhIgikEREQSTCEgIpJgCgERkQRTCIiIJJhCQEQkwRQCIiIJphAQEUkwhYCISIIp\nBEREEkwhICKSYAoBEZEEUwiIiCSYQkBEJMEUAiIiCaYQEBFJsIKEgJmNNbPFZva2md2VY50HzGyp\nmc0zszML8b0iIpKfvEPAzJoBDwIXAQOAT5tZ/1rrjAN6u3tf4GbgoXy/V0RE8leIlsBwYKm7r3L3\nauAx4Ipa61wBPALg7q8B7cysQwG+W0RE8lCIEOgCrMl6vTb13uHWWVfHOiIiUmTNoy5AXSZOnPjh\n84qKCioqKiIri4hI3FRWVlJZWVmQbZm757cBsxHARHcfm3r9TcDd/b6sdR4CXnD3P6ZeLwbGuPum\nOrbn+ZZJRCRJzAx3tyP52UJ0B80C+phZdzNrAYwHnq61ztPAZ+HD0NheVwCIiEhx5d0d5O41ZnYb\nMI0QKg+7+yIzuzl87JPcfYqZXWxmy4BdwA35fq+IiOQv7+6gQlN3kIhI40TdHSQiIiVKISAikmAK\nARGRBFMIiIgkmEJARCTBFAIiIgmmEBARSTCFgIhIgikEREQSTCEgIpJgCgERkQRTCIiIJJhCQEQk\nwRQCIiIJphAQEUkwhYCISIIpBEREEkwhICKSYAoBEZEEUwiIiCSYQkBEJMEUAiIiCdY86gKIiEhu\nBw7Au+/CP/4BmzfX/ZgPc/fClLRAzMzjViYRkUJxh507YdOmsGzenFn+8Y9DK/mtW6FNG2jfHk45\n5dDHU06B664z3N2OpDwKARGRPLnDe+9lKvb6FjPo0CEs7duHx9oVe/r5ySdDixaH/34zhYCISEEd\nOADbttVdiW/efOjrFi0yFXt9S5s2hS2rQkBEpIH27IENG3IvGzeGin3Llkw3TEMq9pYto/udFAIi\nkmju8P77dVfotd/74APo2BE6dap76dgx001TXzdMXEQWAmZ2AvBHoDvwDnCtu79Xx3rvAO8BB4Bq\ndx9+mG0qBEQECJX71q2wfv3hj943bIBmzXJX7NnLiSeGPvlyEmUI3AdsdfcfmNldwAnu/s061lsB\nDHX3dxuwTYWASALs3Rsq93Xrci/r14duls6d66/c27aN+jeKTpQhsBgY4+6bzKwjUOnu/etYbyUw\nzN23NmCbfuCAl11SiySFezivva5Kfe3azPP33gtdL126HH5p1Srq3yj+ogyBbe5+Yq7XWe+vALYD\nNcAkd//5Ybbpe/d6yfTFiSRJTU3oZ1+9+uAKvfbSokWowE89NXfl3r596MKR/OUTAvXOGDaz6UCH\n7LcAB/6ljtVzJcp57r7BzE4BppvZInefkes79+4tnQEZkXKR7n9fsyZU8mvWZJb06w0bQp96t24H\nV/CDBx9cwRf6FEhpOvWGgLt/ItdnZrbJzDpkdQdtzrGNDanHf5jZk8BwIGcIfOc7Ez9sAlZUVFBR\nUVFfMUWkHjt21F3Bp99buxaOPRa6dj14GTQo87xLFzjmmKh/E6msrKSysrIg2yrEwPA2d78v18Cw\nmbUCmrn7TjNrDUwDvu3u03Js09etczp3PuJiiSROTU3ohlm1Ct555+CKPv28uvrgyr1bt0MrfB3B\nl6YoxwROBP4EdAVWEU4R3W5mnYCfu/ulZtYTeJLQVdQc+L27f/8w2/QVK5yePY+4WCJlp7o6HKm/\n806mos9+vn59uMRA9+6ZpXaFf8IJ5XdqpARlN1ls0SKn/yHnGImUr717wxF7dsWe/bhpUziTpkeP\nUMHXfuzaVd00SdakA8NR2Ls36hKIFNbevaEyX7Gi7op+69Yw0JpdsX/845nnXbrA0UdH+itImYpl\nCOzbF3UJRBrHPVxEbMWKsKxcmXm+YkU4ku/aFXr2DBV7jx5w6aWZSr9TJzjqqIh/CUmkWIaAWgIS\nR7t3h6P27Mo9u9Jv2TJU8r16hWXUKLj++vD81FOheSz/2yTpYvlnqRCQKLiHiVDLl9ddyW/dGo7c\n05V8r17wkY+Ex5494bjjov4NRBovliGg7iBpKu5hwtPSpbBsWViyn7dqBb17h6VXL7jgAvjCF8Lz\nzp01w1XKTyxDQC0ByceBA+GUydqV/NKl4Si/bVvo0wf69g2P11wTnvfuDe3aRV16keJSCEhJOnAg\nTI6qfSSfrujbtTu4oh8/Pjz27q1uG5FsCgGJtW3bYMmSg5e33w799Mcfn6nk+/SB667LPNfMV5GG\niWUIaEwgWfbtC5V67cp+yZJwQHDaaZll/Hjo1y9U/q1bR11ykdIXyxBQS6D8uIdz5bOP5tPPV68O\np1CmK/rhw2HChPC8Y0dd6kCkKSkEpKCqq0Of/MKFsGgRLF6cqfSbNz/4qH706HBU37u3LnkgEhWF\ngByRPXtC5Z6u7NOPy5eHo/ozzoDTTw+nWN5yS6j0Tzop6lKLSG2xDAGNCcTHjh2hcs+u6BcuDGfm\n9OqVqeyvvjo89usXZs6KSGmIZQioJVB8W7ceWtEvXBjOzjnttExlf8MN4bF3b13QTKQcKAQSZseO\nULkvWABvvplZdu8OFX26sv/EJ8Jj9+6aJStSzmIZAuoOyt/evaHPPl3Jpyv9zZtD5T5wYFjGjYMB\nA8KlinUWjkjyxDIE1BJouJqacHGz2kf2K1aEi5oNHBjuEfv5z4fnvXrpksUikqEQKCGbN8P8+VBV\nlTm6X7Qo3FZw0KBQyV9+Odx9N/Tvr9MuRaR+CoEY2r8/dOVUVR287NkDgwfDkCHhWvU33xz68HUt\nHBE5UrEMgSSNCbz77qGV/aJFoY9+yJCwfPnL4bFbN/Xbi0hhxTIEyrElcOBAuMpl7Qp/+/bQlTNk\nCJx7Ltx0U3itC6CJSDEoBJrA/v3haH7OHJg9OzwuWBD67tNH9zfcEB579tQpmCISnViGwHvvRV2C\nhqupCdfFmT07s1RVhe6cYcPCcu21ocLXDUtEJG7M3aMuw0HMzFu1crZvj9+M1HSXTrqynzMH3ngD\n2rfPVPjDhsFZZ6nCF5HiMTPc/YhGDGMZAv37O489Fo6eo+Iezr/PPsKfOxdOOCFU9EOHhsezz4YT\nT4yunCIi+YRALLuDhg0LlW6xQsA9XNM+u8KfMyfctCR9dH/XXaHiP/nk4pRJRKQYYh0CN95Y+G27\nhytgpiv6dKXfvDmcc0747ttvDxV+x46F/34RkTiJbQj87neF2dbGjQcf4c+eHfr200f4t9wSHjt3\nLsz3iYiUkrzGBMzsamAicDpwjrvPzbHeWOB+oBnwsLvfd5ht+s6dzimnhIlUjbn0webNBx/dz5kT\nro6ZPWg7dCh07apJVyJSPvIZE8j3DPUFwJXAi7lWMLNmwIPARcAA4NNm1v9wG23dGvr0CdfHyWXb\nNpg+He69Fz71qXDJ43794D//E3buhOuvh5degi1b4Nln4XvfgyuvLK1Zt5WVlVEXIRa0HzK0LzK0\nLwojrxBw9yXuvhQ4XLU6HFjq7qvcvRp4DLiivm0PHw733Reuhrl2LTz3HPzgB+Gc+169oEePULFv\n3QrXXBM+37YtPN53X3ivZ8/SqfDroj/yQPshQ/siQ/uiMIoxJtAFWJP1ei0hGA7r+98PR/VDh4bb\nFfbuHZ5ffjn8+7+Ho37NtBURyU+9IWBm04EO2W8BDnzL3Sc3VcFOPjkEwb33lvbRvIhInBVkspiZ\nvQDcUdfAsJmNACa6+9jU628Cnmtw2MziNXtNRKQExGGyWK4CzAL6mFl3YAMwHvh0ro0c6S8iIiKN\nl1evupl90szWACOAv5nZ1NT7nczsbwDuXgPcBkwD3gIec/dF+RVbREQKIXbXDhIRkeKJ5PwaMxtr\nZovN7G0zuyvHOg+Y2VIzm2dmZxa7jMVS374ws+vMrCq1zDCzQVGUsxga8neRWu8cM6s2s6uKWb5i\nauD/SIWZvWFmb6bG5cpSA/5HjjOzp1N1xQIz+1wExSwKM3vYzDaZ2fzDrNO4utPdi7oQgmcZ0B04\nGpgH9K+1zjjg76nn5wKvFrucMdoXI4B2qedjk7wvstZ7DvgbcFXU5Y7w76IdoXu1S+r1yVGXO8J9\n8c/Aven6lpbzAAACgUlEQVT9AGwFmkdd9ibaH6OBM4H5OT5vdN0ZRUugIZPHrgAeAXD314B2ZtaB\n8lPvvnD3V909fZudVwnzLspRQycVfgV4AthczMIVWUP2xXXAn919HYC7bylyGYulIfvCgbap522B\nre6+v4hlLBp3nwG8e5hVGl13RhECdU0eq12x1V5nXR3rlIOG7ItsXwCmNmmJolPvvjCzzsAn3f3/\ncvhZ6qWuIX8X/YATzewFM5tlZhOKVrriasi+eBA4w8zWA1XAV4tUtjhqdN0Zy6uIyqHM7KPADYTm\nYFLdD2T3CZdzENSnOXA2cAHQGphpZjPdfVm0xYrERcAb7n6BmfUGppvZYHffGXXBSkEUIbAO6Jb1\n+tTUe7XX6VrPOuWgIfsCMxsMTALGuvvhmoKlrCH7YhjwmJkZoe93nJlVu/vTRSpjsTRkX6wFtrj7\nHmCPmf0PMITQf15OGrIvbgDuBXD35Wa2EugPzC5KCeOl0XVnFN1BH04eM7MWhMljtf+JnwY+Cx/O\nON7u7puKW8yiqHdfmFk34M/ABHdfHkEZi6XefeHuvVJLT8K4wJfLMACgYf8jTwGjzewoM2tFGAQs\nx/k3DdkXq4CPA6T6v/sBK4payuIycreCG113Fr0l4O41ZpaePJa+v8AiM7s5fOyT3H2KmV1sZsuA\nXYSkLzsN2RfAvwInAj9NHQFXu3u9F+ArNQ3cFwf9SNELWSQN/B9ZbGbPAvOBGmCSuy+MsNhNooF/\nF98Ffp112uQ33H1bREVuUmb2KFABnGRmq4F7gBbkUXdqspiISILpYswiIgmmEBARSTCFgIhIgikE\nREQSTCEgIpJgCgERkQRTCIiIJJhCQEQkwf4/l2U+zVes2FkAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "ts = np.linspace(0,1,200)\n", "xexact = 2*np.exp(-ts) - np.exp(-400*ts)\n", "yexact = -np.exp(-ts) + np.exp(-400*ts)\n", "texact = ts\n", "plt.plot(texact, xexact, 'r-', texact, yexact, 'b-');\n", "# Note the rapid change at the beginning of the interval\n", "plt.title(\"Exact solution\");" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEPCAYAAAC5sYRSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsfXm8HUWV/7fv8rbkJSHsBEEEBUVFUNxwiYqKOuo44+7M\nuOtPFHR0RnAFxAV13MVxVFAHmUFARUVcQEBElGUA2QQhhB0SXhKSl7d239u/P849t07XPVVd3ffm\nJWCfz+d9br/u01XV1d3fc+p7TlVHaZqikkoqqaSSh5bUtnYDKqmkkkoqGbxU4F5JJZVU8hCUCtwr\nqaSSSh6CUoF7JZVUUslDUCpwr6SSSip5CEoF7pVUUkklD0Fp9FtAFEXDAC4CMNQp78w0TY/rt9xK\nKqmkkkrKSzSIPPcoisbSNJ2OoqgO4A8AjkzT9LK+C66kkkoqqaSUDISWSdN0urM5DPLeq5lRlVRS\nSSVbUQYC7lEU1aIougrAfQDOTdP08kGUW0kllVRSSTkZlOfeTtP0QAC7A3hKFEWPGUS5lVRSSSWV\nlJO+A6pS0jTdFEXRBQAOA3CDPBZFUUXVVFJJJZWUkDRNo6Ln9O25R1G0QxRFSzvbowCeD+BGTTfd\neWekN96INE2R/tu/IT3+eNpevx7p5z5H238Df8ccc8xWb8O28lf1RdUXVV/4/8rKIGiZXQFcEEXR\n1QAuBfDrNE3PUTXjmP54O0loe/Vq4HvfM3pHHAHcf/8AmlZJJZVU8rcpfdMyaZpeC+CgIGUJ6EmS\nBXreBoBf/xp429uAHXcETj4ZmJsD3vUuOjY3BwwP99vsSiqppJKHtCzsDFUXoNvgLv+/7Tb6A4C7\n7wYOPtjo3X470G5v6VYPXFauXLm1m7DNSNUXRqq+MFL1Rf+ysODuAnQN3KWHz9ubNgEPPGD0XvMa\n4MorafvOO4EbVap/m5PqwTVS9YWRqi+MVH3Rv2w7njsDeJ6eNAIzM8DsLG3/8IfAf/6n0TvxxC13\nHZVUUkkl27gs/MJhGqBLMOdjRekbub1mDXDssUbv+OOBVasGfimVVFJJJduqLDy4hwRU+wV3W+83\nvzHgftZZwDHHmGMbNgzmuiqppJJKtiHZup57KGizQbDpG9cxn97dd1MgFgCmpoB99zV6114LbN7c\n/zVWUkkllWxl2XbA3QfarGfTN/J/17Z2jMuenc0GaI8+Gjj/fNpetQr4xS/6u9ZKKqmkkq0k2wbn\nzttpSn+tlg7UDPQ8a6tf+kYCPQDMz9MfAPzhD8B//zdtpynw+tcP5vorqaSSShZAth7nbgMuH5PH\nXXqtln5MbrORKKJn1y/1/vd/jd473wmc05mEOzdHaZiVVFJJJduQbBu0jA/wbQ/f1ndtu+py1aG1\nyVXe/fcD69bR9q9/DRx5pLm+H/zAjCwqqaSSSraSbD1wt+kW/g0FfO2YNhLg/wdhBLQ2zc6aXHsA\neMtbKFALAN/5DtE7lVRSSSULLNsO5w7kA74Gtq5tXxkhoO0yND49ux0XXQRcdx1tX3YZ8MpXGr0b\nMisiV1JJJZUMVLYtzt0FvBrI2oHXMkYgjwLylafptVrULq28iYnsSpcrV9I+APjd7ygNs5JKKqlk\nQLJtcO4+bz7UCJShb/LqCqVvXHq+iVpzc/QHAKefDpx7Lm2vXg285z1Gb9MmVFJJJZUUlW0D3F05\n7bzPR9+EcumhHnlRY5E3EtDaoB3j7bvuAv70J6P35CfTxCsAuOkm4NZbUUkllVSSJ9sGuJfh3EMp\nHe1YiEeutSWPm88zVtJzt49p5wDkuXOA9uSTgf/5H9qemAA++lGjV2XoVFJJJUK2brbMoCmYorRM\nKGhr7cwbbfjKS9PwlS9denfcAfzoR0bvkEOI0gHI07/jDlRSSSV/u7LtB1T7pWrabfoLNQKDSMEM\nNQJ5IwtfXfbs2okJs5TCyScD3/wmbW/eDBx3HCqppJK/Ldm2aBmXV90P5x56rlZ/aKC0aEDVV2/e\n8scu+kb+L5dRuOceAnuWww4D/vIX2l6zxnj7lVRSyUNKti1wl0Dej4df1OvPA+NQ+iZPr2jAV7sG\nF33j8uptvTVrzOzaU08Fvvxl2p6dBT70IaPH6/xUUkklD0rZNiYxDRrIXWDo09OOhdAyRbJliurZ\n/TQoPc3DX7cO+Pa3jd7LX04TrwDgvvvMZKxKKqnkQSFbh3PniT6DBPJ+9/GvNDgheoMoz+WRt9u9\nk6K0sn1l+IyAz8PfsMF8yOSnPwW++EWj97a3Gb3JSbOgWiWVVLLNyNbx3Ntts23/Ft03KK/fN7s0\nlJsf1LIHRde+yaurrBHQjk1NAaecYvTe/GbgvPNo+847zXr4lVRSyVaVhQX34eFwekLjyF0ep70v\n79xBePhl9IqCduhIwFdXnl4ZYyH1JifpDwAuvNBk6aQp8OxnG97+jjtMvn4llVSyxWVhwX101IAD\nAz2QBeiQfaFg7Dq3LNfvMhBl9fLAOM/DHzRoM2WmtUlup6kZfdnH5PZFFxnK5qijgJ//nLZvvjnL\n7zPvX0kllQxMth64j45mc9D5WJIAUaSDTgi4l9HzTaiyy5Nt4nN9elp5efSNRhHJungf9x8vouYy\nEL4Pl9htkR8uCTVMIXpyLZ1rr6UPlbM84hFm2eQ//hG47TZUUkkl/cnCgvvYmAGNZhNoNAywjo1l\ngV8CZhSFg03oSpEuPdurt/W02aWaHv/mlcfHQkcdml6esXC1SbvuoqDdrx5A6ZkM7t/6FvCb39D2\ntdcChx9u9K6+GpVUUkmYLDy4M7g0m/RnA7rcbrWAWg0YGjJgV6/7JyIN2uMsMmEpFGS1tueBsG9f\nqBHQ6t9SoB2q51uK4e67geuvR1de8AJKywSA004zwds0NcsnV1JJJQC2Ji3DnrsEdH7JbQ9fGgHb\nQGjgFeKZDooPLwraeXGALWkYil63rzxXf/v08q7DpwcQN8////73wJVX0vZVVwF///dG76tfBWZm\naLuajFXJ36hsHc9dgrYEdBeAM4jb9A1vA/RrZ+O4smr4/zIebGjGS9l9eYapn3O1+QUa+GvXGmIk\nfSOqPCMQAu6ue7t5czYT5z/+g5ZOBoATTiCwB4j6+fGPUUklfwuy8J47A0Cj4fbINSPA/2v0DdB7\nbGQk3zPth3sO8bTLUDV55ZXdF8fZoKrWB676B0HfuIxAngFzgbvL6NrH1q83k7Fuuim7TPLLXka5\n+QAFcW+5BZVU8lCRre+5+0Cb9Vz0DW8zbzsykj3mMgKcjRMCcmXTHrXVKH11hQL0IPYNmrbqR0+r\nV2ufb9Qh25rXJql3003m04ennEIfNAfIILzmNUbv6quz51VSyYNAtg1w17JlbA8/j76p1bK0jK0n\nAUMek4aDjwP5nmRRkPMBahEevmyQtYiHP8jr7kdPAr1M0fSVZ59XRm9iArjkEqP3zncCV1xB27/8\nJXDSSeaYDPhWUsk2JFsnoMpcOqdC2oDu8/A1cA/VA7Jgz/Wy3tCQH2TzwDBEb1AB1bx9/WTaaPqD\nCKjmxQ60Y0X1tGNavXl1Sb25OfP/9debQO799wPPfa7R+8QnaIIWQHGAe+9FJZVsLdm6ee5FQTtJ\nsgbCNha2HntmGrjLumygZz3J29tUTig33w9Aa3UNYp8P8JOk+AdOyugVDajmGSsXaLvoJU3PVZ6r\nDDkxCwB++1tg1SraPv104GMfM+e/4x1G7847s+dVUskWkK1Hy9gBVelB+/LhfUZAUizapCig1yj4\n6BuNt2cKKBS88oyAzSXnef/2vjLcfGj2jVZH2WyZIgHVEA/fLoPPycsIcgVeXdtl9WZnDYBPTQHf\n+57RO+IIM1Hr//4P+NrXzLG//hWVVDII6RvcoyjaPYqi86Mouj6KomujKDrSqWwHSl3ZMpKi8dE3\ntp7k5jVg5sCrC7RD6RupxwaobABSy2BxldEP4OedWxTw88AuVC/P+BXN5vFx80UCry4jkOfhu/pR\nljc7a2bk3nAD8Lvf0XaaAvvvb67hK18xa+rPzBjKp5JKAmQQnnsC4P1pmu4P4GkA3h1F0X6qpisL\nhoFV49xdKZNF9fjFswO0rqwal1dv69l1DQ/7PVMJAPa21HcB5ZYOshY5V9MvSrfktclXbz9GRWtv\nGSMQx9nF1lzGx7e+j9333KaLLwZuvJG2zz0XeP/7TZte+Uqz4NrNN1czdCvpkb7BPU3T+9I0vbqz\nvRnAXwCsUJVDJzGFePhF9PKMgAu0NW9dbtt6tkGwy2M6JxSUNOD3gVdRfj0P8PNAmwGrjBEIAePQ\nALarD3x6Wvn2OdoKmVrb5fcJ8vpba2+InvxqFkArbE5P0/YJJ5jJWdddl12P51e/Mu2r5G9KBsq5\nR1H0cABPAHCpqiA9aF+ee2g+fN5aNZp3nmcENM69CDefl6Uj2xEKcmWAPHRf6Iihn33aMY2T1/RC\ny2Nj4zNWcZxdSTO0TS4j4Ou7MoZJa3tI0FgC/913Zz+J+M//bNbj+cxnzIdWJieBM8/MllfJQ0oa\ngyooiqLFAM4E8N6OB98jx559NgWXTjsNKx94ACs10G63dTDmB31sjL73aVM5rlUm7e1+UyY1bn7J\nEjc3v3Gj0XMZAbm2fZLoq2C6JkXZ3DD/+kCzH8Pg8+ZDA6p5oG3XEbJMMuCeNKb1Wa3mbx//akZY\nawtnVpUJBhcxYEDvktChRuCuu2jhPYBW3PzMZ4jeAYAnPYly+HfbDbjgArruZzyD6lm7Fth5Z1Sy\nMHLhhRfiwgsv7LucgYB7FEUNELCfkqbpT116x77xjfQ9zgceoHzhdps8iXbbAFqr1evhu7Jgli7t\nNQLSWExP94K0zfVrgMtAKlMhfTNeNSqGjQBzoS5Dwts8xLaP1etZ0HAN4eV2mYBqHhiGnstgnGcg\nuI9DQFtSP5o3HWrA5D4Gdx8Y5/HxoWDsK09rX+goJk8vry7pGKxfD2zaROD+61/Tc/yMZ1D2zqtf\nDfz5z6R39NHABz4A7LgjLdWcJMCKFahkcLJy5UqsXLmy+/9xxx1XqpxB0TInA7ghTdOveLUY+O69\nF3jCE4C99qLhYRQBRx4JnHMOrfB3661mBcCQlElNL4+bt7l+G+g1D98+B3BTNi4jII/J8kJoHu3r\nVT6QDdWzPWNNj7++pIGcXZcNxq418Ps1JFtqX17fucC9DH3julc+vdD2aeeFePhye2rKOB4AOWe8\nHs9JJ5k0zvXrsytznnMOTeQCsl/uqmTBZBCpkIcAeAOA50ZRdFUURVdGUXSYqjwyQr8rVwIPfziw\n557AK15BYHbSScBjHkMPwVln0e8HP0iLOW3YQA+ZDXxaKqQG6L4Fy/JSK13GQnqoLjC2177RDAlv\ny6G27f1zm0ZG3F+v4nXuAQPGocFOG6iLgkhZUMoD8jyDo7Vv0Pv43hX18ENBO3QkUGbEoJ3nMwJF\njYVcgnndOrNEAwAce6zx9k8+2SzYNj8PfOpTRu/uu7MjiEoGJoPIlvlDmqb1NE2fkKbpgWmaHpSm\n6a9U5TVr6HfPPQ14zszQ76JFwHbbAY98JPDZzxJYvexl5NVffz3wtreRB/HLX5LnPztL59Xrpgyb\nvuEHe2goq+cC7eFhGq5z2XaqJutxvQye2qxZn55dr6/tRfUYjHxccT98+JbS6ycOEGosiurFsT4P\nwXc9IZz7ltCT18bnaJ9rDAV3zdvX+i6kvA0bzMqc69cDX/iC0XvLW0ye/wUXAF/+Mm2nKb3rlZSW\nhZ2het55BE7z8wZkp6d7QavZJEDee2/goIOAl7yEbjoP7667DvjBD2hfFJHHUK/3Bl75IeXytbok\n4LJXL42ATbFw+/IoG02vCLj72qcZHHsGLa+TwwaH2xcyuzYUtDUj0U95rnL7MQL9AL7r3JCRjc/T\nzjOm2ojK1495I4sQMLZHbKFGoF+9mRkzoeuvfzXr9kxNZWmef/kXWt4BoNU8TzvNHLv1VlTSKwsL\n7uefT96xBNnp6d6JRXxMeuQcAHvnO4GnPIXonBe/mMD9v/6LJnz89rfECU5PG3CTIKuBtgbGsk0S\nPPkF4DbNzlL9cjVK13o3thFwrXzpAm3ZPq3toXp2po/sJ97nAhsbeH0g4jrHLk+e4+LmbfDy1e+r\nI3RfaB2h4JlXl3Z8EOXlGSbpkQ8atAetNzFBXj9AzhyncaYp8OhHmxHWpz4FXHSROYc/xQhkYwd/\nA7Kw4D4xQWAiwZOBXnqkEow1b5pv+sEHA+PjwDHHAM97HvC4x1EmzuQkTeS47DL6HNuqVVljYXu/\nNmjbbfIZAR+A+zzyIh6+bF+/eq6JWjalZHP4QC837wKlUL08I5AHbKEAPUjD4DI+vm/B5tUfet1l\nYyF5RsDX9q0J7iGzepMky/3/+c/A6tW0/Yc/0Fe5WB75SAPwp5xilnWenaVEjoeYLCy4M7hKUNR4\ncA08NUql0TBgPDoK7LMPDeVGRoD//m9g993p5v/mN5Rv/p73UADnzjuJA8wDbRcw223P09PA0+fh\nF53QpcUE8soDssZNGghbL457P2HoAkjtmEvPBuhQD59/84zFljpX27elv3IVquert6zxscFY8/41\nbt6lx0kBWnsHTfXJtk9MGAro/PPpIywAgfy//qvRO/hgchIBSgvl9X3a7QfV17oWHtyHhrKArnm/\nLnDX6Btbj4/tsAOweDFx9kcdBeyyC/Da11L9N9xAa2/fcQdw6qkUoJ2YKE7f5HHzrg+N5IGx9MJ9\nX6VyLWzmMjgSzBmMXJSSbykG28MvC+4uvbKcez/8ehEjkAfCIde9pfTyRkqh98UH2q4yOIAbUper\nva6228a/DLjb94nLs5d2uPVWyvkHgF/8wgR8r7wSeP3rjd7LX26+5HXppSZbKE3JidzKsrDgPjxs\naBnN+7U5aOlB+zz8vCwYDtDusQcFaV/1Klpre+lSSsmcmAC++10awm3cCPzkJ2SlXZOdfJSSBtrz\n8/Tglwmoarn8Gjcfun5OqIdfdp0daSy0zxnavyFgYwOKtnSAXYaPRgkFfNe5/Xj4WrkhwDeIfrSP\nubxku94ioN1vm0L07JGA/HWV58scCo0/yGWcAfL8162j7bPOAs4+m7b/+lfgMJEN/uY3G7C/9lrg\n8svNMc4g3AKysOC+aFFvoFQLXmqA7uPm5TE7yGlz+LzND+trXkOg//a3U7B2bIxuwO23A9/4BoH+\nvfeSlz8/38vN+4BUix1IT9vW07x1X5ZOHn2jGYEQD18DbTma0LYB/Zi8B4MCgCLLJIcagUGfW4ab\nL9o/oUZAoz149Ba6tEMcL8xHXPrVG5SHn2csfMemp+mdZ7n0UrO+z09/Sh9yAWg5iEMOMXrveY/J\n/LnppuynHkvIwNaWCZKxsXzOvdXSue9Qz72Ih2+D0YoVwPbbk6WdnaWOv+suSuH84Q+Jb/v610lv\nZoYA38eD5xkB3g5drtg1OgF6s3mk3vLlOrjnefjMO/pGCb7F0TQ+n8uT7WVqx/ciai9bntdmnyv1\n7LVl8kDbNxLol+bR2qeVEaqneam+Nmn7Go3efdJ7d53r4s99XHo/YOzTs4E+pDw+L+/Z0+ridXvy\n9Li86WlK+WS59FLCm0c8gmb43nwz8PSno6wsvOdue7Kh2TKufHg7f11u256xLwtGGwksXw486lHA\nsmXAV79KnX7oobRQ2OQkfVHnnntoudVLLiEwzIsdaCDro5R89E3od2e1rBqp56N5eG6Bxs3LbcB9\nTJbHei4DZi+xoFE7oS9s0X390C0+0PIZizygknpFQLuIwSlaXt41avqhoO0zAnlt9umF0jf2yMo1\nKsprUwjN4yvPpopKyMKC++LFvWDnAm1Nz6Y2JPBr21yezKrhbQ2ANA9f0iitFrDrrsBznkNpVV//\nOlE6++1HK+etWkVDqzvuIA7uzjsJnHxZNdLDl23SKBqfYQoNvPpWyHQBs299n1BwL8Lh2zSPvC4f\nUNiemoubd3l0ocDn877LgKZrxU8NnIq0yQfQeXVoYOgDcPnr+9RhqBHQ9HxtyQPPvGfAVa8PcH3X\n6ALtMuWVlIUF9/FxP91ie5QSFLcELRPH2ZUfQ+gbCZ4zMzQaeepTgec+lwD/e98jIwZQzu1115ml\nEy6/nAIr8/M6aHOwRnrJLi5dA8U8vbLGIuRrWIC/TWWNgCsYbOtJYwYY0OTj9m8IaPZrBOz1fUI9\n2H4AP3Tk4GtTv15/0SUbiuq52ly2PNlXoeWxziD1fMdKyMJy7osW0Y33gXGrZbzrzZt76ZY82iNU\nb26OeFf5WTzp4fuoEt8oYaedaN9znkP5sr//Pc2k/fzniU+7+Waib+65h3LtL7mE+H2XAZPgZi8/\nMDpqlhSWvLqmJ5c/lv3iGsW46Bt5TBo9IPu/bRBtI7D99m5w53VItDbJbU5X08rgbVeqpraee5L4\n12Lh35BF2UJBxAVYRainImAcOuooy7kX2afVH9qPGgjz/0XLK0Ib+cBYPlMhoG3rudpeUrYe567R\nLT4wtukbSWf4PHyNSw8xAjZ94+PmNYBkvaEhWvt62TLgTW+iNMwnPhH4p3+il+VPf6IUqssuo5Xz\n0pQMgLYKZh43b4OxppeXh59H39hUib1gmetLWS4joF2HfPh9HL6LvrHPkWvHa6maXBfgBm0fTeHz\nJF2cbVmQ63c0oYFiPwCdR7P49LRyi4Js3jUWrTfP+MnnpMwHUzS9vDTTkrKwnvvixZRH7lt+wM6W\nkcsPsHGwQYwDmTbgajTK8uX54B66cFie3vCwbgTSFHja0+j/t76V8mKvvZbWuL/mGloNb/Vqoow2\nbKDzb76ZvFCNm09TPUtHA/dFi3QjwAC5ebP7emW9POVb1luEm9eAn42D9KCXLTN6MoA8PExlMRjn\nTeiyDVgUZQ1THGcXWyubox8KDnkAHarn2tcvtWJ77nxP7DaFArlWhzYC2lJGIPS+5PWjb0QlPXfW\nC2k7sxUuI1BSFp6W4Zffplv4BZTgnselS+DXtm1w8qVWujx8Dah8ehJIOYCstV2C4tAQxSMe/3jK\nxPnoR2na8623AvvvT4uinXgiTXgYHwe++U3qpwceIJCNonwwzvPI7WvMo2/kyMI2oloZrrXo7W1e\nnrnVcnv4Wttds39dI7BazWy3Woaik3ryuuQLqX0GMZRG0cCjH70iuep5gC/LaDZ79Wxapl8DEtoW\n/vUZkND+drXdpRdaHj9TXDY/67wKq+8a2TFy6ZWUhQd39gj45ZRgPj9PDxC/kCGTnULoFj5v40Y9\ntdL2ZKenjSGSs0u1kYAGMlp5ISMBeY1c5yteQWD+iU/Qehi3305GI4poXez162kdfF4r5/bb/cse\nhILx8DDdi7k5/zUODdEDrBnVZcsM11/EkPioJxdo58ULbD0Gd9bz0XA+qkiWrdE89sscAiIu6iWk\nPBdtZJ+bV4ZW3iA49zIevtYm33WXAe2QEQM7nzYdw3ra3Ilms/dzklpbZGKHrVdSto7nDpiXjrfr\ndeq0ep1eEhuAeZu9Om1NeN5etozO8elps0a5rqkpKsOnpxkSGzxtI2XrtdvZ9mkpkyMjBmTHxugb\nl4ceSqmW3/gGfeVm40ZqL0Be/8aNRD/Nz1OwdvNmQ6NwH7pGEwwsklJyGQvZF671eNasMX1hX7+s\nd/163eCEGoHttus9R3r4tlGt1Yov3gbQr/x2r0YHsZ4sj/e5QDvUkyw7YvDpa/sWIqCaB8y+MrRz\ntespYiRdcwhC+0xOYuLfKArrC26Hi7cvIQvPubvAnQGdO6jZpM62veQoov/l15JsD3/HHUnPBjHN\nw5d6vjVttPL4HG1iVcjiaKxnc/Pa6EQLGjca1KbttqNlFH70I/LmzzjDlBXHtMb15CT1y113EbD5\nPPw4dlNKEph9VJb0eKXe7GxvltL4OM0TsI2ADbKjo9Qe2S9a210UlWyTpGW08rSRgB00tikgNo6A\nG/gl6LMe/4aCXBmu2M49zyvDB+6h3ncZL13bl/dFKV/92nFXGTb/76srr7/51wb30FHZgw7cXZ47\nP0gMlHyc98lz+NfFx/O2racZAann8uo1PT62aRN51pr3r5Xnmqhlr7kjy9MCzzaQshHkY+PjBOSH\nHEKzZ085Bfjc5+hBGxqiEcP73mcCtWvWkHfPHL7LqEiwcxmzIhPENIB0gbbt7WtlyHMmJnr14phG\nOEli4hRzc34jIK9D4+bzllPWypOZPi7gZ2DI49J9RsAGyBA9/uVnXl5PWc+9CJCHGovQ6+agbaj3\n3U/sQkuvDfHcpeOgGYGSsvXAXQNyBgd7Xx64S17YBe7a9vi4aYsLPBmY7fL4HJ6YxTdNBna5DM7X\nt41AHmhr4OkDd1tPesn1OgH+y15GX6r54heB73+fjFO9Trpf+AKtcnfddTQhK0koe4c/YegCd1fb\nNSMQkj5ql+Gig1xlaBPObD323Dk7SDNg0ut2jRh89I2LUpL0jQv4NT0b+CVQ+ECpCNdv77OBUu4L\n9dJDjUDRtuRN8ioK0Nr1FDFMMqAqj4Vco6/PSsrC0zJ2QJW3+Vd68fY+W1969XI4L/W1ICqD7/Ll\nRk/z6m3P3QZP28BwzECe55oUpU2Yylv5cnQ0a3BYT16/a8kGG4zZuO2xBwH5BRfQ7NrjjiMjMDxM\nHs8XvkCzam+/ndbWWbeOlleYnaXyhoZ665VgnBencBkLvn4JkEuW+I2Atm1fv03L2CMkW4/Bd9Om\n3vbZeqOjRH8Beptc3LwvR9+lZwe1JWD41uMJoUdcoD3ogKq2T5uZq7VZK8PnpZc1AnnAK49pAVX2\n3POMle/aSsq2Qcvkee424DPXzPy771xX3jxvs54G2nn0jQbgHDuw0zOZvtFAlrfZEGn0zcwMgZsG\nkLI/XWvlS1B0UUCjo1TWTjtRHv5llwHHHw98+9vAnnsSoF93Ha2QuWoVrZK5YgXVefHFph0aUNn9\n7vJ+uX1y1OEDba2MLRFQdY0EbNDesEGf0GWPLOQsXKaKND1eksJlIGQ9Uo+P5U3assHL9R3bfsFd\n481DQNtVnlaGS6+I9x0KvBrg9xNQfdB77nng3i9Vo5Uny5DctE3fSA/fBvpFi4yei5u3y7ONwPh4\nL8jadWlGWjLZAAAgAElEQVTer22YGg1qK18Teyxam0KXYrDBM0l6DdMee9A1/OEPwFe+Ahx9NHDg\ngQQ+d91FQdxVqygff9ddib+//HLyePPmK/g8fJtSsfU04F+0qPd6bb28gGrRkYDdf/U6edQu+mbt\nWnpWNECXerwkrMuAsd7mzUbPFci1YwIS+HmZBl+wTwuoaqmBGhhqx0LB2AfuRY1F6L6iRsB1PaEB\nVV8dJWXrLD8A6B65C6BDgNy1TzumBW0lb2/PjM3z3G09V0DVBbIaReMD45Dvzrq8c61NTCnJl53L\n0ECWwWPpUlo0bZddgH//d5ph+453AE9+Mj3Uf/kLUTif+ASlbN53H33bNkmoXP74iWvdHh+Ay2Na\nQNW+Dk0vJKvGBl/Zn7Z3XmTJBvmC+9bPkS+7zwi4yvPFBFwGbGQkC9ry+eFz5W9ZMA7R8xmLIuXZ\noxOXIQsB3tCRQD960giUlG2Xcx+E5+465tNjILDBko/LzBeN35Z6Li/Z5tIl4DKvLsFOMxayPBlQ\nnZ4m0HUBpJbuGboUA2fSuAKqQ0O0INjoKHn5H/kI8O53A69+NeWx//KX9Ds5aWbZNpv0ZZr168nz\nt/vZBaS2wRkdpdGCBDsbtFlv8+Zezt1lwHxGwAfMIeUBbqAuq6e1SepJ0LYNouwL7nNZXrvtDvJK\nrn9r0C0+42OX4Vot1GVIQgA6jvsLqPrqKClbN1vGB+55XHoRcOcgZ56eZgR43Qduk4ti4IAnn7tp\nkw6es7PZeIEN2kuWuDl325DwiyXbG7L8sV3vhg1Zj1QzbrzNL7tPzwaMNCVaZ+edgX/8R1pW4S1v\nIc7+vPOIxpmeBk47jdIym03g2GPJmPzlL/Qy8iQu2zjOzZkUTz6WR6Ns2EAvopyMpIE2B1FDZ9Bq\nMQb+X9bFej5uXgNtl6ct1+NJEhN4tvV4aQdfnILbLkc7sr1cHgc+JQVmz3CWbde+tqWBnH2sjLHo\nx8PXAL+fgCpLFVB1ALkWUC2yTwZeXUZFa5PWFslvy4lU09P0UrGey9O2PfzJSTdoa9y8Dcz8Umvx\nAnt1S7s8e339PM/dB2Ia9aSBrCxveJgogNe9jrz6I46g/PTzzwcOOsh8TPjOO8nbXrWKPPyzzqL/\nkySbZmq3TwuoSs5dW/RM6q1b10vlaB75hg0GBG2jJ4F0bi5rfDRu3s6+cXna9jVK0HbFAXz9ZJfH\nZckyJLjnGTqbKuL/tboksGmZPj4wLurhFwV3F/BqgB8aUC1KB5WUheXcmXIAdHDvl18P3aeNDjT9\nUCPAD6i8Do2317JvfKAtXyBtJU3Nwx8Uh69lt8h1dmR5Wl0ajaIZMBu0RkYoe+TAA4niec97gNe+\nltbGP/JIAuXJSfqA8J//TCOBOAZ+8APgyispuMt8vgZA2iJqLj2fMdOOFSnPXnpC0+O+zfO0NZDN\no7I0mievPHtbM2AhVJF2D1hPli9jQLwAF+vxbx5A2npF6Btt1UpNz1duWT15PSVlYT33kNRFLdi5\npQHf5dX7wD1vJKB96JuBVNI3oZOY7DJ47Rt+CF16LmMR4uHbAVXbO3eNLGwA0spzAZrm4ScJgeEB\nB9BL/rKX0czaW26hhdU+8AFg333pA8Nr11L65tQU8KEPmUDxZZdl65JL/oYAmn29moHglS9di63Z\n91RLVR0bM7EDVz+59FxGgA1dXqZPEc491CAWpZ5cRoCXqtb0ZBtkumcofaPp+YyA61g/M1R9+0rK\nwoI7oINnUcAfNLgXAW2tjLzyJG0yO2voG5eH32wSQPCSoZqBkC+WpIqY6/cZCxvQ82iZ7bbLNxYu\nzl0rTwM03mbqTvPwNWDhZQQOPZT67P77gX/9V+ClL6VA7vnnE+CfcQb1+bvfTSC8ZIlZre/OO3uN\nm80fA72g7RqduK7RpnZcGUG2Hj8L9ijBtcSCrcfLM/vq4glYNljK+yjbwFSW6z7yNi/ZwDER+7mw\n1+NZvDjcCGijCXskwN4/U2AauEpAdQG9ra8d8wVU88rV2sTvdknZOuAeRQa0eB//hgRUJfCG6uWB\nMS8M5muLC/j7oXk0gJRenZwU5QJtLX+ft3kBMH5Y+IXkXHkXMMtjO+2kGxUN0HxA5TI4mlfrWojM\nVZ4dAGy3gb32ovXw99gDeMlLKE3zq18FvvtdAjyeSHT00eQFj48T0E9NEeUzMWEC3CFt54lqWvwh\nj6LSvG7Zt3Z5chVMl7F0GViXhy9X5rT7nZ8zn1Gx65XGwqc3P0/3QR6zjUWeh28/c3Z5dtl2edL4\nMAUUQvNwv0hwt/V959r77Dp5TaASsnXAXQM5ID+DxrcvD/DlvlDAd7UzBLR9ei4j0G5nX+ixMXPc\ntbY9B1Q1kJ2eNssf+8DYB9oyYJmXWsmrW9ovkKtejau19WTbeY35+Xk/4Dab+gqezSawww7UJ8uX\nA29+My278MUv0gJrSUIzbi+/HLjiCpqNe/HFwM9/Ti/Z179OI4MlS8z6+Zqn7Vr+2O4zl0Hk/nNR\nVK7yXODJ5dl1uXL0XeW5AqouD5r7KE9PTtRyGaaxMXJIADp3fFx/zkZHqb+4PFeMxTWHwNazjYB8\nvpkCYgDWjIFGy2jr4vgMSUlZ2IAqEAZyrn2DAGiZVZNnVMqCe56H7zNgrrbYk6xs4NOCrbZHrnHk\neeBu0ze87aIYGNzqdfpjDloL0HJdnKLH2Uc+0LbBs0hcQdPjto+O0mSsgw+mNh9zDE3Qes1raILW\nyAhRRpOTBPxHH01pm5/9LBmByy6jGAAHfF1AKu+Ba4KYdg9CyrO9Vxu0ZF0uGkXWq5WXN4dA25Z6\ndraM5PPbbernELol1CMvEsj1XQfgHhUlCT3D9sQv1sv7iIu9LX/7CKouvOdehMvm3y3FuZfx8LU2\nlTUCZY2FDNbyQyCByvb4bDC2AYMnPtVqWTC2DQRv12rm60uu8prNXkrJNgKaYfJ5q7J8nhDnWxzN\nN1GL22a3b3aWjvF6860W8LCHGV5/40aie/bdF/jwh4E3vIEWXHvgAVptc9Mm4B/+gYzV+vXADTdQ\nSuXtt5OHqn2cxQa+vFGMC4x5HRrJM8vytkRAVQug+4yKSy9vZGFz8xrw26CdJNmcf1d5Gmhz2SMj\ndB8ZtH3ZVuyNSz0ZuOe67JngXJ4rGFxSti4tI7+6xMc0fpu5LAl2GkDLD33YdZUxFj5PW9av1atd\nj63n8tbt63W1k0E2ikw/assfc5aSRu3w9g476Px+yFr52rbsO18cwJVJEwLaNoDbxkJrE+tJcA/J\nWNIMTqtF3v3ZZ9P3b5/4RPowyjHH0Feydt/dfA3rxz8mz/7662nVzQceAE46ieqbmDAA4Io/aMAc\nQmWFeuQykCl1ZBk2uMs2uUYCGs2zfLneBnl/bDDOo1HyPHIbwF0BWi2YXquZlEhXvRLc7REE68nR\nia2ntV1OBCshAwH3KIpOAvB3ANakafp4r7IEKP7fFRTlIKd9nguMGcTyjICtp5XXaPQGWRmgbYOk\nlaeVKwElRC/UCDGvZ+vxQ8dt13Lli4J2v3rT0yb7poinLcHEB8Y27cEpo7YeB/Rdbc+jvCRg2DRP\nq0Ujocc8hkDz6quJ1z/6aOCZz6S6Tz+dgH9qij6Gvn49nffZz9IaPL/9LdFA7TYZAmlUXJ6xvL+8\nJHNeQJV58TwaZWYmy7m76KBQo5Kn5wqoSg86zyN31cXeM2BA2zWr19d2u8+k5y6Bukj6qNb2kjIo\nzv27AF4YpClBC3CDsWYE8kCuiF4euMtfl/Fx6fkomyJB47J1uMqTgdfQD5xoQMp6LpD16dnpnjz6\ncNUrgdQHxraxcFFUtrHwGT1fva4X3NazPcp6ndbF33574O/+Dth7b+Dww4nX32MP4J//mfjbyUma\noXvXXZTSee+99HGVn/2MKJ7rrqNg9+Sk3nZf9o3PMy4D2i6O3B5Z+Lj0UCPQbhvaTPN+h4eza8KH\nTlrTQNvXF66EAbuuMnrAQMB9IJ57mqYXR1G0Z1iNAtwAPxj79LSgqHbMp5cHuL429RMH8J3rap+r\nPElb5Z3ruh5+6HifvS6ODYJcrivYKsHddUwzAjZoj43lg3ZeGqc8R+rZtAzPnpbzBOwX0l4jyF5v\nXtPjz/jlBVSnpsj4PuUpdN3PfS4d+/nPKTf/E58Anvc8YPVqAvo//Ykydz76URO8/sxnqI5rrqHf\nqaks3cB1cdvz1uOxwThN9ZiABMX77/cbFQaxTZvCgI/TM/MmiPFyFD4DJid05YG7do2ugCrr8ShA\n9pk2Wc5Fm8mspFbLBJdLytbl3O3/88B90EHMsmAYUi+XVRTI88rT9GxaxqfnyuWX58pJV1qWDutL\nusD2ul1r4LuMha0XkgWTF1D1zQ2wA6qcz67lnmt18WJrcmSh0Tf8GT/Na8wLqLJxm5+n7cc8ho6v\nXUtfzHrTm2jC1m23AZdcQkHeK64wSy1v2ED/z81RGuf8PP1deqmhGG0w4vbxfgmQ7bYfqCSQ2ouy\nufRGRw0o+oCZy9Cei+2370271MB4dLR3QpcrKMuzf/lZlSMhec7kpB5XkF732Bi1W9PTRha2XknZ\nuqmQ9v8SbCQw8v8aeGnglAdyvnp9gKqVJ9sbCsb9APkg9pU1apylw/v65eZDPXyfZxwaUPWtx+Mz\nAq56Q9bUt0HGlT6aB2i+gGqSUAxjhx1oHf3ly4G3vx34+7+nUcB3vkMe4OMeR576pk1E78QxfST9\nzDPpAyynn07XdMkl2bxzzau1+9kG0iJZOs2mSYPle+/qC1e/axSQvQ6SlmEkDXMofeNKrdSMtK99\nWnmuuQslZUE992OPPZaGlOvWYeWFF2LlypW9gK6Bi/2/vc28rU+vqEfu8qqLGAsZeB0UD6/tcwVU\n8+ooQ0fJXxeHb3vkoUbAnjXLLzvTCBo3L5dHyANj+cLIgKoG2i4jsHx5vocfys2HLtkQwqW7jMDO\nO5PH/YQn0Ghq9Wri8N/3PgrennUWvZMrVlDZv/0t6WzeTEDfagH/9V/kyS5ZQrOVkyR/tq42d8EF\n7j6DyNu2XogR0OglzTOWZbioIq3tvlTIkIBqkvQut5AkuPD883Fhu00U14knoqwMEtyjzp9Tjj32\nWOqAD30IePjDaeeSJTTbDHADiu9YET35W8TDD6WDfF58GSD3Aa7U18A91DCEgrs2QtKONRpmpi3v\nc3nGMzM0fNf0XMskNxpZI2AP0+3USs2Ll5y7T883EgiZ+NVsmvX7XeX5vFpfm2xv2mcEXCOGxYvp\nb/lymqx18sm08ubvfkfe+8qVwHvfS9/T/c1vKGXz17+m+/CqV1FZ4+NkHKamaGbvmjVUp/2heBtk\nbeopz1hoMRYXyLoMrM+DlmXwGvVFA6r8cW6px2sf2c+IJ6C68qlPxcrFi+m7xW9+M4779rdRRgYC\n7lEU/Q+AlQC2j6LoDgDHpGn6XVW5VjPADpB3wOC+yy46kNr/u4DH/r8IeGmg5atrkNx8ESDXjEVo\nQDWkXtav18PomyK0FXtbEgR32MHou2gP28MPpW8kLz49Tdws68n0VVfWj2YsbD2Xh89r5Wu8un1d\nsgzNW81bcdNnBPLWtOEyhoboOWJ6ZGiI/ppN4DnPobX063XK5nnDG+hLWqedRqO3FSvoPb7pJpq0\ndfXVRPO0WsBRR1GMYHqa+nxqikYFHGh29bsLjLWVSbW+CPXwtWN5RsA1yhpUQNU2AiVlIOCepunr\nS5/MwA7QUJFl553JY2CRHLy9LY2ADTA+j7MMGIfq+YyAdkwDw1Aj4JvklWcsXO2zKaVB9g/fPw7c\n8j4teMsvNOcl2yBWJKC6yy5ZcOdjrq9m+UBRA20JLJK+sUFLu0YfsLi4eV+gmfV4Eg6DtgtI7bri\nuPcaeUTG92rpUvL8n/IUWozt6KOB//gP+rLWs59NSzEfdhhl/KQpTeCanAQ+/WnKgrnsMqJ5pqaI\n+tmwgVI/x8ep/NCAqovmcVEvUo/TcWXWkz1iaDSMoZdlcG48G2IJ7vYkJvsZ0e6jbdj5PpaUQdIy\ng5WddgJOOcX8/6pXGY/fBSj2/7ZXycf5N8TDD6Vv7GOhYCi5ZJ+er/6ynLsL8EONinasXz1+meQS\nCzyrlfVcgVL+5J6Lm5cvjLb8AJfHXrdrJOACXO2Yy/gw8DOlEsfZ9DdbzzYkNqft0pPGiOsKXd3T\n9n7lwnZ2umeaZmdejo3R//vsQ+/t3nvT30030UfSjz+eju2wAwV4eaG3a66hNM+JCfpQS6tFhmD5\ncqrz3HOpfM78yTN0si9cs3/lNbmMhQRj/vyiK2jMRmB6mp7ZvPRR7dvEsrySsu2Cuy3HHGO2X/Qi\n03E+cLe3eXYpkO/hh4KSy/t26bnqDeXDNb1Qzr0o4G8NcLf7UbsOHxjJ/vRlt7iWH5BBUxfw+Txy\nV71ssOTaMps2kffHevx82itauoK3miFxgfb0dLYuCdo+sJPlyRRZ38QvH8jK0cnICPDoR9MI4LDD\naGbu055G7TzxRAr4HnccsN9+1HftNn1zd80a4FvfojJqNeBLX6LUz1/+koze7Cxw442G97dn62pU\nVih9MzZGawVpRsAOqPIqlsuW6eWxweH0TJcRKCkPHnCXcuCBZnuPPegTayyjo/TQAG5Asf/3gbGr\njEFx7j49H32jgXueXqgB0a7RRW/ZAGy3eVA0WEh5cdybX6+BkRZQ9YGRRsVMT9NzVq/3grZG32ig\nHRJXYCOggacGxnlGIIR64nXpbQPm+26A1HMZHx+XrsU6dtuN/t97b1qn57bbaELX175GqZ6bNtFH\nWF7xCuA//5Puwz33EPh+8pMU6P3kJ2nUlySEExMTtMAbT3qamNC9brt9/MEUe3Qi+88VUM2b1Ztn\nVErKgxPcpdRqFPBhOfVUCvAAlCK33Xa07QP3PNAOoW9cej5ALcJlu8ptNPSFxrT6QwC/iDdvc/NF\n+qfsSEm2M+/a5OcHmc/XaBktQMtrGjFoa0ZgfLx3UTYJdvIcbpMN/Ha9Uk+WNzamjybyMn24HfKY\nzDbSPNnxcXOOBD77ewPatqtebrsr1TAvrsD9JNtUq9F9euYzidp51rPo/1/8giidT36SAsA330wB\n3/32o4yeW26hUcL991NweH6eYgUMury085//TF51FOkGTMuWsUcnoQHVJOmd+MV1lZQHP7jb8rCH\nme03vtF8pmrJEkrtYikCIj4wdOkVNQIaaNvH7Nml/OsD9zKA77u2PG8+pB/LjJTy2hlqLDjGkaaG\nm/fNwnXRIz5PW6M27PJcFIhLj4PB7EHys5C36JkGnj4PX4Inz0C2OXdu09AQbcuAtK9epp7yjICP\n+3YFVHnbHnUsXUp/Y2Pk7Z95Jn1Ufd06Avq3vAX44Afp9//+j4A/iqifzz6bjMDkJPDKV1IdX/2q\nCcIuWmQooDjOzvgtE1C1P+fI/VlSHnrgLoWXwgWoE3/2M3Psla+k6dyADiLaapRFPO1+PVOtPPsD\n4/wbki1T1CMP0Qvd57vuvJFSyPLH3Dcu46ftY6Nv67Va2fqnpnTw9BkB2yAsX67raXRQnh5nvbjW\nyndN6JKgnbfYmlav7blPT/u/d6vRYRs29FIgPt7e5+HnceShAVXOxNprL6JxNm6kGb7nnAMccQRl\n+Nx2G/DylwP/9m/Ai19MX+dav56AP0loFvC999K9ueQS2vf1r5NxXLSI6mm3KZjMcYu8a5T9VFIe\n2uDuk3e9y2y/8IX0rU2AvJHR0WyQ0gcULirGBXIuD78fDr9sQNVVfwjgu7x539r2dps1MLbLDV0m\n2W5nCGWTp2cf03LP+bgrnVLq+cAzT08bMWh6/L1bH3gW8fBtcNdA2zU3wAbZvCwdG/hDZrIW9fBd\n/eIDXJ4DMDwMHHQQBW/Hxyl1+9e/Bj7+cZrZOz1tJn49//nAr35FlM+qVfQcf+5zZBBaLfpw++Qk\nBYbXrCHjMDND9fDooQL3Acjuu9MfQMO6a64xx5Yv13n7omDsAhGfh++jaDRwz+O0tTbZ69O7QFtL\n2XS1xS6Dz+1nRBAK2qHnujx3rY9DDHKel8x68hN8ttctZ2q7uHQfN59Ht8i6li7VQdYVeA0NqGog\ny7NQQwOvbEhk8LIIuNttt3lwBnCXUdECqr50VL6u4WFDCx1yCC3LPDREvP6ll1IM4Kc/JW9+//1p\nRc9nPpMYhakpmgAWx5QZuHYt0UZnn42ysvALhz0YZJ99zPb3v09WmPc/+cm0rXmoGvDkBTt5W/sI\niLZtA6UL+LUPiAwiXhCqlwfCdr39gLZPL5SO4oCd1o+u2bp5bWdvUO5z0SMhNI9GB7mMAH9yj9vu\n88hds1xD8uFtvtwFsvZ2HjevGSlZl2y7b5mCEArEvkaN+w6djBaiNzJC6ZGtFrEGK1ZQBuDhh5Mj\n+dnP0icd3/xmoohKSuW55wl7xQB9PPngg83+L30pO3NyaMhs8wuo0Qratv2/C3C1MqQBsRdR04xA\nPxSQb3RQZnQySNAuAuShnnvoPu1YXj9KMHJx7howy3TPUCMg0ylD8+Z58hi31Q6o+gAyD/g07pv5\naQ2YZT/ZbbfXrLevUYK2vSZ86GjClSIr29tq9bZdPguyTfb3bu16h4aoPzheU0IqcO9H5HIJxxxD\n38YEKDNnyRJzzAZtDQR9xzQ9F7g3Gm69UFAKBe2yIFeGPvLVUcbr94G7b1Tk8vC1Y1r7QvvHxR/P\nzOiLrbkAl/VCUyZtLl0GaGVAVdJLPnCfmtJ5eh/dEmIsJPCPjPTONQhd0ydvAhZfv2/UweVJcOfy\n5D32zf71xR9KSgXug5JHPtJsr1wJPP3p5v8Xv7i/lS8laGh62nLHvvK03zI0igu8XHp5wCuXDnDp\n5bUzlJsPCagW2acd0/ovRM913c2mASne58u+ca3HM+iAKnvd3L4Qmseetm9/LpC3bdCenKQlC5rN\n7FLQdl/keeQ++kbqcYDaBlwNjDXaSD6Xdvsc5aXTM4g623++oxPvKyEV574lhL/xyHLGGeb/17zG\nTLrilfdYfC+9D7Q1TzOvPNa3zxsUfZMHzHZdRcDQvu68NoXq5YGra1+jER5w1o5pBsmn5yqP0zgZ\n+GQ7XYFSXpKY6UOpZwNk2SwYezuEc282TbzDtX6OK7gc4pFrQK9RRVr7PKOO+yYa3e07p5bj9qkd\ngEYDs/VF+P7lj+ke+9T5T8OaOs3A/dHEs3HiLS8Emk2sj7bH/l9+W1fviLNfgLJSgftCy6GHmrTL\npz+dUqFYVqwwHr6PvrEBpowRAPRc/jKgnQeoIeXlefhlAdplrDQ9rY5+6vXRNyH9WFZPq1dSKmU+\nrOILqLq8X5dRYUMoQdu3FEORJRtswNXaJMqYmjT9ctfmZdgwTyOQyfoynLNq3+45n7ziMEzXFgPN\nJr5772H40Qb6xu0t2AfPOevIbhmHnvVu3NzeG2g2cfLav8O37vk7oNnEbdFeOPaPL+iWd8ZfHovV\nNdK7ZWYFbtq8G9BoYLK2FGs3j3X14lZ5iK7AfWtKs0lTolmuvprWHAeAZzwDeNKTjJ4E7VBu3gZ0\nWy80l991bEtx82V5+Dxuvqjeltin9d0g9Ir2N/cBpzfycXuNnDzQlno2aLsAl+tygbHUs+MAApg3\nbq4jadA6Upvq2+G69bt19b65+oVI67T9nTtegEvbTwLqdVyJg/COq9/VLePRF3wDG0FpoZ+54w04\ndeYVQKOBq9In4Pgb/qGr942bD8Ua7Aw0m7h26hG4cXZPoEme9trp8a7edNzEdJ2Wb5hv1Qicm00k\n9WEkLZPBlrSAuEFpk3Grhjghvbgxirhd75YXJxW4P/TkVa8iLx8AHv94So9iGR/PLoHrAnef52cb\ngTww5P/LgE2/RsC+DnsClI8ycZXbjwHxefh5bdGO+frbdY7vekL7Ma9enk3JoC1pnyK594reXGMR\nJqcI0NJGExc/8Niu3pnrnoO7sQIA8KfoaThxzT8CzSZa9SE8/qbTu/fg8NX/jrPbLwYAnJseio/e\n8fZuGe+57QOYjsgLv3D943FN7UAginBPtAI3b96122cb4kWYqZHefKuO+YSeq7g+QsDaud64XUdc\nH+l401H3WFwf6QI4e9pxfaQD4BEBOpfXEqDdqiGuj5pzOmV0y+N6W96P23mlAvcHg4yMAC8Q3NvF\nF5u17fff3yyjIIEH8IOnDe4agPs4ch83X4TDdxkO+9cHhnkGIq+8kH391CHLzZvPUOS6+zUCtVpv\nW/g817kd/c1Y1G3T1a3HIqlR7Oja2gG4ZO4goNHAFBbhiKkTunqvn/wmrpnbFwDwg/QN+PCmo4BG\nAw+kS/GiDad29b616bW4MnkcAOAGPAYXTx0INJuYq43i2vn9kDboGjfHw5iqEY05XxvpAnPaaKKF\nhgHZBEjqFPOyPeg4bQjQriFuWyDLYNyuI64NdwC3ntVrGzBO2jWqiwGcQbs23D2ne6zWq5fUh5G0\ns0agrFTg/mCUZcvM9utfD7zudbS9997AJz5hjo2N0R/QCwDSCLjAOBRE8oDNpVfUk7S3QznyUN48\nz8P2tbMfIxBi1IoYPzu462n7usbO2FinXOrZ2hjOwYu6esfjo9iYUkrv9/AmnIS3AgBun90ZT8Mf\nu3pvxUm4aupRAIBz8CKcjlcDzSYmWtvhh51tNJtY3d4Ta9v0WcXNtSWYSgh852sjiNEUQFpDXCO6\nJa4PdykLNiDSM47rpJfUhpC0kNFL6sPG064Nm/IkaKcNA7LtWpcS6QHjtJGtt+OFExhbRqDOdEud\nymMPPwPaUq+GpB2ZetNG14Al7QrcKwHIw3+9+OLhz34GPI48IOyzD7AveU0ZMAL8YMzeWp6e/VtU\nr4iHH2osiuzbEoBf9FxtBGR72gK000YTcxjqlnElDuweu7D2XNyQEODeNrcrPoxPdfWeiCswEVOu\n/Ofxb/j2zD8BAK7bvCc+hM909b6PN+KuOZq7sSraB6tq+wAANqeLMInxbv/MY8iAbDSEBB2vOBoy\noOALXT8AACAASURBVN1oIEHDeNC1ISSod8E4o8feL4CkNpyhNgBkPGMJ+JLa4H1oNhEnxgjE0XAX\nmNv1JtoQHnkSIY6GgSjKgnGjgSStI+Z62/UuoJMREKAtPfy2RbeI8mLh4SftGrWpYwQAoNVgw1SB\neyWa7LqrAeZXvxr4f/+PtlesoO9Ysvg8/FAjYP8OMhBog3GopxvCzQ8a8MucG0WYqo2jVSeg2pgu\nwTV4XLeME3AU2nXS/1x0FC7Y9EQAwO83PBavwQ+7RuB5+C3WzRNN8b94Hc6bprkWq2d3wfl4bre8\nu7A7NrVolutMtAizNbr3cTSc9aDRNB5vNGS2ax692jBi0JorSX24C/S6XrMLkG3U6RpVvUZXj/d1\nwbMu9BggbT2mQIAMyMp9BO5GL6kNZT3ytIGkNtwF5izdUs/oxazHnnuz4+Gn9Qx9kx0x1BzX2Mlo\nKyEVuP8tSrNJ61awnHYafeQYoMlY7O3b4C49SXvb/vUBWj80TyjnXkTPVX/JgOokFndB+8rakzCZ\nUvD7+qk9cWbrFQCAVq2JQ3Fu97x34Fs4Zy1lR517/wH4OI7v8tsfwacwC/Lo/hwdgFVztMjdxtZi\nbMKSbtvnMYQ46niyNQHGDNrMM6Op60kPWgPZUHBnbzrXCBhwB4RHrulJ0I5oHXnSGxF6lhGIhnrB\nXYBxl76pUXlJKzIjAdZrNNCqDyFFLQO4XbqlQ6N0PXJJ37SzdEuCZpducRkBM8LgtlfgXkk/sssu\nZpbri19M61gDtK7FN79p9MbH3TNt7V8NIMvSLXlGYFDcPAc7O8da9aEuPZA2mrgUT+6C9sm1t2HN\nLFEbv1h7ML60kYzlZDyC/XBjt9wP4dO46N5HAgAu3bAvfhL/HQBgDsP4LQ6llx3AZLQEk23ypucj\nA5Bd6oDBOBoW224wTqIhJBGVTRSIAKCMRz6cLZuNAOuByohrQ0RZIAu4vcZi2Ohx+1xGhUFbXI+5\nDkuvma2jp7zIrZegkQHyBAEevgDtXr06GZLO6IT10kYTCfetNAKiTV0DluH66z2jia5eH5x7I1+l\nkr9ZqdVMsBYATj3VLFq1//7Z6dXy4+OhHrkN/KF0S1FjYWWoTLSXA2mKHQDMtIfxQ7wRb+ocOyL6\nGo7evBQrAHzztsNw30QDxwO464HFeAV+gnuaVP538Fbste5O7Axg1fQuuGGe+mW6PYIN2A5oTgLo\nALUENAk8IA++ASCJhCcZDSGRgGf9ZoC0A2hpo4kYQxYY60YgC3ZNxLVWZ7tXTxqBVuf+qkYlMxIw\n7WO6pdYDxr2euwH31IChqKvrVXf0qH0StBuqR+40FjZoi/JUPa63Jbh5YcCYVktqQwa0ldFEU4J7\nnXQSzpaJmpk2JWnluVeyEDI+bjz8pz3NfPBkdJQ+WsCy445mMlYZj5zrcKVJsl4UYX1th+7/10/t\niVundwEArJ1ahGPnP9xt0jNxEe7eSPTIibe8ECeuOgwAsHr9UhyPj3Vn6/4ufTbu2kQZIhtbi7Gp\nTefMYwjzHLwEskAq6JAuyHTaK49lgJnP1Y5FQ4jZ67Z+ifu2QLteR6tmlQdFTwPjaDhjSGyO3Bic\nJpKaaYscCWSMRdTslidHDqoeaBSUAbSOXs8oQXq/DJ4wdWk0T9bDt0YTqlGReqY8G3BjcR1xNIQW\nKLtF0ltkVBqI095RQlpvIOG+jSIyYJ10ykzbG42Kc69kG5DnPc9sn3yyWQP/oIMomAuQ9zw6qi+N\nbIH2DY3Hd/8/9/4n4Np1uwIA/rp2GQ5fe2y3qiell+G+jaMAgJNWPw9nrDoIALB63RL8eP4lXb27\nsDs2zJLeHIYxh85LJKgHAIijZhdY6WVXqI2uXi8YU3lDXbqFwI63s3p8brcMsU/z2E158hwCiFxj\nITz8FhpWXc0ePZX28NFBDr3ur8/D10C7Zvo5cQVUpUcuPGhvvXLEIPvM0pOjmEz7eCRgjcC6+oC4\nJ0NdykYaNfbwM33RpWWstqflyZUK3CvZotLe+5GY+0eTnnn+t1d1wf2U9htwVfsAAMAfb9kR/3Tj\nR7t6hyQXYsMU6Z1x11Nx4W17AQDu2rAIV80+uqs3iXFMJ/xCmBe26xl2JIZ5USSgS/ClYx3A7JY3\nlD2nA+4JhEdn6QEQnrQAApg2JTUbALN6pg2KHpeBBpKokSkvyYCdAnw9RsCut5nl5iPd4PRktxQF\n9yJ6Gmh326d45HJ05DICCn1DIwtzvTwCywBuvd5TXldf0FFUrhmdZGMSBtS7bWJ+X1xP2miiVdEy\nlSy0xDF9/wGgWepnnGGOff7zwFVX0fYvfwm89a3m2AveuCtmZ2n7V/WX4Ko6ZYjcv9sBuHe3J3b1\n5qIRzDc6KXoysGd52gkaKmhLYCI96alZ3rksL6KXnI4N9dbLdIsog4DetCH7O5Rtn+W5c12+9vWU\np3n4ll4SNTKUTtf7tcA9sQwJA1+r1qQMEdkmzeAweGYoKpsikkZAXseQW89uUwdUneWx0VM88gQN\nY1RlX0QKaAuj16VbuoDb2+/SGHjjI5qRVPq9e/86MZM6EpSVCtwrycjMDH20HaClRL73PXPs6KOB\na6+l7TPOAN7/ftqenc3GXX//e+Cvf6XtTZuMEWi3qUyOw/JaVAAQDy9GvMjMvI3ro4iHKEiZrHg4\nkj33pv2N0e7LAXS8YQWoyVOXHnkjq5fq4J7x6iV9wXRLrZ75v+ccBZQ1Q6L9JtJYePR6PG0oHn7G\nMDW6ANnL4QsjxQDZkH0lyrOvrTOTMgvovaOYjKctR0UOI0Ccu+mLuOP9yjZ39cB0hyOg2gPa2kig\naWgUBXjtvkgy19gL5IkG5A76rVuH1XautxlV4F5JjmzebEB7fp5ocZZ3vcuA8X//N33MHQDWrweO\nPNLoXX45cPvtprypKdqOYwJt/v5FBrStbfuXt3m5cIDK4e8rAEC80wrEO+5G24/YF/Fue3bbFEdD\niBfTBw2yHjQFtLp6EoBTG7QFuKdZADZefIduaYnyHCMGaotjNOEC7UjX63r2tofqKk9tkwIo4pg0\nJEy32Ho0oskahiRqol0zIw/7HmQMCXvQigHr9pXTCPQaLtbrMexNq0868YgeAyv1Utm3vaOJ3r7r\n9cgTrZ0ynuExAl2jWqtl642GKnD/W5bNm+kPIK/7lFPMsTe9CbjjDtr+1reAE06g7XvuAT7yEaN3\n1VXA3XfT9vQ0/QFZwAV6wdgF2vYxuV/+uoyAPN5zrF3rggHQAeNmh7555KMR7/tY0htf3p0QAvSC\nsaFDGpmgVZzx1s12V5/bljECDQMQadYYJKgLUBRlW3q2UYmh0zxEMRjwdJYXalQ0w2TFLnxlSOPE\nv5peNwdcGjOboooM6Epv2TawSdfTtsrr8cjlKMaAcdeAKaBtGybpaUsDa5wIAcbqPbCAXI4gtfYh\nW28F7g9BmZkx39CdmgJ+8hNz7F/+BZiYoO0vfQn4yldo+9ZbgeOOM3r/93/AfffR9vQ0lQlkwdL+\n37WtHdNAWDsW4rmXMQI8Wugeay5C3CQqJ95uJyRjS03boyHEO9EysvHQIsRDHYOQNrLGQr6kqQKy\n3EbUDehoehpQp81eTzxjIOSow+eRN7r77PJMufroRJ5r9BRQyjEWRfepdJVtGNh7dRmmEANWq/Xq\nZwyOCG73pJv29jtRRL1xjV6jJgKqyigmm4FlA7lOUcXREJoo/w3VCtwXWCR9MTUFXHSROfaOdxgv\n/Pjjgf/6L9q+/vrsYo+//z2wdi1th3raLgrEd4xBOo9u8XnaPnC36yprBPL00hRotWvmZd97P8R7\nU8ZNvHxnpKh16ZYkrSMRwN89p5Nv3C1fAKD0prv6nQWfMmVY2/I3kSCbZmMCmd9UjBK0EYPSpiSt\nk7faBuK2PbLoHZ3ESeQAUmMQzD6L2kgVcE/NAliZc1PrXDTRboNm5CpA6jUQHiOQ6Qtt1FPQWMkR\nSOI7N1L6Jxry15Fm29dA58EsIRW4bwGZmQFuvNH8f9RRBhSOOspQJ5ddBnzYzLPBWWcB69bR9uSk\nAfpBeNouqiSvPD6Xf31grB2zz7F/B2EEZPsKjwQ4d5r/Fx5+suOuiLfbifY3RrN6DtDuAmRXr6mD\nsQ2oArQJjOrd/fRbN+WFjBhk+7iuRABkqhkS03avnjAovaCtU1SyP7rn2m1PG4L2kkZPH8VIg5ig\ngRQ1tNvKKMZlcIJHHb4Rg55e6ytPM5zaiC5Bo/LcF0rimIKMLF/7mvFqP/QhSvsDgAsuAN73PqN3\n4okmY2TDBhPYzANjDYzK6rmoEi4jBGRDwbiM3iCMQBFjoR1rt4F2OzL7lu6AZAnNtOVMnoznvnxn\n2pbeGAMPGxTUzT5tJCAMhAR0c069U45C32gjgbZmBAxo9xihjCEx1JPdPqKh6plzEzTEdcoRQ6ie\n1Za0IfrXc43Sw7cNTgxFXxoS0089nn4G8OUIwx6xyNiJGAHZejlGpcvvawasYxgqcB+wnHee2T7h\nBODKK2n7F78A3v1u2k5TyiRhXnz1agpUApSNUsbT9nmooWBcRE/SLS6vW/76RgLyt9+Aaj9GQNZB\nYJ2v56tfa3uaAkm7jniUlimIH7YX4h1oBi2nb2Y8/K4HWe964D00TgbE6r0AJcA+EQBt6BbLWEAB\nRQW07Xr5GrvXrYwsMudqHrlPL0M9KX1glycpL0kl2e0Tow4JkGkK0FwDywgkkdJnzUzf9l5Pr6ed\nGXWIazP3IDtC6+bZZ/pMj4nEaQNNzKOs/E2BO2eEAESN3HUXbf/kJ8C//ittT08DL3qR0fvDH4Bb\nbqHtqSnDb3dT4kp40C4wLqIXYgTyQFteQxnPeEsHVPv13H3tDPXmtX1dfp7LbTe6S7NKIAU6QNBJ\n3YyXbG8ye1gvibr/d/lwr0fe6/0nifTw+ZgD3IXRcdYVqpdEYXouWkaLO3RB21yvZsB6jUXv6CRu\n10Xqaq+e9JwB6sPYAl7pkSfovcZE9eZ7RzEZw90zOpB9kfXwm6l4oQvKQMA9iqLDoii6MYqiv0ZR\ndNQgyiwrf/wjujMgzz0X+Pa3aXvDBlrmhOV73wOuuYa277/fBCgZDEKDiIUzOhxeaF5dIUYgtDwN\nREM84zyPPK/tWnlFOfwiRsBXv88IhHrzeX3WpXmYlx7fHglz9/yCcxlRE/H48s6xIQP0IjALEGAZ\n/t0NxuTh17vbPX0hjnFdiQQ+i/smQPPVW3frtetGTzNSkpqyOXfZJo2+8YxOMnSQZfySJAukaUrU\nV2+MQwHtzMhK0GA9tIyDqlGuuwVqQ8bQpQ00tiYtE0VRDcDXAbwQwP4AXhdF0X79lmsLc9YAcPXV\nwCWX0Pa6dcAhh5hjH/gABSoByjK54granp01k26AcODLDFG3YHBQoxA04Mmry2UEfCAb2vayeoPo\nM7lPM5i+ayy7r0x5ZQwDAFojvMPrJzuvQLsTHOzqZzjbTrocA4sstwuedTESMHpytMG/BqA85bV7\nRwxZPW3EUO/q94Jxvbe8tN5bnsMjt6+jS6ckeUag3lueYlTk6CSxjZXDMIXvs8sTAW/h6SdooJlu\nXVrmyQBuTtP09jRNYwCnAXh5mYImJky2yAMPAMccY449//nADTfQ9jnnAD/6kdFbtcroSb47FPhc\nLyf/LkTOdh6w8HaaZqfw23psJLTyi3rk2q8rM8W+1qJGIJRLDwXUvBhCKGiHjA5coG2f66qXjZUB\n+lpPGfyh5BhNxDvQssZdT1/Wz3nS0sNvN3r0klQxAope1wtOa+pIoMczzowYxChC0etmt0gqqQvG\nisFRuflQPeUaixqVtE7vnxbPSBuir3LA3TZIcLQ9rW91cF8B4E7x/12dfarMzwMXX2z+f9ObDKB/\n7WvAN77RKeQu4H/+x+hNTtIfUB60pV7mwXS8iL7y+gGqUCOQB1qann0NRUF7EIZp0IYuD7T7AeMy\nAdVQw9BPvdoxMuwRfdgBQCKCt1098Zm2uDOJq7tP6gmvOklr3e2eaxQGQoKhs7y0Lrx9j54Gsm0F\nZDN6nesWYJwEl6e0XY5KrPKSQL0uULdrwus29XO8IGPo0OjRi9Du1NVAHQmBe6v+4MqWueqq7Hol\nF15o+O7ZWZN9YoO2C4xtntkH2mnq9mqLApUGpEVf2FCqpAilwr893mCOgRgEaPdbHsc6thSX7tu3\nJbl5qeeLwYT0j6ZnP/MA6GPP9Q6/r3n4I5TpI2dgJgpoJ+1A0JaAGqjXPdblwzXQFgZH0Ewh5al0\nEIOxoEDkyEI1PtYIJEJbeOl11NDqXDdlt8TtejdVsxElnbbXMYzZbv+MYjozShirzXbrGsM06SVR\nX8sPNEqfaeRuAHuI/3fv7OuRY489FnfcQV75hReuxMqVK9WHm7dd/7u28/T4d3h4cIDWr1cb6uEz\n8GmAwR8ukvv4w0V5AGlz/b5Rhw2AITSKzzDlURdFvXQXVeSq39dmX70h1+F7Bor0Y5lnqt3OGkn1\nXnQDtU3EQx0+XxiBbpt4dNCuo9WOuts99QqQ5Rm/weDO2y3JpRs6qNumtNeD1uIKnLmUoYO6Xn+t\n4yXXhJdMRqAZxZn2DWMWSTKCJO2AcUdvrD6HmfZwV28M04jj8Q5oTyJuE4ffiBI0ay3EccPopYsQ\nJ1Fn29S1tLkZcTyGuF1DHefhFxdchuU3jeDW9GaUlUGA++UA9omiaE8A9wJ4LYDXaYrHHnssLriA\ngqErV9K+bQ3c88CmqBHIA4LQ8nz1M7iXAbnQfT0rNRZoeyh4Fb3ufvaxDLreQYze+old5OnZI9mM\nfqvWW4ageXgifKzRQRzwbUnQVsrj4/IY75OplW1DFRkw7my3a13w7FI17RqGonkkyRCSVHjJMTBa\nm82ANnnJYwZw2wa043bT6EWziOMRxK0a6bVIb7ieYKo1Sk5RWscoZhHH44hbNYxipltes9ZCM2rR\nNaZc11KqCzOI20PdVM3hGnn4SVrHUhyE5z7rBXjkgYsxe+bPcXMi+OkC0jctk6ZpC8B7APwGwPUA\nTkvT9C8ufY1GcT3ooaDNXm2env1bBrQHwaWHtM/10g+CTsjb5zvX55H7QCmPbumHS8+jnvLaV7Ze\n7RnMu1daudpxW28QjoV9zJ6r4SsjSSIknZz87jdSJR3U4fe731lNgKQ+nNGL0O564knbUBZJu4Yh\nzHXrHanNZcB4tD7fMQJZkB116fWA9rwC7r3ljTWkXo2MQkfPBu1GrdU1TGPRTLftY5hBktYMuLMe\nt6lFx0ajWSSdrJ9GrW3K63r4pq6yMgjPHWma/grAviG6Gmjb3gSQfZjt/13A0my69bQXzQcoeS9i\nnrEoawTKAEvno/SZfSG0zJYAQ1eaop3h4yojr66ilImrH+VzVqR/8mimMiAbqqcZ2LLl+d4Nuww2\nzFFk7mO9bvSaTaHfqpuRcTSM0VGjNzaaIh7tBHnRwOhwJ4iYNjAazSCOhw1o8yiiVeuA9gjiNoPx\nmAHj1lC33rHaHOJ4sQHt9mhHbw6bWsNCb3MXtEmPyhutx0jmjPEZZb0kIu+8vdgCbQpGj9Vmuu/8\nWDTT9c4luCftGkYxi6QVIUl4NOEzAgzunSBhCVnwgKoG7q4XUXr4rmPag+7atvXyAEDT873YZYyA\nltqo6fnKHYRXm3eN/QCudsy+7tDrcQFpkXNdz0XodWvPnKYf2hf96mU8aMd1RJEpo1Yz1yGBOUmI\nsuTyRkfNdr1udOMYGBvr3ea6tGNxDIwtrme4/rGlQ3QsGsLYTouN3pIGkiYFg5O0gbHhVsfDJ+BL\nOqA4Vp9D0o0D1DDW6Hj/XY+8o9eYFzQPGYskAZJWhFHMdD3o4XqCFFFnNmzNePhxB7Q7MYFmrY1m\njdpExkLo1YRePas3Gkm92Sx9U+sYulbWw+/Hc9/mwF3bztOzfwf94ri25W+Z8soAiw9Qi4J2kbYU\nHVnk9U8IQEqqI6RNZfoxr0023eLqR61+7b74+iLkGWWQ5n2aBy1BmuuXQF0GmJvNfHDn4PzIiF9P\nrWtx3bRv+zEzezQawthu25Feu46xPXc05ywRs0fbdYw1YzqWMGgTNz/aiE1AtV0z9E0MjEZzBkjr\nLTS6IFsj48Ftqs10jUqz3kKzzqmLNRppcHk1Ux7RLWlXb6wm6JvarKi3jWa91a13FDPdYHDjweq5\n583K9E3IcYFnv3qD9sbK0i2hQBUK0GU807x6tfvnAqqylIW8/yFtt/u9H8DPA+28vssrT3rTDNJ8\nfGjI6GWojVAv2dKTZcjtYWIruqMol14IuIfqsbHUrkvuV9u+/aJuP4/uvoOZ3Vkbwug+K8yo4+G7\ndCcxjS6qZWiesYYBd/K0Jchmwb1bb8YIZPVGa8JY1Gw9A9pZD7+3PBqdEH3DweVm/UEE7j6uz/Ui\nakP4UHAI0fPxqHYZoV6WzT0PggqxqZ+i5ZX1vgfV3yHA5yovxMN30TxaXaF9xvp5fWEDNZ8zNKQD\nFXvQ3BYXGEuvW9MLBXeXXqNBf11wUsCd6RuOadkjAbkdotduEzVkU0B2vXZfuLZJL8Lo4oa5xr12\nRhwTQA7vsweAiJ6P2jDGHrOX0dtzx+7CZs2hCM16x9NuZT38sfqs8MhTywjMmnsqQbtmGQvp4dfn\nuyOLXg/fGIEHFbiHAqQ9NJZ6LgPhm5qfZ0jKglJeXf0A75Yqr2gdvn6UdYUaSZ9e3n0pSoPZ7ZP6\n2iSvPNrD9qZZzwWePqCyz5Ft0cqQAMnX4yvP1z4bjH1lxHGYESiqJ+vVroPvT9lrlPfN1BVhbMdF\nRu+ARxkw3nsPNBc1qW9qwxg9YF/Tvt23N0agmaLReS6SVg1jHQ4/jk0cII6BZsOAO9FBhuYZrc/r\nI4HkQQ7uvhc2FMg1PRdnauvxQ5Onp7XXB3J5QOUrI+QaXSBW1DDkGYt+ZlEOQq+oMZXnuPqCg4ic\n/lerGUqJPchQwLBB1gcsofRIKFBpxkLzkn10iwZ8fJ4EfjmyGCQtw3rSCMi+5e1Wi4wrj37sNoXe\nA7tNrj5rNCI0mwTMcRJhbI8dTL3PeGJ3wbLmHruhuf04tSNqYuzg/U29uy3rLkTWbKRoNjrtTaIM\nfTPWmKdZxDHQrMuRQNQdCdCxBxm4FxlCb4l9oWBoT9zJ80zLAn4RmsC3r6iH72qLph8y2nLV4Ssv\n715o9bF3TS+kOYdT8FgvD9A0YOGXnZ/REF7YRZ3wdqhHLttehm7RPGMXPZLXF6F9Nkg92X+2UdH6\ndniYnoUk6b1XeVSR1mcug9PTpuEamkMEzHFSw+h+e5ryDnt2dwZtY9cd0dh9FzqGIYw94yBT3q5L\nTSpkMzV9kRCH3zU49U7ObgnZKpw7oAPGoEBuUPu0SR4hqYsaiOXt00DRV56PCunXgGj62nFNzwXa\neeXFsfGu+UXj8jj7gvVcgGYDru3lhQCQBFybOgjxyF00gtamssbC1RchnrHdFz56ZCGNQGi9Whxg\nZkY3ZvIafaDta7v2/AT32XAdzRFKd4yTCGNPfLQp+xUvNLnsO22H5r6PoLZHDYw+/xmmjHGaIVxG\nFhTc7Zc6FGwGva8fo6Id2xLB0y113Tbd4vPIQ+gR7Z66eGumC+wXgfvMR0WEDL9tKkIb6ktQyHvZ\nGRQ0wJDUAdB7XfY5dttdQM96roCq6xpDgMqlt9ABVU2vjPHh+8Pg7gJwrmtuzp+q6TM4oSO/0HuQ\noYOG6miMNjvHIowdcqAp459fi7KyoOA+NRUOSqEefhm9fgDSLjcPIOU5Id/09PVFESPkapMr6yiE\nHuEhMB/nl4lB2/ZyWE97qF2AluetaiArh/AuoJJ6EghcoO162bX2hYwEQoKDIfSNz0MtCu62ofMZ\n2IUMqA4P0yhudjZ7T13gOT3tv17p4Yfo+QwO95mrTa44RVlqrEmZnqVkq4F7KPCGArkLqEIBX/O+\nXRSDXa7dZq0u17VpbWfqJ4S+cV23ts/Vdulp88MrXzT7oeNzQ/OUXS+Ji2LwURtchm1INI98S4G2\n7U2HAqk9y7NMQFUzYD7PWKtLa7u8fg3cFjKgysdCQbuMntZnbHB42XF+9kP6wn5+7NGJfQ9GRsjR\nm5vbcuA+kLVlQsX23F3epb1PWztFA0of8OXt4zLyDI12TDuu6bmuN88IuI6xMbJXhYxjelBtesQG\n7TxKQAM02eYQ0PZ5l6Oj9G1brtfl4fOLwKMelxHwBd9sPR6FaOVt2KC/xIsW6dfr09NASzs2PGy8\ne406mJzMB09fQLUsaEsjINNCx8aAjRvdQJrXPh9Aap62bXxcevbzI++PRt9oAVU2Ava9chnz0FEM\n98X8PN3jRoO+IjczQx8suuIKwsg//pH2T0wA112H0rLg4B4CnnJfHIeBe7/7tGMhei6jwr+h3rTr\nHK08Bm2eDMIfIYljN2gnCbB4sfsllnUsXeoHSNYbH9e9En7xWW/JEnd58lqXL3eDonyZtt9evw4f\nxcCgzeXJSUe+4XcoN+8bzruG6RIgJX/sM1I+YJH956Ooli3rbbvUm501z4Eszwb3iQlzf13tk20f\nH8/SLa7rsMFY05P3x+WRy2u0+3Z62t23dnn8Xs3PA/feC2zeTFj2pz/Rc3TppQTCMzN07vw88MUv\nAn/9K72X115Lx97/fvoWdJIAp55Kvx//OJ0zNgbcfTd58rfdRvsWLaJ7VVa2queugZe2Wt8gwD2P\nqtHK8IG2T69MeRK047h3OrpNj0jAkOuv8MNar5tt9no14Lc9ctcxuV/Ts4FvYqJYeT4PX4KnDQRj\nY/TiaEBle6tr15o+c4F70YCqC0jn591esrwHDJAu79JFt8zOmlmerjaVoVFGR80owQZSOfO2TEDV\n5ZG7+kmjW4oGVOfnDabcdBNt33orcN99wJ130pfgNmwATjsNuOceuh9HH039+6lPkfGdmwN++Ut6\nhz7xCfqNIgJ5ALjxRiojisiBAoB996X6R0aAQw4Bjj8e+PzngR/8gJyTV7wCeNe7gM9+lr4HIUvK\nnQAAIABJREFUvWIF8LjHAR/4APDv/071PPKRwDOfCfz4xyglCwrumzeHg7vc5/sYRR5ob+mlC7R6\nXcbFBu041rNH+OGcnjZ6EtDtB1wD90ajWPqf7CvXMNMG7bKcu7wOGevwcd++ANbkpE57+IyFzJXf\nEtw8OzKSW3UZn6KA5qIitDbNz5t7atModp/JvrX7udHIeu42f2wDbl7b82aoNpvUXgbR1avp/xtu\nANasoe8wz8wQqH7hC/RcnnYa7ZuYIJBcvx744Aep3QDwyU9SWRddRH2yaRPwwAO0PT4O7Lor1fuC\nF5D3fPjhwJ//TPV98IPAP/wDgfFll1HZL30p8J73AO99L3nic3PA614HnHkmcPDB9Aw0m8Duu1Mb\ndt/dTJ7TRhOu+1hWtprn7qIzbHBPEv0zcnmAa5fh2rb1NWOhGQ27LRq/7fI0ZTmuQBd7UAx8Pu6S\ndeQLyeDuogTsNs3O0n4fAPG2a+En+3q5L1wBRh8HPTpKL18eyPoMmC+rRnru0rvMA202JKEBVb5G\nFy9s1xUSHLS91TwvWQNSu+2a0bLrtWkZ3razW7iu+Xn6feAB+rRmHBOw/eQnBJhXXEFgvXYt8JWv\nALfcAlx/PbVlchJ45SvpGbrnHkOrXXwxPRdr11IbhobI6126FHja06iOyy4DPvxh4GMfA97xDmrH\n734HvO1t5EG/971Eo1x+OXDEEQT2z3wm0SI33gg86lFU9sMfDvyl89khSWna3LzsP210oo0stPst\nDfGDGtw18GTg4OP8q9EyGuD6QDvECMhjXJ6dSSJvkB30kvW6AG101Azn7BdNnjM0ZIbBEgg0AJIz\naVmv0SjmheYF7KRB4GGpNuqQ7eM+K5OuNzpK3lFem1y0h8Z9Sz0toOoCPlkvtyk0S8d+iTWv1kVZ\n+PqsiIfviwnIDBEO5Mrrj2Pg9tvpmW21CCTn5oDzziOvFiCABoD3vY9AvN0GTj+dwOrDH6YyFi8m\nOiJJCDCZJtxpJ2rDS19K5zzpSURHfP/7wDHHAN/8JvD859O7cPLJ5H0fcQTwvOcB220HnHEG8PSn\nA1deCTz+8QTQ9Tqw117mW7I+kA29VzwbVhow7R5oxle7p1oZ2r1q9IHQ21RANY6Lgfsg1z+p1bLB\nE/vms57L0w7lo30ZIr4goiuDRXruWps2b9bL8w2/XZ5cqJfMbWi3CRAkYIyNUWYA1yU9z9C+COWP\n5THOYOHybM6duW9bz1deq0UA5gNPF2hr/e7yBvPA3RccnJqicjduJA/43nsJlKemiMv9y1/IEz7m\nGOqTj3+c+mJ+no63WsQ9z82ZNV5aLboe5pgPOgi4+mrg1a8Gfv978qSf/GSiMn74Q+DTnwb23pu8\n4499DPinfyIK5YADgB13BM4/H3j2s4Gf/YyCszvsQO+gBLu5uTCDKEcx8nkMpUBc94DbMT3dS41J\nWkvq5RkBF7jbHn5Z2WY5d2kEtFQ/1vGVlyRZblW+WCFpfRq4u4BlbIxekrzyXB6+/UDKobQGfNK7\nkuDO4MSeO6f1SSCN42wmhJYOxvWOj/uBlF+mpUtNlgWPOmZne6/DZyDkSEhmxOQFVLU+42PsUNh6\ncqmDMgFV6aFpHrlGj4SAtiugOj9PowYO8F1wAR3/4x+Jorj9duA73yG649OfJr2pKeDlL6fn46tf\npWeu3TYgzRkZzSZRGldfTRzyRRdRve96F/HMJ5wA/OY3dE+f8xyiM975TuDLX6bnY+VK4H//F3jY\nwyhjZ9ky6l+mGPL6IsQgMnYUDajGMRmMkNGOBFz7WdLa5NMLBXfZ9slJXa+sbLPZMnKfFlDl/20A\nd62fnQfa2jGbYnABAeutW2f0JChKb1pOH0+SbBqZLwNju+2y5TEfPciAap5eCB/tehF8dXHmh/xA\nRWhANY8eeeABvbx+AqqyrjyOnIODDATsRU9NEYWwbh2B5fQ0tfULXwDuv5+yKqanyTi/7W0E7l/6\nEvXR3Bzwq19R+exQDA8DT30qcM01lImxYQMFBd/0JuBDHwI+8xmiVFavBj7yEeAlLwFe9jIC9Pvv\nBx77WDKu++xDBoPfJQlA9boOpC6QTZLeYKtNP2jvkqtvXUagn8lOMl6wbJl7lCXLWLzY/XxLvZER\n2t64sdeYudIzHzLgLkGaMzVcAVVflknZTA1ZP4On5iXLnO3ly/2cLuuFZla4aB7tYZVtZ2pDeu6y\nTUUDqj6A5PsmwU3zlOxjvpcu7yXWDBPzwsPDJrBpj2J82TJSzxVQHRkxdEseYDQaRHPEMVEZV1xB\nAH3eecDNNxNoAuRxv/3t5FVfd50xYqeeSl73qlXkDDQaRGksW0accrtNAcQTTqAg4DveYbz2t7+d\n+OdXvYqCkJdeCrzoReS9r1hB7bv/fj1LR3suQqiDer0YxSCX63Xp8WgvdFQ0KHDXRhM77JAP2jMz\nhBUhAdXx8V7vXGvT8uUPMXDnzpAeVL1OLxZPyOEbHkXZTI0kIespwdPeBvRj9s2XetqDIT1yjcuz\nXxhfeS7wlO1jT9NlBDTvxebcGdyL1OsyPtymJDGjDm3KdNnMD5eeBkChwUbNQGgBVVkG50NPTNBz\neP75lAu9YQOB86ZNwNe+RkB88cVU/vQ0URdJQrrLl1PfMB+9886UcbF+PQUBTz8deOELyYCcdBJR\nJe97H/DiF9O5P/wh5URfdRWw335kOADjUIT0GW9rKZMMaENDxlvVOGJXXUXA3RdslAFLbTGvvOcn\npC/swLALZLV+0gKqsu18zrp17mssElDdeWe38XnQBVR9E3JkYJNfXAZ3aQSkR6ABachCTeyR2x6g\n7SVLI+BLXZQ52yFcutam++/PBmY0vlcLqNojkhDQ3rhR56PlCx4C2i5P207Xs9sUkuubR4FofTs7\na8B7YoKM03XXEZCcdhp5unNzBNRzc8CRRxJ4Nxo0qaXdBn77WxodNBrkgQFEZZx9NmViPOtZxG1/\n/vPAKadQpsZuu1Hw8fDDgaOOopS6xz6WjMF22xG/7RpxhgZUeajfz2JZWn+G6BUBdxd/rNVrP3Na\n8NJ+HvM4fDZgPHIJzVph2szW04xgFVAVogVFbFDgQJcEnVrN7YVKPbmtfcVFPhijo/Tiaznbtp4E\ndx8HbQfBZHkyrU8DOxdA+gwYP3TcNnnMToXUQJYzRHwB1ZDhrQtYihgB3uZhutb2NWsItG+/nfKR\n77yTZg6yV33HHfSMveENdF+POsos1jUzYyYTLVlCIPmsZxHof+UrBPrj48A//iNRHocfDvz0p+RV\nHXAAbT/72eTRj40RWA8N+b1L2zBJL9EHxlpANYTrl307POymxrQyZExAApUERebcbSDV7mMZPd9z\nMTdH91QaR19AVbvGTZv0uux+0oyK1qZBB1RdfVFWtgotEwLu8litViwNj2fihQBkq5XN2XZ5yYBe\nBk8z97VpbCw79V0zTK6X3dV2+dD9//auPUivosr/eh5OXmYwk0mA8NAkPGKK8IpAMIHwUkSIPKSK\nsiq4oAILi5ZlyZYVLMGVxy6Kq4KgJY8ySkm5UApERXaL1BoNorySsEjAtSJgwSoBTDKTTDLp/ePM\nqe6v5/TjPub7hi/9q5r6eu49t7tv39u/Pn3O6b6uQ/X1143GKeXnltvbS8SXameuojXaMwteObhh\nAx175hkyb7z4IpktXn2VFp288QaR/rp1lOeDDzaa9QDSnPv6qG7nnUdOxDvvpNWCW7cCV1xBkSNn\nnUU2a6UoPK+7G9h7b7OKOUaePpL12fp9ZJc6ILr5SWYEm9DZ3MJOv5CJYXDQ2HtjJNvbS8/rzTfD\nDlX2HcRMIBylE3NeskNaaj9bjp2/ZUk2NIspMzD5ynX9Ulu3xutXFk3d8pc1d0nTlpx0vkFAchCV\nkfM54iRHqftRhhRzi6SdxwamUJ0kE02ItEMO2pADy0fGrsYitRnbnIeHzRLxxx8HVq0i0r75ZjKH\nXHMN2Zw3bqQIkYEBOr9lC5HIwoWkXV9wAcVFL10KfPvbtOhlxQoKyZszh7RsADjySLJv9/WZjp7S\nie0ZjU8ztqN53OddJLwutW1TBthQ3WPapVv3ohp+s+QkO7hSo2cQ0rPq6EgjWfddSCXtWJuVWaH6\ntneouuYC98V144937TJmGXs5umtikNKunG93QrdBe3vNognWhmwnos8e3dVltvNM/dpLquY+ceLo\nwZHlqjpUpRfcrtPQELX99u10b6tWkZlpzRqKid60Cbj6anIAPvWUyfu++yjPwUE6NmkSLSm//XYi\nbLZxX3IJ2aovuogW1axdS/btu+8ms0hnJz2HVCKVoiwkkpDyK2ozdUlm9+60gdM1E8YGgZTBwjVT\nuPULmYrqjpZJdahKMwu+D/aD2HJAuCxXLkSyto3cfY59ff66223G4ZN220oOVZ5Jxdqsr8/MaHmB\nXFk0ldy3bKGbdJ1+rkPVJtmdO41D1Q7/sxvKjgFPJdJYfDSbkKSO4A4+vs7udmKbPHgfiQkT5NV3\ndThUY6GVQ0NkwuKojJ/8hFYx/v3vpG2/+iptU7p5M/2/117UDk88QSQ2dSpN6195hXayu+UWip3u\n66NFLZ/9LC10WbaM7PsvvwwsXgx85zt0352dcYeqT2vic7YpwkdUkv3YR+4SGafYTLu7Ry/akgYL\n3wxRsrvys5P2bkkl2TfeSBuYYmaPusi9u7vR2Wr3KzcEtQhpu6QKyCRrX9PZSe+B+3ykmSnX0bbb\nS0pEHQ5Vu+58z2XQdM29p6dxP+mJExsXv3C0TFWHqiTHK/18jkKXSKWR3dc5Qw8yFEnS3W22CEiN\nh485VHfsoL+XX6a2u/9+WryyaRPtA/LWWxSC9+ab9DymT6fyn32WXvb+fmD2bFoEc9NNwK23Asce\nC8ybR2R98cVkE1+yhOTXrgXmzzdrFXxk7A5SHPkRI3f7GaQ4sGL2464uapcUua4uOXrCJgIfQfq0\nZPdLP75y3ffYzv/NN6uTrGS+8TlUfYNPXQ5VNrfYA5ikdQNyndx+CoTbwh4seJbga7PhYVnDT30f\nizpU3bqXRVPJfWCAHC5uhMjmzeZlqsOhGtLcbaeuuzeIj7Ql+5pvIPGZNlxTUaq2au/bPTREA8Fr\nr1H6Zz8jm/XgINmwh4YonnrLFiLOJ5+kdt+4kTrO9Olkm37uOSLpO+6gsL4TTyQH5PLlFIs9fz5t\nf3rffWTHtp2N7tTUjW7xhYOFno9NLO7AzteESJvbs6dndLmSFsqbQPGWrL7OaQ8+7kzKlps4Mc0G\nGwpjdOtnv2c+ErPlJGLxmT0k0kmJWomtULWdrTGHqo+MfQNEGdIuKycROGv4qTM6qSz2BfjCM906\nsbO1LJpO7i6JuYTGRBIid5c8WSN3ySOFtOvSyKXpd8ih6namHTvo2sFB4Ne/JkL+618pWkRr+oAA\nR8FwdMmLL9L//f3kcHz6aQoJvPdeKvP884FPfIJCA1etIrn586mN999/dISIb8CxPfy2XKjNQgNY\nSlifvWBq6tTR5OnWaWCgMUIkpFFNmUKdNeZQ5RWGPrswn3P3LvERqb1i0e7stra6eXPj8naf6YAj\nTqS6u4TBKy9dMrYH1YEBEznj1p1nLmyWYTtzjNC4XUJOY34eQJyMJ05Mk+N3p7u7caGRROAs5/NT\npA4q7Gvj0NuQucV+3u6CQLssngWXRVPJHYiTe8ih6nOUTprUOHC4L5PPBl1Eg5ZskpLjzPXc79pl\n7Lvr11OkyKuv0j4hg4PAV79K9ujt20nT3raN7Nvbt1P+hx9ONu4VK2jHvYEBCus76yzStB99lMj/\nkEOozBkzwn4A29af4lAtEvkRajOpbUMafkiz8RELX8Maj2TaGGuHqm9QsU07Ptuqr1wgXCde5dnR\n4TdT8HvrG1QHBsxKSd/zZnLv708j976+eHQLDxbcFu6AY6enTjVt4TOjuIOFz6E6OEi+uljbDg6a\nRWxdXY1at93vmYiVov+1NnJcBv/afgU+Zqft364KDN1UcmcnXmgVpbuIaWiocYWq/TJJhGF/nNju\n4DYBpUYn+GymdrlDQ+SEfOUVeuHuvJPMInffTVrDW28RWe/aRZs98cOaOZOuX76c7OEvvUR7X597\nLmnb69fTIHDkkUTEs2aZqBi37lrLMwaJgFwCD7XFrl3ysnCf09gXMsnEUkTDD5GnpBnHBgF3sBgL\nh6p0HzY5+Ug75PRLIXf7nnznpLpv3izL1elQdevumopcOdfZGmoz95w9sAFGS+a2dddx2Pnt3m3O\nhYiZf+00z8AkeelaHvAkObvu9nVl0FRynzjRaKQ2gfuiZThqpaPD7J4Yc6j69j9x7bihRRRMVDt2\nmC0THn6YyPuRR+h30ybaHnXrVuCuu0hbGBw09T75ZDq3YQOF+X384xTy99prtBhn8WJa4n7QQeT4\nHB5ONwe5xAKMlmPCcUk2VVvdsmX0NzpDdUqJbuHnwytFJYL0aeTd3eHIj6IOrI6O4o4uV66rq7FO\ntm0aGB354RI9IJOxa7dmuVQN33a2ujNJOwbcp9Xys5Ls57ZJLmTWS7Wlh0jb7cM28dn1YzKWiJR/\nQ4Rrt7GdTskvJp9yLDaQlEVTyX3uXDIx+EISudO5Wl4Rh6r0AvlMAp2dpCFv20YkvG4dEe2NN9JH\ndJ94wswE1qyhl6i/HzjgAHrxrruOHJEXX0yE9cgjFJ+9YQN97PbZZ+mayZPjIZN8H244mDT4uB1D\n69Fy/MK4JBvS5IrIcd5uXH/MoSrdl21GCZkOQpqxjzx898g2d5+crUHGHKopRBUitBSnX8jcEtOS\n7edhk0hoQJRs+LZD1Ta3SIMAa8xcd96nh9uQ0/ylJK4T/8YIcjzJFT0maeeSHA/EZdHUFaqnnkpk\nkLI/i0vaqStP+dy2bUQ6mzcDv/kNaVj33EOLbO64w3xx5uabyXyyeTOtQOzpoT2w3/1uChdcsYK0\n60suIVI/4QQKE+zqoogSJm2JSH0O1Zg2mGJiiHVOn1br05JD5drOS1d7S72P0P2XucfYIOCag0Jy\nsUFFcpS6JGsTKSCX5aYluVTTRiiEzmeycAcL34xu+/ZGR27I6eeaM+omzyoaeVE5tyzJfBKqk3TM\nJfIi13K6LCqRu1Lqo0qpDUqpYaXUUTH597/fOH5C03n7nO28tKeBnZ0UTTI4SFEizz9Pmzrddhsd\nO+ccIvK1a+l31y5arj5zJn3d/CtfIQK/9lraue/MM6l+73gHfSJMqcYX16dpu4QmaUA+k4XPFGF3\nXJ+cpK0WJW1fuTE5H/GFZh1u5EdKWe5AkupQlcxBobrb+YUGlc5O2aQkkXGIPOt2qNqDiqvhp0Rj\nsGnDJqNUAgL89uNUknXLdAkwtU6hcu32kTTnUJ0kuZRryxC5VEZZVDXLrAdwDoDvpAhPmEC/Gzca\nW25o47ChIXIqsu1bKTKDvPACaRerVtGL/NRTJLPPPrTY5i9/ocU3P/4xLbw58EAi8fPOow8q9PY2\nbnDkkrG0iMnnHEw1t9gDlts5be2f6/Sud43OTypryhTqmBzyZc8YbLne3tFkLGmr9ocDbBusS6Qc\nMWFP0yWimjEjTFTuNbzAx97NUXoGKQOOS9p2WWUcqr5YZ1eLl+5RInqWcwdz911iOYmou7uNEsJy\nkgadSrj8HCU59iW413K7pJQVIk3fudS6S1oy/5ZpC/dYqE72oFGGyN0y7DqXRSVy11o/DwBK8d58\nYQwO0ouwejWF7g0PU2jg9u1ExOvWkf31978n0rn6aiL7KVOIdDo76UvofX0UPbJ4MYUFfulLtLf2\n/vsDRx9NC3HslZKhWYJLHj4bvjQ1T4lf37mz8eMIHHMcsnHapOiG9UnaKhB3qMaIzz0XskdLDjvA\nT4qS6cBnRok5Gzk/Dul0idV3j25ZRVaoFo1aYTk7xtquO5M0y9nlupEaMcJIJaoUuZB2W5W8bEdu\n2bqH6lTEfOMzt/DzsY/FtO/YIJRyLNa2ZdFUmzt3hPe9j0wpv/0trbIcHiYSmzaNvkDzmc/QA1i5\nEjj+eJK/8EIi3YMPJjnJKesj7RC5Sxp0yorSouYWn1bmI20fyUomBql+Pj9AirNRKte1M8e0UCk/\n+wX2mXZig4Cdlhb4+OoeKksyt/jMHu6zC5k9qmjQEgHwb2p+KaaIIqaDGLmmHIuVW/Qey8i5x+og\n6NRjvmcRurYsopq7UuoRADPtQwA0gBVa6weLFLZy5TXQmjrS/PlL0dOzFFdeSREm55xDH/AdHiby\ntvdn5hHf7mjd3fVuxSpp5D7ik/KTiFQacPicHT0QGiBSNG07zp3lfJp7im3eHojsJdPuOTsN+MnT\nPg74z9n1c/OT5Oz8XPMIp4eGSMZdOchmGTYpsaxL1N3djWYPH8k2U9OOadNSHiFtNfWYzyxT9Zh9\nD66GX0QjT20L3zFfWXUOiCl1Wb16Nf7859XYto22FSmLKLlrrU8rn30jzjjjGjz6KO3Rfe+9jdN7\n7qx2/HFXl1mhGiIWl0hdLXnrVhOzHYqeKGpumTAhrPHFloW7+blyTIIxMuZQSFuO0zyQpJZrDxyS\nturmkappDwzQIMxyvJ2re83AgFmr4BtwOG/24fjKZfONXT97JaEdhsfvUBmtUersMU071ulTCc2X\nXxFzQqq5ZazIvcwx6R5D/oKiZqYYGadcW2QQsNNLly7FoYcuxcaNRO7XXnstyqBOs0zU7p5iYpAI\nKIUgi67K5IcrOfPcT3nZZOeSoktAqfHh0iBga55F77FoFAwPnLGYciD8fGKatntPrlxIw3frzku/\n+R5tOSnaQ9K4+H8feTFBhOSkcyG5WGfnc3WXOx6PxZ5Pihwfcx25bjtWJdlQ3WPPr2oZdp3Lomoo\n5NlKqZcAHAfgIaXUz0PyEkFKZBdzdPnkpEHAJd8Q8fm0ZJeM3E2CQnXy3WMoUsNnj5YGxKKkLYX1\nhZ4BIGvNbn4slxJj7eZh78xpb5bU3S3brVNNID5NLrWz+zpsaPl4LN8UOYmgQvcYM7ekmgQkrbZu\ngixzLKR1p5g9qpijisyAqpQRes5lUYnctdY/0Vrvr7WeqLXeR2v9oZC8S0Axk4BNfNI5iYwluZBW\n6yPSkEPVDqeTyI7LtSM6fAQpDThFo1t8M4aQXGwLU8mWLrW7K+d+sILz2707TlTSeUkuhSBTCSP1\nWqkTx47x/3XcY10kG6tfyrUxkg0dK1OuNEinvgP2u+6Ti80SpGdR9hn4Bu7QtWXR1GgZn0bu01Zd\nuRTiKytn10HaijVmsvBpv/b0McVU5OZXJbrFJyeZUdwZjkvavhhr3/Lx2BRWOiedD8m58r5pfWrH\nLto5fffLv3WE/4XaNkSeUp1aqX2nPgOpzkXbp4ycUqOfla9ORe4xlchD70pZjAtyt8nJJRlXgy7j\nUJXs1hJ5FonFlgjSTtv3wXKhxSqxDZ3cNpNMIJy3vVLUZyrhOrnfhuWyON6a5fg3puXwb4xQpHPS\neSm/IpEfsSlvUTnftUXl3E7v3k/Z6JZUDb+otjoWx1IdtNK5onKxATHl3av7GcTy4HzKouLlxeCa\nWyS7sLRbnx0K6bNvhyJTJC05lB8w+rqYyYI7SSg/e0/v0EzAFyEird7s7m4MheQprB3W59uoqS6t\nqKimbf/GOhHnn0pAITn7HC+7q6tzSnVOvZYHrFaRbOqxuqNleEab8k6lmmVcs0dZM5jvnG973yrK\ngS8PPl4WTSd3rrRvYYy0cjCF3N2ZgBRm2dk5Wqt1TSopGnnIiShp01J+qTMGu9wdOxrDCbn9gNGk\nvWtXutMvVYMuo2lz3aVy3c6c2sFC+aZcy+/ZWJbhI6rUWQT/2vLS1L3O8qVyY0RZpS78f+o7EKuf\nVH7qO5XSN6T6xZSTonu3S3Jl0ZIVqkxMIRu5a4pwCThVziZj3/4nbhoYTawxDV+quy8/N+1Gptj3\n4UaLxEwl9jHXkVSnFlPG3BIjh46OuPMr1tFSySbWsUJ1ryJX1g5fRAst0hZV2sxHsqzhhwaGVDIe\nK9L2ycXe0dRnwGG1Uv5F+2FZtMTmzpplsx2qgP9cEY1cckpK5huJtPkcR4+kPuiyWlvsBQ7Fdlcx\nt5TpsPxbt4ZYR5tJnbPMfaS0TyoB+ciz1WaeotEtsZlFHW0We1dCg0+VviHln9of7DzKoKUOVV9c\nehlHacpMAGj83y7XTrNcqkZuOyXtbVT5vP3re9Ap0+W6jqWSl9thOzvllz40myii4afUs+g9ppoT\nqpgdipRRpH3qIiqpvVmm6pYEddZ9LOViA2LR/EIDmJ13rO48EMfyKIOmm2VcMwrfkM8e7XOU2po1\nT+ft/U9SSTtkI3cdsSxXJPyvLu3Xd6zK9LsMuabU3c3X3uAr1Ga+tvMRVlEyLmpaKdI+deyJEpsR\nFCWg2H13dxfT8JsZ3VJWzp1ZjdVMgPOXynXlQmW5+fnyKIOW29xtDVpKp6xQBfznYnL2Ck13X2z+\nLdpxinRsqcPwrzRVTCGbIqSU2rGLaL++ZeGhOtVlRkmtZ93tI8lL9yblkUJUZQgoVj+fKcIlSI4B\nD0W3pNYp9E5LMwupLdy0K1+lfqEBzK6f1IbSNbG2SNHwy6Ki4l8MKWaU3bvlUEjbJOKmgdHn7Jhy\nex+Sri55H5KyHcJHGKkvldthYlqjr3z32K5dademkFwRrT+kKdnnqvoQqs5EeCpcVJsv8mzr1BrH\nUgtNlXPrESLZmAbt61MhuWbco2tuceueam5x+6UvP5+ce64MmkruKUvubXLv6hr96TY37ttujG3b\n0h9ujIxjL1VRbSym0ZUl3u5uGqzcYxK5p9YlldDKTtNDMxuJCGy7cF1kHIpzd+/RnolIhFb0HSir\nCEj5SG1WNL/Yc0whILutUtoitc2aISdp0FU0clcu9HxS8iuLlpB7dzcRs70ggM8NDzfuf9LZ2Whz\nL/sR3lhHTNGKiphbUjX8UL4horKPMbnbx5jci5o7fBqidC1rOaH6+covooUW2Rsk9X6c+YvlAAAO\nk0lEQVRccq+iGZch7dgMLaXcGFGFTCD2NTFzSxlCi9XdHlCK3mMZuVh4Zp33mGpucfObMMFsd93b\nC0ydikpoulmGX7Dt2xsX5DCBS4tLtB69Kix1mlf0JUjtEKEOy/chrWTzEVusfr5j7ktkt2foPnzH\nfOWW1VaLtrevk6S2d+r9pK5QreMeY3I+00Zd5Up9w3727v9lCK2o3Fhp+D5zUCw8s8w99vQYMn7n\nO026v5++KMf48IfpE6AA8IEPkPIK0Ledv/c9I7dyJR0DgIsuIkW3CppK7vZ3TQH5wTC5S5pp7OFW\nWSBSR8cJHYvVKfSSptQzRS5GkEXrLN13SqdLbUdfm5aV43NsPkj9sHOrNPKxNrf4TBE+gozVvQxB\n+t6fri4zcw/VKXVABMzX3QD6IM+kSZSePp2+Wcz44AfpA/UAcPrpJr1gAfAv/2LkHngAmD2b0ldd\nZY7PmgU8aH2j7qabTPqYY0y6s5PKYnBefI7vvSyaSu4ANW4oNG737vgx9yXgc6lkKMVsp3RYHxmG\nyqr7mK9891gVh6rbPj6tX8qjrGmjqlxoJufK2c7UMu9A0TrVMdspYzYaC3NLiGTtQaCMnK+vAUTM\nPNOfMoX+AHpPDz/c1PHEE4H3vpfSCxeafgDQ95rZ1HH99eZLXkcdBdx/v5H74Q9N+owzTHrSpMb/\nuRyu33hDS8gd8JNhCrmXIUOfrc13ra3l2edc211qnWIafkgDiRFLyOZepi5SWXWScR0afhFtugyh\n8SxQqnvZAcz3Dhate+ydLkPaYynnasy2LXn2bGDvvSm9YIHpX5MnA9/6lpG74QbggAMovWwZcOqp\nlFYKePppI3fFFSa9337ARz9q/j/xRJNmbbyd0TJy9xEL25nqIveig4B7zI2UkLQhPlfnsVjdY2ar\nKg5V6VxRMo5Nl+vQ8KW8XVNEKvGVIUhbASha9yqDRcz05GrGdZL2pEn0YXgAmDnTmDZmzgTOOsvI\nffrTRrO94AKKZAOIVJ95xsjddZdJL1pEfwANrJdeas7Nn994T1WdjXsCmk7uMROH61Dt6jLknmqC\nKUNoqddK55j0fQNErPwqx1x7KP9KDlUprC9E2tJ91K3hpxJVjCA5XYcpwiXIUH6umaeOtogNsEWd\ng5yfrUFPndpIkIsWke0ZIK2Y5ebNA775TSP30EPAvvtS+rLLzPGpU4HbbjP/L19u0qyZM6ZNQ0YT\nMC40d/vllhyqIZt7KmnXRa7uuSL2aGmTriJ1D5FxqpxETilyZdvCR9SpbZYyCNj3EcsjlbTtdi1D\npO41tnMwtW1Ds4lJk4ydd8YMQ8xz5gAnnWTkbryRtGoA+NznjNa9YAHw8MNG7t57TZpNHlzuaaeZ\n/2fNQsbbBOOK3CUi9x2TOn1MCwzJFfnosU+rrUKQdh6SwzdWlntMsrlLdU+tUyoZ1zlYSM8vNgiw\nTJ1yobZw5aS2sNvclbPNHF1dZH7g6047DTjsMEoff7xxIgLAk0+a/+1ojLlzgWuvNf9/8pMmzeF4\nDHYoZrQnWkbuPi2vbodqiNxig0FIrqj5JlXDLzITABrDpWw5H7n76h4izyIzgbJtYZ+T6iAtC08l\n2TrlQtr0xImGqPfZx5Dp5MnA5ZcbuauuMvboZcuApUvNPW7YYOQuucSkZ8xo1KjnzEFGRhAttblL\n5FW3QzXV5s31CF1bxtkYq0uR/CQ5n73XdagWrbt7zld3H9nFyi3aFrE8ypCxT85eKQgQUff1Ufro\no41pYsYM4Pbbjdzdd5t46bPPpj/O++tfN3InnGDSPT206CUjo2601CwjkWfdDtWU/U98+RbVVstq\n+GUGgZj92Nbcy3wVpqpzsEqbpd6jVFf3nB0f3dtLf4xTTjHOwWXLaA8jgFYJ/uIXRm7VKpO27dkd\nHY2hdpxXRsZ4wLi1uac6VGOEX8c2pTGSDQ0WRQktJJdqikjVVt20G9tdxnxTVMOXtmlIuUd7heE+\n+1BUB+PKK412fdllJu85c4Bf/crI3XqrSR9xhEkrBRx8MDIy3tYYt+QeOxYjZvc35lCV8qhT05by\ni80YpO0UUknb1txtuSKRH6l1rzJwAo0henvtZbRrN/LjlluM2eOLXzQa+dy5wA9+YOQ+/3mT5r06\nGLa5JSOjndH2DlXpXFGNXNJW63CoppJsmcgPH7mHNPIUEwhfn7Lxkx3xs99+JiTvsMOAc84xZf3y\nl7TxEgB87WvGlHbQQcA11xi5Cy4w6b32QkZGRgDjzqGaSu5FNgnzRVnEzB6t0PBD9uMimrbtUE0l\nbV9+dh49PRT9wZg3z5hAFi0yC2M6OoD1642Z5xvfMNccdBD9MZYsaSw3IyOjOsadQ1Wyr0t2eP4/\nhbRjZNxsh2pspWiqEzHVoZpC2pMmNZL2woUmlO+UU4BDDqH05MlE2oybbzbpuXPpj2HbwTMyMpqL\ncWtzjzlU+f+YQzVFrqhG7tO0pXJTNe3UcL2imvbkyY2kvWiR0bTPPBN46y1KT58OPPGEkbNJ29W0\n3cUwGRkZ4w/jhtx5Vab0sQ7JLCPlUVTTLivn2rfLxGz7okdSBwE7WgSgFYy8h8e559InDAHaSW/N\nGiNnR4jYHxQAjE08IyPj7Y9xQ+6cljZFssm9zB4dVeViC3JSNfwUc8uECY3Lwo8+2my0dN55JkJk\n3jzan5ph767HS9YZeZFMRsaeh5Y5VG2yizkRbXIPyUkzgZCcZMaR6pc6CLj709h52ItpAIqj5giR\nZcvMxk8LF6LhYwE//alJH3usSStF8d0ZGRkZEsaN5m5rtS7J2nb4FHJ35XyOV/6tayYAUBw2x1Kf\neqoxfSxZ0rj16eOPm7T9qa3OTuA970FGRkZGJTSV3EPfrEzV3EODgGQHd/+X0ikx2z095rqjjjLf\nO1yyxOw7AgB/+pMh9wsvNMcnTqRtVjMyMjKagaaSe4g86zbLlJGz98U+9ljjsOzsBJ5/3gwCK1aY\na/r7Gz/flVdAZmRkjAdUInel1L8BOAvADgB/BHCR1vrv3sJKLn2vSu72Vqzz5pnP0k2d2mjTXrnS\npA85xMR2A6OXsWdkZGSMZ3RUvP6XAOZrrY8A8AKAL4SEbZNI2c2tfHLTppmokhkzGvcX+f73gcWL\nKX3SSeYjuko1fmWmWVi9enXzCx2nyG1hkNvCILdFdVQid631f2qtR/RqPAZgv5B8yAHqs6XbUSb9\n/cBxxxm56683duxPfQr48pfNNVddZeT23bfxoxatRn5xDXJbGOS2MMhtUR112twvBvCjkICP3KdM\naVxFef75ZhXkpZcas0xfX+O3Hu0ok4yMjIwMgyi5K6UeAWCvXVQANIAVWusHR2RWANiptb4nlFd/\nv9lY6swzzWKdQw+lnQEZ111n0rzFa0ZGRkZGOpRm72LZDJT6BwCfAnCy1npHQK5aQRkZGRl7KLTW\nqug1VaNlTgfweQAnhIgdKFe5jIyMjIxyqKS5K6VeAPAOAK+PHHpMa3154JKMjIyMjCagslkmIyMj\nI2P8oWqc+ygopU5XSv1BKbVRKfXPHplvKqVeUEo9rZQ6QpJpB8TaQin1MaXUMyN/a5RSh0n5vN2R\n8k6MyL1PKbVTKXVuM+vXTCT2j6VKqaeUUhuUUo82u47NQkL/mKqUemCEJ9aP+PfaEkqpO5RSryml\n1gVkivGm1rq2P9Bg8SKAAwF0A3gawKGOzIcArBpJHwsy5dRaj/Hwl9gWxwHoHUmf3o5tkdIOltx/\nAXgIwLmtrncL34leAM8CmDXy//RW17uFbfEFADdwO4DMv12trvsYtcdiAEcAWOc5X5g369bcjwHw\ngtZ6k9Z6Jyju/SOOzEcAfB8AtNa/BdCrlGrHz0RE20Jr/ZjWeuRbSHgMwKwm17EZSHknAOBKAP8B\n4P+aWbkmI6UtPgbgPq31KwCgtf5bk+vYLKS0hQYwsjE23gngda31ribWsWnQWq8B8EZApDBv1k3u\nswC8ZP3/MkYTlivziiDTDkhpCxufBPDzMa1RaxBtB6XUvgDO1lrfBlpH0a5IeScOBjBNKfWoUup3\nSqnlTatdc5HSFrcAeK9S6i8AngHwmSbVbTyiMG82fT/3jNFQSp0E4CLQ1GxPxL8DsG2u7UzwMXQB\nOArAyQAmA1irlFqrtX6xtdVqCT4I4Cmt9clKqTkAHlFKLdBab211xd4OqJvcXwFwgPX/fiPHXJn9\nIzLtgJS2gFJqAYDvAjhdax2alr1dkdIOCwH8SCmlQLbVDymldmqtH2hSHZuFlLZ4GcDftNbbAWxX\nSv03gMNB9ul2QkpbXATgBgDQWv9RKfUnAIcC+H1Taji+UJg36zbL/A7AXKXUgUqpdwC4AIDbQR8A\ncCEAKKWOA/Cm1vq1musxHhBtC6XUAQDuA7Bca/3HFtSxGYi2g9Z69sjfe0B298vbkNiBtP7xUwCL\nlVKdSqlJIOfZc02uZzOQ0habAJwKACP25YMB/G9Ta9lcKPhnrYV5s1bNXWs9rJT6J9BWwB0A7tBa\nP6eUupRO6+9qrX+mlDpDKfUigG2g0bntkNIWAL4IYBqAb49orTu11se0rtb1I7EdGi5peiWbhMT+\n8Qel1MMA1gEYBvBdrfX/tLDaY4LE9+IrAO62wgOv0lpvblGVxxRKqXsALAXQp5T6M4AvgRaIlubN\nvIgpIyMjow1R+yKmjIyMjIzWI5N7RkZGRhsik3tGRkZGGyKTe0ZGRkYbIpN7RkZGRhsik3tGRkZG\nGyKTe8YeDaVUr1LqH1tdj4yMupHJPWNPx7sA5K+HZbQdMrln7Om4AcBspdSTSql/bXVlMjLqQl6h\nmrFHQyl1IIAHtdYLWl2XjIw6kTX3jIyMjDZEJveMjIyMNkQm94w9HVtgPuWWkdE2yOSesUdjZAvZ\nXyul1mWHakY7ITtUMzIyMtoQWXPPyMjIaENkcs/IyMhoQ2Ryz8jIyGhDZHLPyMjIaENkcs/IyMho\nQ2Ryz8jIyGhDZHLPyMjIaENkcs/IyMhoQ/w/t+ukiQqreNMAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def forwardEuler(h=1/100, tstart = 0.0, tend = 1.0):\n", " n = int((tend - tstart) / h)\n", " ts = np.linspace(tstart, tend, n+1)\n", " xs = np.zeros(n+1); ys = np.zeros(n+1);\n", " xs[0] = 1.0; ys[0] = 0.0; # Set the initial condition here\n", " for k in range(n):\n", " xs[k+1] = xs[k] + h*( 398*xs[k] + 798*ys[k])\n", " ys[k+1] = ys[k] + h*(-399*xs[k] - 799*ys[k])\n", " return ts, xs, ys\n", "\n", "# Choosing h = 1/200 leads to a numerical solution with a spurious oscillation, as shown\n", "# Choosing a larger step size leads to exponential growth (instability).\n", "ts,xs,ys = forwardEuler(h=1/200)\n", "plt.xlabel(\"t\"); \n", "plt.plot(ts, xs, 'r-', ts, ys, 'b-');\n", "plt.plot(texact, xexact, 'r-', texact, yexact, 'b-');" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEACAYAAABVtcpZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmUVeWV/vHvZlIpEZlkHlQmmUSZjUohDjhEkTbRaEdj\na2J+ia0ZOh3TaSPplY4xK51lJ3a3MTF222mnxBkHcCoMRhEHHBAUZEYsRCYZCopi//7YVRaSKii4\nt+45dc/zWauWVbcO577eBec5Z7+TuTsiIpJNzZJugIiIJEchICKSYQoBEZEMUwiIiGSYQkBEJMMU\nAiIiGZZzCJhZDzN71szmmdlbZnZNPcf9yswWmtlcMxue6/uKiEjuWuThHDuB77j7XDM7FHjVzGa4\n+4KaA8zsTOBod+9nZmOAW4GxeXhvERHJQc5PAu7+obvPrf5+MzAf6L7HYecBd1YfMxtoa2adc31v\nERHJTV77BMysDzAcmL3Hr7oDK3b7eRV/HRQiIlJgeQuB6lLQn4Brq58IREQk5fLRJ4CZtSAC4H/d\n/eE6DlkF9Nzt5x7Vr9V1Li1mJCKyn9zdDuTP5etJ4PfAO+7+7/X8/hHgUgAzGwtscPfy+k7m7viO\nHfif/4z/8pf4Cy/gu3bF6xn6uuGGGxJvQxq+9Dnos9BnsfevXOT8JGBmnwMuAd4ys9cBB/4J6B3X\nc7/N3R83s7PMbBGwBbh8nyf+/vfhiSfgpJPgtttg+3b4+tfhiiugQ4dcmy0iIuQhBNz9BaB5A467\ner9OXFYGd9wB/fvDNdfA//0f/OAH8E//BMceC//yL3D22QfYahERgbTOGN66Fd59F4YOhdGjoX17\n2LwZNm6Ee++Fli3h3HOhXTv4h3+I44tQaWlp0k1IBX0OtfRZ1NJnkR+Waz0p38zMfdYs+Na34Gc/\ni7LQK6/89YEbN8IPfwh33glbtkTZ6N/+DUaMKHyjRUQSZGZ4wh3D+TV7djwB3HUXXHxx3ce0bQu3\n3FL7dPDRRzBqFPTsCTffDFVVhW2ziEgTlM4QePllOP54ePBBuPDCvR9rBhdcAPPmwbJlMG5cPD2U\nlMBFF8GqOkeiiogIaQ2B2bOjzn/ssdB9PyYW9+wJ990Xf/anP4VZs+K1oUPh4bqmL4iIZFs6Q2D9\nenj7bZg8+cD+fPPm8J3vwMqV8VTRvj1MmQKHHw7f/nbRdiSLiOyvdIbA4MFxAT/yyNzPNXIkzJwZ\nfQdf+UoMO23TBkpL4dVXcz+/iEgTls4QKCmJWn6PHvk756GHRofxhg3wxz/C2rWf7UjetSt/7yUi\n0kSkMwRatYoQ2J/+gP0xZUqUm5Yvr+1IPuSQeH3x4sZ5TxGRFEpnCLRsCZs2QadOjfs+PXpER/K2\nbfCLX8Brr0HfvnDUUfDrX2uYqYgUvXSGQFUVdOkCzQrUvGbN4O//HpYuhQULYNiwmIl8yCExM3nB\ngn2eQkSkKUpnCOzc2XiloH3p3x8eeihGEN18c8w/OOYY6NUrnhZ27kymXSIijSCdIbBjR3IhUKN5\nc/jGN+D992HJEhg7Fq6/Pp4OJk2Ct95Ktn0iInmQzhDYvj35ENhdnz61k9B+85sIhZqJbD/5SYSW\niEgTlM4Q2LYtXSFQwwz+7u9ihdPly2HCBLjxRmjdOhawe/bZpFsoIrJf0hkCW7emMwR216MH/OEP\nscT13XfHSqannhqzkq+8Ej78MOkWiojsUzpDYMuW9IdADTP4whdieOmmTfC1r8Gjj0K3btCvX6x0\nqqGmIpJS6QyBTz5pOiGwu0MPhZ//HMrL4fXXYcAA+O53ozP5tNNiHSMRkRRJZwhs2tQ0Q2B3xx4L\n06ZBRQXceit88EGMMOrYMbbLXL8+6RaKiKQ0BA46KO6ei0FNZ/K8ebHxzUUXxX7JHTrAkCHw+99r\n3SIRSUw6Q6Bdu6Rb0Dg6dIg+go8/jr0OOneGr38dDj4YJk6E559PuoUikjHpDIHDDku6BY3vhBPg\nmWdiOOx//iesWRPLW7dtC1/+shayE5GCSGcIFEspqCGaN48hpW+9FXseXH01PPccHH109Iv84AfR\nUS4i0gjSGQIHH5x0C5LRpg3867/GhjqLF8Mpp8RTQtu20X/wm99ouKmI5JVCIK2OPBL+93/j6aCs\nDLp2jVFFBx8cZaOnn066hSJSBNIZAlkqBzXEySfDU09F/8FvfxvDS08/PeYlTJkCc+cm3UIRaaLS\nGQKtWyfdgnRq1iz2SX7jjegn+Md/jO+PPz5GVF16KSxalHQrRaQJUQg0VSUl8KMfxVLXa9fCFVfE\nAnb9+sXQ029+E1avTrqVIpJyeQkBM7vdzMrN7M16fj/ezDaY2WvVX/+81xMqBPZP+/ax4c3KlbG6\n6eTJ8Mc/xvpFPXvCdddF34KIyB7y9SRwB3DGPo553t2Pr/76yV6PLCnJU7MyqGfPGEW0Zg3Mnx/9\nCbfeGuWifv3gpz+NvgUREfIUAu4+C9jXYjjW4BMqBPJj4MBYomLDBpg9GwYNiiGoJSUx5PSXv4wN\nfEQkswrZJzDOzOaa2WNmNmivRx56aIGalCGjRsHDD8cy3TNmRKnohz+M0tvgwVFOqqhIupUiUmCF\nCoFXgV7uPhy4BXhor0crBBrXqadGEGzdCtOnRwnpRz+KQDjmmFgOe+vWpFspIgVg7p6fE5n1Bh51\n92ENOHYJMMLd19XxO7/h/PNhWJymtLSU0tLSvLRR9sI9lqv4xS9g5szoN+jfHy67DK69Vp31IilS\nVlZGWVnZpz//+Mc/xt0bXnLfTT5DoA8RAkPr+F1ndy+v/n40cJ+796nnPO4PPwznnpuXdskBmjkz\nngjKyiIQ+vWrDQT12YikipkdcAjka4joXcBfgP5mttzMLjezq8zsa9WHXGBmb5vZ68DNwIV7PaHK\nQckbPx4eeyz6EGbOjKeCn/401jfq1w+uvx7W/dWDnIg0MXl7EsgXM3P/y19g3LikmyJ1+ctf4gnh\nuedi1nK3bjEv4Xvfg969k26dSCYl/iSQd1pALr1OOAEeeigmn731VnQy33cf9OkDnTrB3/4tvFnn\nnEERSaF0hoAWkGsaBg+G//7vmJi2YgVccgn8+c8wfHhsDPT5z8cTQ8qeNkWkVjrLQUuXqrTQlG3c\nCDffDHffDQsXQqtWMGYMfOMbcMEFsRCeiORNLuWgdIZAeTkccUTSTZF82L49lr++447aMtGQIXDx\nxbG/cps2ybZPpAgUXwhs3JiNfYazpqoK7r8fbrsNXnwxhp726AHnnAPf+Q707Zt0C0WapOILgR07\noGXLpJsije3116NsNH06lJdH8J98cuyzfPrpYAf0d1okc4ovBFLWJimAtWvhV7+KkUYLF0Lz5nDs\nsbFRzpVXarCAyF4oBKS4VFZGp/Lvfgdz5kS/Qu/ecN558O1va9CAyB4UAlK83GMZ7JtvhmeeiSeG\nww+HE0+Er30Nzj5bo40k8xQCkh0ffAC//nVMWFu4MF4bMACmTIktNbt0SbZ9IglQCEg2VVXFHgm/\n+10sZ7FxY2y1edJJcNVVcMYZekqQTFAIiAAsWwb/8R8RDIsWRQAccwz8zd/ERLVOnZJuoUijUAiI\n7GnnTnjgAbj9dnjpJdi0CTp0gNLSeEqYOFFPCVI0FAIi+/L++3DLLTBtGixeHAHQv3/sW3HVVbEA\nnkgTpRAQ2R+VlfGUcOedMXN5/fpYvmLkyFjO4uKLtZOaNCkKAZFcrFkT6xs98ADMmwc7dkDXrjBh\nAnz1qzGLWbOXJcUUAiL54g6vvQa/+Q3MmBFLZDdrBgMH1paOevVKupUin6EQEGksO3bUlo5eeqm2\ndDRqFFx4YZSOtB2qJEwhIFIo5eUxL2H30lHHjrEd6iWXxFabrVol3UrJGIWASBLc4Z13oj9hxoyY\nwVxVFfsun3wyXH45nHJKLIYn0ogUAiJp4B4zl3//+9hWc9myeP3II2NewhVXRBlJncySZwoBkTTa\nuTOeEP7nf2DWLFi9Glq0iPkJkybFEtkDBybdSikCCgGRpqCiIvoS7rorOpnXrYv+gwEDYhOdr3wF\nBg3Sk4LsN4WASFO0aRPce29sufnqq/Dxx7GjXr9+cOqpcNllMHy4QkH2SSEgUgw2b4Y//Qn++MfY\nTGft2igf9e0bHcyXXRazmhUKsgeFgEgx2ro19k24997YWGfNmhhpdNRRMZv50kth7FgthCcKAZFM\nqKiARx+Fe+6JPoXVqyMA+vSJndYuvDDKSC1bJt1SKTCFgEgW7dgBjz8eTwovvggrV8KuXdC5M4wY\nERPXvvAFaNs26ZZKI1MIiEgEwOzZcPfdMU9h0aJ4ejjsMBgyJIalXnJJlJOkqCQeAmZ2O3AOUO7u\nw+o55lfAmcAW4CvuPree4xQCIvmyeHEMSX3ySXj77RiR1KpVdDaXlsbaR+pXaPLSEAInApuBO+sK\nATM7E7ja3c82szHAv7v72HrOpRAQaSwbN8aQ1IceihFI5eUx2qhnzwiDyZPhnHO0KF4Tk3gIVDei\nN/BoPSFwK/Ccu99b/fN8oNTdy+s4ViEgUiiVlfD003DffTGredmyeK1t29ifeeJEuOgiGDxYQ1NT\nLJcQKNQzYHdgxW4/r6p+TUSS1LIlnHkm3HFHLIC3fTu8+y5897vxu//6Lxg2LEpIRx4JX/wi/OEP\nUVaSotAi6QbUZerUqZ9+X1paSmlpaWJtEckUs1jb6Prr4wtg2zaYPj2WvHjpJXj44RiZdNhhsfbR\nhAkxPFWzmwumrKyMsrKyvJwrqXLQAmC8ykEiTZA7vP9+DE196qnocF63LjqXu3ePlVLPOSd2Ymvf\nPunWZkJa+gT6ECEwtI7fnQV8s7pjeCxwszqGRYpIRUWsmPrAAzFnYdmyKC21bh1lpHHjIhROPRUO\nOSTp1hadxEPAzO4CSoEOQDlwA9AKcHe/rfqYW4BJxBDRy939tXrOpRAQKQarVsGDD0Yp6Y034IMP\nYi5DTRnppJNgyhQYPVob7+Qo8RDIJ4WASJGq2Ynt/vvh2Wfj+7Vr43edOsWEtokTIxgGDFD/wn5Q\nCIhI01RZGZ3NDzwQQ1Tfew8++SSeDLp1g+OOixLS5z8PvXopGOqhEBCR4rFlS5SQHn00lsFYtixW\nVG3ZErp2jSGrEydGMBx1lIIBhYCIFLsNGyIYnnwyZjovXRph0aIFdOkCQ4fGUNVzz40hrhkLBoWA\niGTPpk0x2/nxxyMYFi+OYGjWLFZSHTIk1kc699yi37ZTISAiAhECzz4Ljz0GL78c8xk++SQC4Igj\nYvmLE0+Es8+O/oYWqZwvu98UAiIi9dm2DcrKIhhefBGWLInyEsDhh0e/wsiRcPrpsY3n4Ycn2twD\noRAQEdkfO3fC669HMLzwAsyfHyuq7twJBx8MPXpEP0NpKZx1Fhx9dKrLSQoBEZFcucfubE88EZvy\nzJ0LK1ZEial585jLMGBALLl91lkwZgwcdFDSrQYUAiIijWfLFpg5M0YnzZ4d/Qzr1kVotGkDvXvH\nsNWTT46SUu/eBX9qUAiIiBRSVRXMmwfTpkU/wzvvwOrV0f/QvDl06BAlpBEjYk7D+PHQrl2jNUch\nICKSBps3x8znp56CV16JPRo++ij6Glq1isluNSWlM86IkMhDSUkhICKSVu7xlPD00zFKae7cmOxW\nM0KppCS29xw6NBbVO+006Ndvv/Z9VgiIiDQ1VVWwYEH0NcyaFeWlVatiiQyILT579YpJbyecEGso\n1RMOCgERkWJRUQGvvgrPPBMT3hYsqO1vgDrDwQYOVAiIiBS1bdtqw2HOnNpw2LoVA4WAiEgmbduG\ntW6tEBARyapc+gQa3v0sIiJFRyEgIpJhCgERkQxTCIiIZJhCQEQkwxQCIiIZphAQEckwhYCISIYp\nBEREMkwhICKSYQoBEZEMy0sImNkkM1tgZu+Z2ffr+P14M9tgZq9Vf/1zPt5XRERy0yLXE5hZM+AW\nYCLwATDHzB529wV7HPq8u5+b6/uJiEj+5ONJYDSw0N2XuXslcA9wXh3HHdAKdyIi0njyEQLdgRW7\n/byy+rU9jTOzuWb2mJkNysP7iohIjnIuBzXQq0Avd99qZmcCDwH96zt46tSpn35fWlpKaWlpY7dP\nRKTJKCsro6ysLC/nynlTGTMbC0x190nVP18HuLvftJc/swQY4e7r6vidNpUREdkPSW8qMwfoa2a9\nzawVcBHwyB4N7Lzb96OJ8PmrABARkcLKuRzk7lVmdjUwgwiV2919vpldFb/224ALzOz/AZXANuDC\nXN9XRERypz2GRUSauKTLQSIi0kQpBEREMkwhICKSYQoBEZEMUwiIiGSYQkBEJMMUAiIiGaYQEBHJ\nMIWAiEiGKQRERDJMISAikmEKARGRDFMIiIhkmEJARCTDFAIiIhmmEBARyTCFgIhIhikEREQyTCEg\nIpJhCgERkQxTCIiIZJhCQEQkwxQCIiIZphAQEckwhYCISIYpBEREMkwhICKSYQoBEZEMy0sImNkk\nM1tgZu+Z2ffrOeZXZrbQzOaa2fB8vK+IiOQm5xAws2bALcAZwGDgS2Y2cI9jzgSOdvd+wFXArbm+\nr4iI5C4fTwKjgYXuvszdK4F7gPP2OOY84E4Ad58NtDWzznl4bxERyUE+QqA7sGK3n1dWv7a3Y1bV\ncYyIiBRYi6QbUJepU6d++n1paSmlpaWJtUVEJG3KysooKyvLy7nM3XM7gdlYYKq7T6r++TrA3f2m\n3Y65FXjO3e+t/nkBMN7dy+s4n+faJhGRLDEz3N0O5M/moxw0B+hrZr3NrBVwEfDIHsc8AlwKn4bG\nhroCQERECivncpC7V5nZ1cAMIlRud/f5ZnZV/Npvc/fHzewsM1sEbAEuz/V9RUQkdzmXg/JN5SAR\nkf2TdDlIRESaKIWAiEiGKQRERDJMISAikmEKARGRDFMIiIhkmEJARCTDFAIiIhmmEBARyTCFgIhI\nhikEREQyTCEgIpJhCgERkQxTCIiIZJhCQEQkwxQCIiIZphAQEckwhYCISIYpBEREMkwhICKSYQoB\nEZEMUwiIiGRYi6QbUBd3MEu6FSIiydu1Cz76CN57DxYtgqVLYeVKWL06Xl+3LrfzpzIEtm+Hgw9O\nuhUiIvnnDp98Au+/DwsWwJIlsGwZfPABrFkTF/WNG2HLlrgWVlXFn2veHA46CFq3hsMOg/btoWNH\nGDoUFi8+8PaYu+fn/yxPzMzXrXPatUu6JSIiDeMeF+8FC+KOffFiWLHisxf2Tz6BrVuhsrK22tGq\n1Wcv6kccAV26QM+e0KcP9O0L/fpB584RAvUxM9z9gOonqQyBVaucbt2SbomIZNmuXfDhh/DOO3Fh\nX7o0LuyrV8PatbB+fVzYKyriwg7QrFlUMUpKoG1b6NAhLuA9ekDv3nFRHzAAjjwyLv75kksIpLIc\nVFGRdAtEpFht2wYLF8K8efHfJUuixl5eDh9/HBf2bdtg5844vkWLuLAfeigcfniUYI455rN36zU/\nt0jlFXXvUtnkbduSboGINCXucWc+b16UZBYu/Gydff162Lw5auy7dtWWYkpKoF27uLD37QsTJ8JR\nR8Xd+qBB0LXr3sswxSCnEDCzdsC9QG9gKfBFd99Yx3FLgY3ALqDS3Ufv7bwKARGBuLiXl8Nbb8H8\n+bWjY2pKMjUdqDt2xPE1nadt2tTW2EeOjDv2fv3iwj5wYNTgJeT6JHAd8LS7/9zMvg/8oPq1Pe0C\nSt19fUNOqnKQSPGrqIB334U334y795rO1PLy6EitGR0DcXGv6UDt2DE6T4cPj7v2/v1h8GA4+mho\n2TLZ/6emKNcQOA8YX/39/wBl1B0Cxn5MTNu6NcdWiUhiau7e33ijtu6+dGmUZmru3isqYuij2Wfv\n3Lt2hXHj4oI+aBAce2x0ouri3nhyDYEj3L0cwN0/NLMj6jnOgafMrAq4zd1/u7eTfvJJjq0SkUZR\nVRUdqa+9VnuBX7EiRtGsWxd195rSTIsWcfd++OHQqVOUZCZMiHr74MEwbFj8TpK1zxAws6eAzru/\nRFzU/7mOw+sbb/o5d19tZp2IMJjv7rPqe0+FgEjhuUet/bXXokTz7ru19feaUTM1QyFbtYq79w4d\n4u79hBNqR8kMHx5378XeoVos9hkC7n5afb8zs3Iz6+zu5WbWBVhTzzlWV//3IzN7EBgN1BsC9903\nlaVL4/vS0lJKS0v31UwR2Yf16+MC/8YbtTNVV62KEs2mTbX195YtYzhkTXlmzJiouw8dCscdF0Mh\nm2nVsUSVlZVRVlaWl3PlNFnMzG4C1rn7TdUdw+3c/bo9jmkNNHP3zWZWAswAfuzuM+o5p994o3Nd\nXT0LIlKnqqoozcyZExf5996r7WTdsCFq8O5xd14zLLJLl9oJTEOGwPHHRy2+KY51z7okJ4vdBNxn\nZn8HLAO+WN2grsBv3f0copT0oJl59fv9X30BUGPz5hxbJVJkduyIGvycOfD223GRX748FhDbtKm2\nDn/wwVFn79w57tgnTIj6+4gRUao56KBk/z8kfXIKAXdfB5xax+urgXOqv18CDN+f827ZkkurRJqe\nbduiVPPqq3GxX7QoZrF+9FHcFFVWxkiamo7Wrl3j4n7MMTGCZtSouOhr9V3ZX6l88FMISLGpqIgy\nzezZ0elac5FfuzYu8lVVUWcvKYlafLducWEfNCg6WkeNihE2IvmWyhDQPAFpatxjJM2LL8Yd/YIF\n8fOHH0a5prKy9iLfoUMsKHbSSXE3P3JklGvatEn6/0KyKJUhoGUjJI02b447+TlzYhmD99+P0TXr\n1sXfWfeoubdrF+WagQPhggviIj9mjO7kJZ0UAiLV3GPY5AsvRG1+/vy4m1+zJsbIV1XF6JrDDosL\neu/eMHp0lGvGjYtJUBo6KU2NQkAyxT2GUj7/fJRtai70H30UZUj3GGHTvn2UbEaNipmto0fHVz7X\ngBdJA4WAFJ1du2IjkN3v6Jcvj07Ymv6m1q1jIbJevWL54BEj4HOfixq9ZrpKlqQyBGpmLorUp+ZC\n//zz8Mor0RG7YkUsb7BtW+1wypqyzVlnxYX+xBNj9qvKNiIhlSGgJwGpsXo1PPccvPRSdMYuWVJb\nujGL0TadOsVaNZMnx4V+/PhYrExj5kX2LZUhsGFD0i2QQqqoiKGVf/4zvP56zIZdvTqWHN61K0bc\ndOwYd/Snnx61+QkTYi15XehFcpPKjeYPOsi1sUyRqemQfeaZGGL5zjtRp//441jyoHnz2Ji7e/co\n1xx/PJx8clzwW7VKuvUi6ZbL2kGpDIEWLZx16zR5pinasSNKN88+G52yNXf1mzdHEJSUxJZ/Rx0V\no27GjYu7+o4dk265SNOV5AJyjeKww2IXogEDkm6J1GfLlrjQz5wJc+fGMghr1kR/TrNmsb5Nz55x\nRz96NJxySixFrA5ZkXRJZQi0aRMzMRUCyVuzBmbMiOGWb74ZHbM1JZyWLWM8fZ8+cNppMcTy1FNj\nfL2INA2pDIGSkngSkMJZuTIu9rNmxcV+2bLYhKSqKjpma0o4X/xirHkzcaK2BhQpBqkMgUMOUQg0\nlnXrYPp0KCuLMs6SJfFaVVV87l26RMfspElQWhoXfK1BL1K8UhsCq1Yl3YqmbevWGInz7LOxPMKi\nRTG+vrIyLuqdO8fF/soro4Rz4okahSOSRakMgVat9CTQUDt3Rr3+qafg5ZdjNE55eYy9b9kyRt0c\nfTR84QsxCmfixNg/VkQEUhwCehL4a4sXw7RpsVTC229HHX/Llhhj365dzJo9/fQYXz9pkoZdisi+\npTIEmjfP9pPAjh1Rxpk+PSZWLVwYI3KqquIuvlevWL74mmtiTZw+fZJusYg0VakNgdWrY8mAYh9X\nvnw5PPJI3N2/+Wbt3X2LFnEn378/XHpp3NmffHKUeERE8iWVIbBzZ8wV+Pjj4tmNqaoqRuQ8+WTs\nTrVwYXTUVlXFkNiePWPD8GuugXPOibt9EZHGlsoQ2L49NtpetapphsD27fDEE/H18stRy9+0KZ5w\nOnaEfv3g4ovhjDOis1Z39yKSlFSGwMaNsZDYBx9E7TvNtm6Fxx+PC/6cOTHufvPmuLB37QpDhsQF\nf/LkuPiLiKRJKkPgnXdiZury5Um35LM2boTHHosO21deiVm1W7bEaKZu3WJBtCuugPPPVzlHRJqG\nVK4iOnCg86UvwRtvwP33J9OO9etjOOaTT8Zkq2XLYnG0Vq1ibZxhw2KS1eTJ8dQiIpKUoltFdOTI\nWJfm6afjTrukpHHfb/PmuNhPm1Zb0qm54PfsCccdB9/6Fpx3XiyrICJSLFIbAvPnw5gxUWu/4IL8\nnbuiImbXTpsW694vXlxbw+/ePfogvvlNOPdcrYYpIsUvleWgWbOcb30LvvrV2F/27rsP7FyVlTH+\n/qGHYvvChQtjlE6LFnFHX1PSOf98TbgSkaYrsZ3FzOwCYCpwDDDK3V+r57hJwM1AM+B2d79pL+f0\nzZudTp3g3XdjI5KFC/c9VLSqKi70Dz4Ya+m8917sVdysWSyDPHRoDMecMiVG6WhvWhEpFrmEQK7z\ncd8Czgdm1neAmTUDbgHOAAYDXzKzgXs7aUkJ9O0bG5pcey2ccEIEQo3KytiU/Hvfi9917BjlnPHj\n4Q9/iFVIr702lkqurIyhptOnw3XXxQzcphIAZWVlSTchFfQ51NJnUUufRX7kFALu/q67LwT2dlkd\nDSx092XuXgncA5y3r3OPHg033QSXXRbDLo87LtbNOeSQ6LAtLYU774yL/9e/HpOyKitjBc3nnoMb\nbohyT1O54NdFf8mDPoda+ixq6bPIj0J0DHcHVuz280oiGPbqZz+DX/wCRoyIC/+IEbHd5LBhEQBD\nhhT/ukIiIo1tnyFgZk8BnXd/CXDgh+7+aGM1rGPHCIIbb2zad/MiImmWl9FBZvYc8N26OobNbCww\n1d0nVf98HeD1dQ6bWbqGK4mINAFpmCxWXwPmAH3NrDewGrgI+FJ9JznQ/xEREdl/OVXVzWyyma0A\nxgLTzOyJ6te7mtk0AHevAq4GZgDzgHvcfX5uzRYRkXxI3WQxEREpnETG15jZJDNbYGbvmdn36znm\nV2a20Mxv3+w0AAADLUlEQVTmmlnKF5Q+cPv6LMzsYjN7o/prlpkNTaKdhdCQvxfVx40ys0ozm1LI\n9hVSA/+NlJrZ62b2dnW/XFFqwL+Rw8zskeprxVtm9pUEmlkQZna7mZWb2Zt7OWb/rp3uXtAvIngW\nAb2BlsBcYOAex5wJPFb9/RjgpUK3M0WfxVigbfX3k7L8Wex23DPANGBK0u1O8O9FW6K82r36545J\ntzvBz+IHwI01nwPwMdAi6bY30udxIjAceLOe3+/3tTOJJ4GGTB47D7gTwN1nA23NrDPFZ5+fhbu/\n5O4bq398iZh3UYwaOqnw74E/AWsK2bgCa8hncTFwv7uvAnD3tQVuY6E05LNwoE31922Aj919ZwHb\nWDDuPgtYv5dD9vvamUQI1DV5bM8L257HrKrjmGLQkM9id1cCTzRqi5Kzz8/CzLoBk939v9j7LPWm\nriF/L/oD7c3sOTObY2ZfLljrCqshn8UtwCAz+wB4A7i2QG1Lo/2+dqZyKWn5a2Y2AbiceBzMqpuB\n3WvCxRwE+9ICOB44BSgBXjSzF919UbLNSsQZwOvufoqZHQ08ZWbD3H1z0g1rCpIIgVXA7psv9qh+\nbc9jeu7jmGLQkM8CMxsG3AZMcve9PQo2ZQ35LEYC95iZEbXfM82s0t0fKVAbC6Uhn8VKYK27VwAV\nZvY8cCxRPy8mDfksLgduBHD3981sCTAQeKUgLUyX/b52JlEO+nTymJm1IiaP7fmP+BHgUvh0xvEG\ndy8vbDMLYp+fhZn1Au4Hvuzu7yfQxkLZ52fh7kdVfx1J9At8owgDABr2b+Rh4EQza25mrYlOwGKc\nf9OQz2IZcCpAdf27P7C4oK0sLKP+p+D9vnYW/EnA3avMrGbyWM3+AvPN7Kr4td/m7o+b2VlmtgjY\nQiR90WnIZwFcD7QH/rP6DrjS3fe5AF9T08DP4jN/pOCNLJAG/htZYGbTgTeBKuA2d38nwWY3igb+\nvfgJ8N+7DZv8R3dfl1CTG5WZ3QWUAh3MbDlwA9CKHK6dmiwmIpJhWoxZRCTDFAIiIhmmEBARyTCF\ngIhIhikEREQyTCEgIpJhCgERkQxTCIiIZNj/ByBqpCt6Uy+yAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def backwardEuler(h=1/100, tstart = 0.0, tend = 1.0):\n", " n = int((tend - tstart) / h)\n", " ts = np.linspace(tstart, tend, n+1)\n", " xs = np.zeros(n+1); ys = np.zeros(n+1);\n", " xs[0] = 1.0; ys[0] = 0.0; # Set the initial condition here\n", " Delta = 1 + 401*h + 400*h**2\n", " for k in range(n):\n", " xs[k+1] = 1/Delta*((1+799*h)*xs[k] + 798*h*ys[k])\n", " ys[k+1] = 1/Delta*(-399*h*xs[k] + (1-398*h)*ys[k])\n", " return ts, xs, ys\n", "\n", "# With the backward Euler method, there is no sign of instability,\n", "# even for large step sizes.\n", "ts,xs,ys = backwardEuler(h=1/200)\n", "plt.plot(ts, xs, 'r-', ts, ys, 'b-');\n", "plt.plot(texact, xexact, 'r-', texact, yexact, 'b-');" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Implicit methods can have better stability properties than explicit methods;\n", "# but they are less straightforward to implement for non-linear ODEs, as a method\n", "# is needed to solve the implicit equations" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "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.5.1" } }, "nbformat": 4, "nbformat_minor": 0 }