{ "nbformat": 4, "nbformat_minor": 0, "metadata": { "colab": { "name": "SolvingODEs.ipynb", "provenance": [], "collapsed_sections": [] }, "language_info": { "name": "python" }, "kernelspec": { "name": "python3", "display_name": "Python 3" } }, "cells": [ { "cell_type": "markdown", "source": [ "#Solving Ordinary Differential Equations (ODEs) using Python" ], "metadata": { "id": "TzR4FNN5XfxL" } }, { "cell_type": "markdown", "source": [ "Simple differential equations can be solved numerically using the Euler-Cromer method, but more complicated differential equations may require a\n", "more sophisticated method. The \"scipy\" library for Python contains numerous functions for scientific computing and data analysis. It includes the function\n", "\"odeint\" for numerically solving sets of first-order, ordinary differential equations (ODEs) using a sophisticated algorithm. Any set of differential equations can be written in the required form. \n", "The example below calculates the solution to the second-order differential equation, \n", "\n", "$$\n", "\\frac{d^2x}{dt^2} = ax + b\\frac{dx}{dt}.\n", "$$ \n", "It can be rewritten as the two first-order differential equations, \n", "\n", "$$\n", "\\frac{dx}{dt} = \\dot{x}\n", "$$\n", "and\n", "$$\n", "\\frac{d\\dot{x}}{dt} = ax + b\\dot{x}.\n", "$$ \n", "Notice that the first of these equations is really just a definition. In Python, the function $x$ and its derivative $\\dot{x}$ will be elements of an array. The function $x$ will be the first element $x[0]$ (remember that the lowest index of an array is zero, not one) and the derivative $\\dot{x}$ will be the second element $y[1]$. The index indicates how many derivatives are taken of the function. In this notation, the differential equiations are \n", "\n", "$$\n", "\\frac{dx[0]}{dt} = x[1]\n", "$$ \n", "and\n", "$$\n", "\\frac{dx[1]}{dt} = ax[0] + bx[1].\n", "$$ \n", "The \"odeint\" function requires a function (called \"deriv\" in the example below) that will return the first derivative of each of element in the array. In other words, the first element returned is $dx[0]/dt$ and the second element is $dx[1]/dt$, which are both functions of $x[0]$ and $x[1]$. You must also provide initial values for $x[0]$ and $x[1]$, which are placed in\n", "the array $yinit$ in the example below. Finally, the values of the times at which solutions are desired are provided in the array $time$. Note that \"odeint\" returns the values of both the function $x[0]=x$ and its derivative $x[1]=dx/dt$ at each time in an array $x$. The function is plotted versus time." ], "metadata": { "id": "rP1rs_5MXfxM" } }, { "cell_type": "code", "source": [ "from scipy.integrate import odeint\n", "import pylab as pl\n", "\n", "def deriv(x,t): # return derivatives of the array y\n", " a = -2.0\n", " b = -0.1\n", " return pl.array([ x[1], a*x[0]+b*x[1]])\n", "\n", "time = pl.linspace(0.0,100.0,1000)\n", "xinit = pl.array([0.0005,0.2]) # initial values\n", "x = odeint(deriv,xinit,time)\n", "\n", "pl.figure()\n", "pl.plot(time,x[:,0]) # x[:,0] is the first column of x\n", "pl.xlabel('t')\n", "pl.ylabel('x')\n", "pl.show()" ], "outputs": [ { "output_type": "display_data", "data": { "text/plain": [ "
" ], "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEGCAYAAABLgMOSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2deZhkZX3vP2/tVb1v07NvzAzDMMAAwwDKFgFBRfEaF0xUNBDjNWbRGxNu3CKJiUZvNHFJRAHRBDdQQUURWWSRbVhngdm37umZ6em9u/aq9/5xzqmu5ZxTpzeK6f59nmeeqXPqvF1vdc2cb/12pbVGEARBECaKr9YbEARBEE5MREAEQRCESSECIgiCIEwKERBBEARhUoiACIIgCJMiUOsNvJK0t7fr5cuX13obgiAIJxTPPPPMca11R/n5OSUgy5cvZ/PmzbXehiAIwgmFUuqA3XlxYQmCIAiTQgREEARBmBQiIIIgCMKkEAERBEEQJoUIiCAIgjApREAEQRCESSECIgiCIEwKEZAJcmQoyd0vHK71NgRBEGrOnCoknA7e9o3HODyU5LwVrcxrjNR6O4IgCDVDLJAJcngoCcCW7qEa70QQBKG2iIBMgOLpjUeGkzXciSAIQu0RAZkAQ4lM4fHxkXQNdyIIglB7REAmwEC8SEBGUzXciSAIQu0RAZkAw0UWyGDRY0EQhLmICMgEGE6Oi0Y8la3hTgRBEGqPCMgEGE4YotEUDTKWFgERBGFuIwIyASwLZH5jhLFUrsa7EQRBqC0iIBMgnjZEY15jWCwQQRDmPCIgEyCZMQSkrS7EmMRABEGY44iATIBEOodPQXMsRFxcWIIgzHFEQCZAIpMjGvRTHw4wls6WVKYLgiDMNURAJkAikyMa8hML+8lrSGbytd6SIAhCzRABmQDJdI6IaYEAEkgXBGFOIwIyASwXVixkCogE0gVBmMOIgEwAy4VVH/YDSC2IIAhzGhGQCZAwXVgFC0RcWIIgzGFEQCZA0nRhRYKGBZKSILogCHMYEZAJkMjkiIX8hAPGry2VFReWIAhzFxGQCWAF0cNBS0C8WSC5vOazP9/GtsMyBlcQhNmDCMgESKTzREJ+IgHTheXRAtnSPcStj+3nXd98Yia3JwiC8IoiAjIBkuUWiMcYyPbDwwCMStqvIAizCBEQj2itx11YBQvEm4D0FY2/tRoyCoIgnOjUVECUUlcqpXYopXYrpW6wef4ipdSzSqmsUurtZc9dq5TaZf65dqb3ms1rcnlNOOArBNG9ikHfWLrwuL/osSAIwolMzQREKeUHvg68AVgHvFspta7ssoPA+4Hby9a2Ap8BzgU2AZ9RSrXM5H7TprURKhIQrxZIvwiIIAizkFpaIJuA3VrrvVrrNPAD4OriC7TW+7XWLwLld+orgPu01v1a6wHgPuDKmdxssYAE/D78PuU5iF4sGn0iIIIgzBJqKSCLgENFx13muWldq5T6oFJqs1Jqc29v76Q2CpDOjQsIQDjg8xxE7xtLs7K9DoDBuAiIIAizg1kfRNda36S13qi13tjR0THpn1OwQPzGrywS9E/AhZViaVsMkEwsQRBmD7UUkG5gSdHxYvPcTK+dFKmsjQXiwYWltWZgLMOSFlNAkiIggiDMDmopIE8Dq5VSK5RSIeAa4G6Pa+8FXq+UajGD5683z80YlgUSLhGQ6hZIOpcnncvT2RhGKWkBLwjC7KFmAqK1zgIfwbjxvwT8SGu9TSl1o1LqLQBKqXOUUl3AO4BvKqW2mWv7gX/EEKGngRvNczNGZQzE7ykGkkgbVkpdOEBdKMCotIAXBGGWEKjli2ut7wHuKTv36aLHT2O4p+zW3gLcMqMbLGI8BmIUEYaDPpIeXFhjpoDEQn7qwn5GU5mZ26QgCMIryKwPok8XabsYiCcLxHBZxUIB6sMBGUIlCMKsQQTEI+mcceMvcWF5sUBS4xZIfTggWViCIMwaREA8UpnG6y2IHjddWNGQnzoREEEQZhEiIB6pTOP1e+qFFTddWHUFF5YIiCAIswMREI+Up/GGAj4yOV11XTw9eRfWdx/fzwVfeIDjRd18BUEQXi2IgHikPI036Fdkcl5cWGYQPRyYsAvrxp9vp2sgwcM7J9+CRRAEYaYQAfFIeQwk6Pd5FBDTAgn6iYX9heNq5PPj1s3uY6MT3a4gCMKMIwLikfI03qB/4kH0SMBPOpsvEQcn+sbSZM3rugYSk922IAjCjCEC4pFyATFiIN5cWH6fIhzwEQkaRYheChCPDicLj3uGREAEQXj1IQLikXQuj1IQ8CnAcGV5DaLHgn6UUkSD1iTD6sJzbMQQkPmNEQbjUr0uCMKrDxEQj6SzeUJ+H0oZAhL0+8iZY27dSKRzREKG5VGwQDyk/w4njGD7srYYAyIggiC8ChEB8Ugqmy+4rwCCAUNIqrmx0tl8IfU3agpJwoOADCUM0VjWFmMokUbr6taOIAjCK4kIiEfSuXEhgPFsrHQVAUnl8iXFh+DVAjEEZGlrjExOe87eEgRBeKUQAfGI5cKysEQhUyUTq3hdZAIxkOFkhmjQT0dDGIABGYUrCMKrDBEQj2RzeYLFLixTFKoF0tNFrq/oBGIgQ4kMTdEgzbEQgATSBUF41SEC4pFMThcysGBcQNJVLJBMLl8yRx28B9EbowGao0FgPCYiCILwakEExCOZXL4gGmC0MoHqMZBiC8QSEK9B9KZokJY6wwIRF5YgCK82REA8Ui4gVkC9ahZWzs6F5S0G0hgJUh82hkZ67eKrteZzv9zO76R/liAIM4wIiEeyeU3AX+nC8pLGWxlE9+DCSmZojAapCxkC4nWW+tP7B/jWI/u49panPF0vCIIwWURAPJLJ5Qn6KoPo1WIgxS6s8ARiIGOpHPXhALGwsSbu0QJ5sWuw8HhQ3F6CIMwgIiAeyeR0oXgQigSkWh1IkQUykSyssVSWWMhP0O8jFPAx5rEOZE/veOfevcfHPK0RBEGYDCIgHsnm8gR8NnUgVdJ4M0UxkKBf4VPVYyC5vCaVzRcq1+tCfs8xkMODyYJgHR6UJoyCIMwcIiAeyeR0SRDduklXLSQsEhClFJGgv2oWlvV8zBKQcICxtDcBOTKU5KxlzYAIiCAIM4sIiEeMLKwiF1ZgAmm8Zdlb1eImhSmGZgC9LhQg7jGIfnw0xcqOekJ+H/1jUjsiCMLMIQLiESMLy64S3XsQHYx+WKkq80AS6XILxO/JAtFaM5jI0BIL0lIXZGBMguiCIMwcIiAeSWdLLZCQhyysfF6TzesSAQl5sEDGUjYuLA8xkJFUllxe0xwN0RIL0S9ZWIIgzCAiIB7J5kvTeL0E0S33VrCsCWM1t1ciY4hF1HRhxUL+gqi4MWi6rJpjQVrrQmKBCIIwo4iAeCTrlMbr4o6yhCIcKI2BpKpkYcUrXFjeguiDCUMwWmIhWurEAhEEYWYRAfFIuiyN13JnuVogZXPUrcfVLJAKAQkFPM0DsSYXNseCtMYmZoHc/Og+Ntz4G3YdHfG8RhCEuU1NBUQpdaVSaodSardS6gab58NKqR+azz+plFpunl+ulEoopZ43//zXTO81m9OlMZBA9ULCgoD4J2aBjAfRTRdW2M+ohxiIVXnebFogg4lM1ZG7Fjc/spfBeIafPd/t6XpBEIRArV5YKeUHvg5cDnQBTyul7tZaby+67DpgQGu9Sil1DfAF4F3mc3u01hteqf1WdOP1Vc/CsrdA/FVbs5dbINGgn3Q2Tz6v8RW1lC/HmhnSEgvSHA2itTHZ0Oro68RYKkvPcBKArd3DrtcKgiBY1NIC2QTs1lrv1VqngR8AV5ddczVwm/n4DuBSpZTzHXSG0FpXpPH6fIqAT7lmVFnWSahsFK7XOhCrEt1qA5+qss4SpsZokIaI1YSxuuWy7/gYWhvNHncfG616vSAIAtRWQBYBh4qOu8xzttdorbPAENBmPrdCKfWcUup3SqkLZ3KjWdMNFPKXalfQ7/NkgZS0gQ/6qtaBFCyQ4LgFAtXniIymskSCPoJ+Hw0RYxDVSLK6gBw1rY+zl7VwdDhJ3qPbSxCEuc2JGkTvAZZqrc8EPgbcrpRqtLtQKfVBpdRmpdTm3t7JzciwRKLYAgHDsnALoqdsXFhhDxZIMpMj6FeF15uIgFjzQyZigRwxBeSMxc1k85rjY6mqawRBEGopIN3AkqLjxeY522uUUgGgCejTWqe01n0AWutngD3AGrsX0VrfpLXeqLXe2NHRMamNWiIR8FVaIG5BdEt4wuV1IFUEJJXNEw74C8dhj3NExooExPp7JFm9ncnR4RQ+BactajKOh0RABEGoTi0F5GlgtVJqhVIqBFwD3F12zd3AtebjtwMPaK21UqrDDMKjlFoJrAb2ztRGMzaxDDBcWq4xEDsLJOCrGstIZXOF4VNQZIFUSeUdTWapswRkAhbI0aEk7fVhFrVEgXGLRBAEwY2aZWFprbNKqY8A9wJ+4Bat9Tal1I3AZq313cDNwPeUUruBfgyRAbgIuFEplQHywIe01v0ztddswQIpFZBgwFsMZKKtTJKZUgvECqZXs0BGU+MCYrmwvMRAjgwnmd8UYX5TxDge8t7F91dbeljQHGXDkmbPawRBmB3UTEAAtNb3APeUnft00eMk8A6bdXcCd874Bk0yhZYkpS6sULUguo3l4qWZouHCGl8T8ThLfSydZV6DIQINYSOI7sUCOT6aorMxQltdGKWgd9RbAeL2w8P87/95FoDdn3tDRYxIEITZjfyP90DGpqeVdZzOeqhEL4uB5LUxoMqJVCZXIjpeg+hjqVzBAokEffh9ylMMZCiRoTkWxO9TNEaCnkfhPrjjWOHxi91DntYIgjB7EAHxgJXGGyhP463SlsQujddLBXsymy9YHVBsgXjPwlJK0RAJMOrBhTUUz9AcNYoNW2LBQkuUamw7PISVV7DziLRAEYS5hgiIB+yEAIwguttEQqdmioBrO5NUJlfmwjIeV03jTWapD48LT304wEgVF1Yml2cklaU5Zri8mmMhzxbI3t4xLlrTQcjvY5/MXxeEOYcIiAcsC6QiBuLRAikPooO7BZLK5gkXWSBRDxZILq9JZMZdWGAISDULZNisXm+KGgJiWCDeBOTwYIKlrTGWtcXY0ysCIghzDREQD2RdYiBusQynIDq4WyDJTI6IbRDdWUCsdu/1RQISC1Wfvz6YGO/gC0Yr+EEPLqyRZIbhZJaFzVFWdtSx77i0QBGEuYYIiAcsIShP4w34fKS9tHO3jYG4zBEps0AihToQZ9GJp0o7+FqPq00ytMTCskCaPQpIz5BRK7KwOcqCpihHhrzXjrzUM8wnf7aFvlEpWBSEExkREA9YdSCVLixVtQ7Ep0pboFhi4lZMWJ7G6/cpQgGfqzVhPWd18LUeV5sjMpQYbwEPhgtrNJWtWqvSPWjUiixqjtLZGGEsnfM0dhfgUz/byn8/cZCvPrDb0/WCILw6EQHxQDY/eRdW+RqrLYmbgCQzpZXoAJGAz9WFZVWpF2dveROQcgskWHLeCcvimN8UYV5DGIDekeoWxVAiwzMHBwC4/+WjVa8XBOHViwiIB6xaj4o0Xr97M8V0Nl/R/iRcGIVbzQLxl5yLhvzuAmLOUS+2QKIeJhlalepW5fp4F193Aek3px221YWY12gIyDEPArLr6Ahaw6YVrRzqTzDsoU5FEIRXJyIgHnCzQFyzsHKlrigoioG4CkiuYl006B4Qt+Ij0SIBqQv5SVSZpW5Vqlc2YXRf1zeapi7kJxL0F6rfj41Uj4PsMEfmvum0BQAyf0QQTmBEQDzgVIke8lePgYTKXVgB9+FQubwmk9MlrigwXFPuFojxXLTchZXJobWzlTSazBLwqYJgee2h1T+WorXeiJtYLqxjw9UtkJ1HRqgL+blojdEZeSICMhhPi8UiCK8iREA84NTOPeD3uRcS2riwqlkgVp+scgskEvSTcEn9tQSkJAYSDqC1ew+tMbMBozXo0asLq28sTVudIRxWG5Q+D3NEDvTHWd5ex9LWGKGAjz0eBeTIUJKLv/gQb/z3R6p2JRYE4ZVBBMQDVhZWuRhMKgZiVaI7NFS06kMqBcQ9iJ40b6rRsiwsGK8RsWM0lSupHSlYIFUyqvpG07SZs9aVUjRHvbVAOTKUZEFTBL9Psbg5SteAt86/dz7bxVAiQ9dAgvtekuC7ILwaEAHxQGEioa+8G68ik887uojssrCqWSDJbKUlAYZryosLK2ZTwe72jd2wQMbXNHochds/lqbVFBAwrBAvLVCOjaSY12jETOY3RTjssXX8E3v7OLmzgfpwgKf3zVjnfkEQJoAIiAecRtoG/T60NuIWdri6sBxiJwULpCyNNxryuwpB3MYCsdqauFkgY+lsSfsTS0zcXFhaa/rjpQLSEgsxMOZugaSyOfrH0sw3BcRrAaLWmm2HhzljSRNnLGniWTMNWBCE2iIC4gHLTVUeEA+aYuDkxkrn7ILo7s0Uk4UYSFkQPVAlCytTGTuxxMQtlbe4gy8YIhkL+V0tkFQ2Tzqbp8msGQGjELFaDy0ryG4JyMLmCEeHk661NGAMvOofS3PqwiY2LGnm5SMjVQsdLR7Z1cu1tzwloiMIM4AIiAeyBQuksg4EnK2JqVgg5YWE4aDfNRiezOSIBv2FYDhAndnWxM1yGU1mC9dZNEQCrhaI1YDRcneBUcFerQXKUXNUbmfTuAWS19XrR3aYreJPWdDIqnn15PKag/1x1zVgWC433LmF3+3s5dN3ba16vSAIE0MExAOZvH0WVsgUFKdU3nTWpg6kSisT63xFIWHQT6pKJXqx+wqKguguAfGxVKkLC4xMLDcLxEqlbYwWCUhddQvEmrXeaRYezm8Kl5x34pApFsvbYqxorwfw1D7+pZ4RugcTrOmsZ2v3cOHnCIIwPYiAeCCTyxP0q5Jv9zAeE3EUkFylBaKU0dfKKQsraeOKAjMLy2UUbsK0QIqxBMXN9WW4sErXGRaIs4AMJYznGiPjwtMcC5LK5l2tnUL7E9OF1V5vCEhflRG6B/vjhAM+OhrCrGirA/DU/ffhXb0A/OPV6wF4aGdv1TUWB/rGONAnLeoFwQ0REA9kc/mKTrww7sLKOIy1tSskBKOdiXMdiOXCqiwkzOS0Y8DezgKxXFNjKfubutaasXSO+oiNBeJitdhaIGYzRjcr5NhIinDAV+i7NS4g7i6sQ/0JlrTGUErRFAvSWhdi3/Hq1sSLXYMsbY2xaUUrLbEgW7u8jd3dfniYy//tYS7/8sO8fGTY0xpBmIuIgHggk9MVnXhhvDtvJu/swipP4wUjw8rZheVsgYDzTBA3CyTukIWVyubJ5XWlCys8mRhIdQE5PpqivT5csOSsLK7jVQTkYH+cpa2xwvGi5iiHB6un/77cM8IpCxpQSrF+URNbD3sTkO/8fh/pnJEo8J3H9ntaIwhzEREQD2Rs6jlgPJ7h5MLK2LiwrHWOdSAZ+xhIYSaIk4CkKwUkViULq7wPlkU1F9aw+VxjdHxdi5mR5RZIHyirHYkE/TSEAxyv4sLqGUqwwAy8AyxoitBTpX4kkc6xr2+MtfMbAVi/qImdR0dcW8+AkZL925eOcfWGhbx1w0Lu3Xak6hqL46Mpbn1snydxE4TZgAiIB7I5bSsgnlxYdgIScHNhWYWE5e3c3acSJjI5ImUurKDfR8jvcxQQK7g+LVlYddUtkP54pnCdRXtD2NUCSWfzDMQzhYaNYAyxOjzoHnjfe3wUrWFNZwMAqzrqyeR01UD69sPD9I+led3aeVy2rpOBeIbth6u7sfJ5zZ997xk++/PtvOfmJ6umJgvCbEAExAOZXL4ihRfG60AcU3IdLJBwwO+hlUmZBRKyBMTJcskRDVa+VjTkd3RhWVaGXRZWMpN3/OY9nMwQCvhK4jTWHJGBMWcBGRhL01pUOwJGO3i3IHqvKS5Wy3gw6kdGU1nXxopWixTL9bWiwwi+768SGH++axCAs5a2cPayFgCeOVC9huTxvX08c2CAi9d0sLd3jF9vO1J1jcXR4aRkiAknJCIgHsjkHSwQn3Mar9baSOO1c325WCCFQkKbgVLgboHEyiwJMFq6V7NAyl1Y1vGogxtrOJEtsT5gfCCV2yCqgbF0pQVS726BHDNTfK2Ov2BYIAA9LlaIJSCLWoxrreytvb3uAvLCoUHa6kIsbjFG9S5qjnoqQrzr+W4awgH+6z1ns6Apwt3PH666BuD5Q4Nc9K8PcvEXH+TBl495WiMIrxZEQDyQyebtg+gB5xhIxqEBIxgBcscgeqZyjjqMx0CcLJd4OleRuQXuLVCsFid1ZWm8VlbWqEMm1nAyUxL/AMNiigb9jgKSzuYZSWVpjZUKSFt9iD4Xq8UqMix2YS1oMkTBLdbQPZAgFvIXYjMtdSGaY8Gq9SMv9Qxz6qKmQqB//aJGtve4u7C01jy88zgXndxBNOTn0lPm8ciu4669yyz++Z6XaIgEWd5Wxyd+usWz66t/LM33nzrIwT6xXITaIQLigWzePY03a9PKJO0wQwSqxUAMt5evrGjREgdHF5ZNEB0M95RTL6xRM7233AJprDITZDiRqbBAwGqoaC8gVqNFOwtkIJ52vHEWBKTMhQW4NmLsGoizuCVaUruzor3OVUDyec3e3jFWddQXzp08v5H9x8dcxeBAX5wjw0nOX9kGwEWrO0hkcmzpds/62nl0hKf29fOhi1dywxvWcngoyQMerJB4Oss1Nz3O//3JFq7++qNSryLUjKoCopRaZ3PukhnZzauUTE4XrI1iLKvELgZizQlxDKI73DCTmcpphDAeVHeyJhKZHNGQTQwk6MGFFSl3YbnPBBlOZgtt34tpigYZdLBA+k0Baa0QkBBajz9fTu9wEqUotI4Hwxrx+5SrBdI1kGCR6eqyqCYgPcNJEpkcJ82rK5xbO7+BvIZdR50LFzebMZJzV7QCeI6d/PLFHnwK3rJhIa9bO4+2uhC/eLHHdQ3ArY/tZ+fRUf7hzetIZ/N8/lcvV10DsLd3lCu/8jAb/+m33LddWuILU8eLBfIjpdTfKYOoUuqrwL/M9MZeTWRy+UK8oxi3NF5LIBxdWA6WhN08dCiyQGxcWJlcnmxe21ogMZcgeiELyyaNF5xdWCOJTEkRoUVTNMiQgwVizVBvqXBhGZbF8RF7ATk2kqKtLlzSCdnvU3Q2hDky5Bw76R5MsLglVnJueVsdPUNJR2vCGm51UpEFsna+kcXlVlC47fAQ0aCflea6tvowK9rrqgrII7t62bCkmXkNEQJ+H5ecPI/f7ex1dWPl85rvPr6fC1e38/7XruD9r13Or7cdobtK6nA6m+fPvvcMvSMp2upC/MX3n2VPr7dhXj946iBv/fpjfPqura5tcYS5hxcBORdYAvweeBo4DLx2Jjf1aiOb0/ZZWG4CkrWPZQCEAn7nzK1sriKFF8Zne9i5sOymEVrEwoGqdSDlabz11VxYSWcXllMMxGr1XmmBmNXoDtMMj42kSgLoFvObIhwZtr9pDiczDCUyLG4ptUCWtRmC4pTxtNe8oa7sqCtaU0ck6ONls6GjHdsOD3PKggb8RV8yzlrawrMHBhxnxSRNF9emFW2Fc5eeMo+hRIbnDg06vtYLXYMcHU7x9rMXA3DNOUvRGu58pstxDcA9W3rYdWyUf3nbaXzv+k34leLL9+10XQPw0+e6uOEnWxhOZvjvJw7wF99/znVEssXW7iHe+c3Hed2XHuLWx/Z5WgOGQB7si4tQnSB4EZAMkACiQATYp7WeliR3pdSVSqkdSqndSqkbbJ4PK6V+aD7/pFJqedFz/9c8v0MpdcV07MeJTN6+kLAQRLepA0m5ubDcWplkKhswwnhWlt23Z8utZZeFFQv6iTu0MhlNZokG/SU3PjAq0cF5KuFwMlsRRAdojoYYTNhbEv2FGEhZGm+9ezX6sZFkSfzDYkFTlB6HWSLdZRlYFlZK7wGHwPOe3jEaIgE66sdfz+9TrJ7XwM6j9gKitealw8OsW9hYcv6sZc30jaUduwY/d3CQTE6zaUVL4dwFq9vx+xS/2+Hcs+u3Lx3F71NcsmYeAEtaY7zmpDbufLbL9Sb9o82HWNIa5bJTOpnXEOHa1yznl1t6XGfSp7N5vnTvTs5Y0sxv/voiPnXVOh54+Rh3v+CeYband5R3f+sJ9h8fo7UuxGd/vp3//N0e1zVgCPibvvooF33xQc76x/v49iN7q64BI7vvWw/v5Z/veYlHdx33LFZgZA16SXYQ7PEiIE9jCMg5wIXAu5VSP57qCyul/MDXgTcA68yfWx5vuQ4Y0FqvAr4MfMFcuw64BjgVuBL4hvnzZgSnSnTLrWVnTaQLXXWdWpk41IFkc+4uLBcBsYuB1IUDzi6sdGUnXhifi26XxpvM5Ehn8xMOog84uLDaq7mwhu0tkM7GCEeGkrY3C0tAyl1Yy8xU3gMON/U9vaOc1FFf0TRzdWe9o4B0DSQYSWVZt6Cp5PyZSwxheN7BmnhqXz9KwdnLWgvnGiNBNixp5tHdx23XAPx2+zE2LW8tmcVy9YaFHOiLs82h4PFgX5zf7+njnWcvKSRnXHfBCoJ+H9/5/T7H1/rpc110Dyb468tWE/D7uPb85Zy2qIl/uedlx5tuPq+54c4X8fsUP/nwa/jxh87nzWcs5Iv37uCZA86TJIcSGd53y1McG07ymTev48LV7fzTL1/iaw/sclwD8NzBAa74ysN87p6X+M5j+3nPzU/y8TterCoKD+/s5cqvPMwZn/0Np3/2N/zNj1+o2lJHa80ju3r5uzte5PrbNvOle3d4qt/RWvPCoUF++PRB7nq+29MQNYu+0RR7ekddm5TWksq7RyXXaa03m497gKuVUu+dhtfeBOzWWu8FUEr9ALga2F50zdXAP5iP7wC+poz/3VcDP9Bap4B9Sqnd5s97fBr2VUE2pytauUMVF5ZLDCTkd+uFlbd1YVmV6HbrLBeWXQwkGnIeRGXMQ7cTKx9+n7INots1UrRojBodeZOZypTi/rE0DZFAhRA3RgKE/D6O27iwcnnN8dFUSQqvxYKmCPF0jpFUZU1K14Dxn7rchdUSC9IQDjhmLe3pHeWCVR0V50/ubOAnzwl1SJ0AACAASURBVHYzFM+U3LjBiH8AnFpmgazprCca9PPcwUGu3rCo4mc+vb+ftfMbC/UzFq9d1c7XHtjFUCJT8dzBvjg7jo7wyTedUnL+ilPn84mfbuXnLx5m/aJSIQP48TOHUArevnFx4VxbfZi3nLGQnzzbzcevWFvxWtlcnq8/uIfTFzdxyRrjd+LzKT75plN4101PcNvv9/NnF59k+1pP7x/gX99+ekHAP/+203j2wAA33LmFX/zlBRVfkLTW/O0dL3BkKMmPP3Q+Zy5t4X3nL+dvfvwCX/rNThY2R3nbWYsrXmvHkRGuveUpmmMhfvEXF7C6s56vP7Cb/3hgN90DCW5639mFL0PFr3Xzo/v43D0vsaKtjo9fcTJHhpL88OlDPPjyMb78rg1ctKby30D/WJpP/mwL92w5QlM0SGdjmAd3HOObD+/hf198Eh/+g1UV/+a11vxuZy+f/9XLJS5QpYxMvT+7aCXnn9RW8YUllc1x77aj/OCpg/x+T1/h/DnLW3jrmYt48xkLbb/A7Ts+xq+29nDvtqPsPTYKyojnXbSmg9ev6+TUhY0VrzVVqgpIkXgUn/veNLz2IuBQ0XEXRrzF9hqtdVYpNQS0meefKFtb+b90mmiIBAqV1sVYLizbNN6scxpv2K2QMGNvgQT9Cp+yz8Jyi4HUhYwuvnZtVexmgYDRcr4+HLANog/btHK3sH5HQ4mMrYCUxz+s12qrD9laIH1jKfIaWxfWfLM31pGhpI2AJIgEfSWZW9ZrLWuP2bqwRlNZjg6nSuIfFmvMQPrOYyOcs7y15Lnth4fx+xQnm9dYBPw+TlvcZBvPyOTyPHtwgHecXXlTvGBVO/9x/y6e2NvHFafOL3nuty8ZmVOXr+ssOd8cC3HB6nZ++WIPN1y5tuQmkctr7nimi4tWdxTqZyze/5rl3PFMFz/efIjrL1xZ8txdzx/mYH+cT121seTnnbuyjYvXdPCNh/ZwzaalJcIzkszwr7/ewTnLW0reW104wD++9VT+5Dub+caDe/jo5WtKXuu23+/n3m1H+cQbT+HMpYbl5vcpvvCHp9MzlODv7nyRRc1Rzl05Hi/qGohz7S1PEQ35uf1Pzy2I1cdefzLL2+v4+B0v8u5vPcF3PrCpYOVmc3lu/MV2vvv4Ad542nz+7Z0bCv9O33v+Mj5y+7Nce+tT/PklqwpWF8Bvth3hEz/bymA8zd9eeTLXXbCCcMBPz1CCL/zqZf7jgd3c/cJhPnXVOv7g5HkoBVu6h/jivTt4ZNdxlrXF+MIfnsb5K9sZSWW4d9tRbn/yIH/07SdZv6iRP71wJRuWNHN4MMl924/y0+e6GIgbMbyPXraGpW1R9h2P86stPXzip1u58efbuXL9fM5d0UYo4GPXsREeermXHaaVfMaSZv7w7MXk8prtPcN89YFdfPWBXTz9icsKv4vpwosFckKjlPog8EGApUuXTupn/PhDr7E975rGm3MOoluFhFprm28feerq7G/qkaDf1jRPpt0skPGphOUCMuogIGCIpp0Ly80CaY4aN+zBeIbOxlKrYSCernBfWThVo1sjcO1cWFZzxZ6hZKHflUX3oJHCa/dta1lrnW1hoBVAL87AsrB+/o4jlQKy7fAwJ3XU2Yr3mUuaufWx/RVuyW2Hh4mnc5yzorVizYYlzcRCfh7bfbxCQO5/+Sir59UXXHHFvOm0BXz8jhd5sWuIM5Y0F84/squXnqEkn7qqIhuf9YuaOGd5C999/AAfeO2KQiwsl9d87cHdnLKgkctOmVex7uNXnMxVX32Ubz28l7+54uTC+a8/uIe+sTS3fuCcit/969Z2cvWGhXz9wd1cdkonpy02LKUXDg3yz/e8zKVr53H9hStK1oQCPv7rPWfztm/8nj/772f46Ydfy4r2Oo4MJXnvzU8RT2f50YfOr3BVvu2sxTTHgnz4f57lHf/1OP/vnWcQ8vv4l1+9xGO7+/jgRSu54cq1JbVWazobuOvPL+Af7t7G1x7czYM7jnHBqnaePzTIk/v6WTu/gds+sKkk1rWgKcpXrjmTPzx7MZ+5axvX3baZ5liQkN/HsZEUTdEgn75qHe85b1nJ/71TFzbx4UtO4qfPdfOtR/byVz94vvBc0K+4fF0n15yzlAtWtZfs8aOXrWZL9xA/3tzFXc93c5fZ7SDoV5y9rIXPvHkdV5w6v9CpwaJvNMVzBwenXTygtgLSjZHdZbHYPGd3TZdSKgA0AX0e1wKgtb4JuAlg48aN3qNrHgj6PGRhOdSBGOs0oUCZgGTyBXdVOdGg3zaN17JAnFqZgBHvKHe/jKWyheFO5dSHA4Wuu8XYNVK0sL6NDtrUdPSPpW2FAIxaELuxtsdGDF9xh40LyxKoIzbFhF0DlSm8FkvbYvxm+xFyeV2SPLCnICCVN+eFTRHqwwF22cRBtvcMF+o/yjlzaTPffDjP9sPDhW/WAE/vM2IBm5ZXrgsFfJy7opVHd5XGQYYSGZ7c28+fXrSyYg3A69fN5+/9W/jFi4dLBORHmw/RWhfislM6bddd+5rlfOT253hoxzEuNa/5+QuH2Xd8jP/847NsRXj9oiauOn0BNz+6j/e9ZhnzGiLsOjrCLY/u421nLeL0xc0VawBufMt6ntjbx4dvf4bbrz+PVDbP9d/dTEdDmC+94wzb12qOhbj1A+fw1q8/xjv+63EuXzePe7cdJZnJ8d0/2VTotlzO69Z28t/Xnct1t23mbd/4PWD8m/7CH57Gu86x/yIZDfn5wttP5zWr2vjPh/Zw62P7Wdwa5ZNvOoVrX7Pc1psAcOHqDn711xfy661HeGJvH5mcZsOSZt58xsIK16BFJOjn3ZuW8q6NS3hqfz+H+uO01oXYtKK1wu1moZTi9MXNnL64mc++5VQODyXI56GjIVwxC6iYtvowl62z//ynSi0F5GlgtVJqBcbN/xrgj8quuRu4FiO28XbgAa21VkrdDdyulPo3YCGwGnjqFdu5ic+nCPiUrYC4ZmEVNWEsfz6ZzVX0wbKIOMxFj7sE0aMuLd3dLJDGSJDRlF0MxBCVJrssrCIXVjkDY2nH/+zt9WFbq8DNAhkXkErh6RqIF77hlrOsNUYmpzk8aAypstjbO4bfp1jaVik8SilWd9YXXAQW/WNpeoaSFRlYFhvMQPpzBwdLBOSp/f0sb4sxz0G8X7uqnQd3vFSwpAAe2nGMbF7bWgQATbEgF67u4Jcv9vD3bzwFpRTHhg2XyPvOX2777xCM+Mn8xgi3PrafS0/pJJnJ8aXf7GDt/IYKC6iY//P6k7l32xE+9sMXuOENa/noD5+nMRrghjesdVzTFAvyzfdu5L03P8nr/t9D5PKa5liI71x/TkWHgmKWtdXx/Q+ex40/384vXuxhw5JmPvPmdaya1+C4BmDj8lYe/vgf8OCOY+TymktO7ijUHblx9YZFtnErN8IB/6TW+XyK81a2cV6Re87rOqcvSa8kNRMQM6bxEeBewA/corXeppS6Edistb4buBn4nhkk78cQGczrfoQRcM8Cf661rkmaQsCvCn2virHcWnZZWJZbK53NQ9m/Z6c0XjCyt2yzsFxjIMZHbJeJ5RQDAaMW5KjNrHJPFoiNgPTH07TW2X+zam8I0zeaJp/XJSa7ZZV02AhIKOCjvT5cUQsylsoW/Md2WAJxoC9eISBLWqK28ScwAunl1dtbu60Aur1YzW+KML8xUpKJlctrntzbx5XrnW/OF6xuB+Cx3cd550bD0L5nSw/zGsKF7C47rjp9AQ+8fIyn9vVz7so2/vuJA2Tzmveet8xxTdDv4/2vXc7nf/UyP3r6ENt7hukaSPA/159b0U6nmBXtdXzuf53GDXe+yFVffZRo0M/N799om/BQzIYlzfz8Ixdw+1MHCfoV7zlvWUVsxo618xu5/U/Pq3pdOU2xIG89c8bCo3OemsZAtNb3APeUnft00eMk8A6HtZ8DPjejG/RA0KGmY7yQsPKGFHZpjGgUEtrfxCIBewsk6ZKF5TZUym4eukV9OMAeuyC6WwzEskDKUnkT6RzJTN7xW2Z7fZhsXjOUKJ0XcmwkSVM06Pj7MAZLlYqcVZHt9O1seSGVd4wLaC+c39M7Wqgkt2NNZwM/ePpQYaoiUOh1td5BQMC4YRYLyPbDwwwns7x2VbvjmpM7G+hsDHPf9qO8c+MSxlJZHtrRyzXnLHG9qV+5fj7/9MuX+Pf7d/HvHfXc9vgBLl07j+XtlW65Yv7ktSu4b/tR/vbOFwF43/nLXPdn8c6NS1i3oNHoKLy6w9Z6s2N5ex1//8ZTql8ovOqZ9UH0mSbk95G1GWlbyMIKVP6HL7FAyki6WCARJwuk4MKyr0SHSgskl9ckM/mJB9ETWUJ+n+0e68MB/D5VUUxY6IPlGEQ3zveNpUoFxKEGxGJ+U6QiD99K4S3vg1VY0xghFPCVdLHN5zX7jo9x4Wrnm6YVSN95ZIT2VcaetnYPsawtVhFbKubMpc38etuRgvD8fo8R2zjfxWWhlOLNpy/ktsf3c3w0xX3bj5LK5rnqjIWOa8CIgX308jV86mdbueALD6A1/N2Vzi4li1DAx/9cfy53Pd9NUzTE6yfgL1+/qMk2dViYG0g33ikS9PtsK9HdsrAKMZAyAdFaOxYSAo5ZWAUXls26OgcLpNDK3SbwDoYLy66VidXK3S7gqZQy+mGVubAKRYQOFohV+d1blsp7bCRVkc1VzPzGSgvkUL9hgSxptRcQn0+xpCVaksrbPZgglc27WyDzjeeK8/m3dA9VvXlaqacP7zSqy3+3s5dV8+od4x8W7z53Kbm85tN3beVrD+xm/aJGNi5zdl9ZvOfcpXzijaewaUUrt37gHFZ3uscJLCJBP+86ZylXrp/vauUIQjEiIFMkGLAPortlYVnf3suLArN5TV7bx03AFBC7LKy00T/L7j9+IYhe1s7EqZGiRUM4QDqXr3CzDScyjlkiAM3Rymr0AYdOvBaFhoplqbzHhu3bmFjMb4owlMiU1MYc7I8TCfpK2pGUs6ytrmQy4V6zQ+9KF1fPvIYIC5oihbqOvtEUXQMJTqsiIKcvamJeQ5hfbT3CsZEkT+zt440u8Q+Lkzrque6CFdyz5QhHh5P8w5tP9VQEppTiTy9ayfeuO9eTG0oQpoK4sKZI0G/fmt21Et1BQJIuwXAw03gdminaxT/AOYg+LiD264rbmYTrx68ZTmZtiwgtmmwaKjp14rVot+mHlc/rqhaIVQtyZDjJCvPmf6g/zpKWmOvNdmlrjCf29hXqcAo1IPOcLRCAs5YZDRIBnthrpuI6pPBa+HyK/3XWIr718F7S2TwaPAd1//6Np3D+SW3Mb4w6ZnoJQi0RC2SKBH0+9zTeCbiwrDVOabyOWVgOw6Rg3AIZK3NhOQ2TsrDOl7uxhh1auVs02VkgY+4WSEsshN+nSgSkP54mm9eOdSowXo3eU1QLcmigND3XjpUddcTTOY6YWWYv9QzTWheqqFwvZ9PyVroHE+w7Psaju49THw5wugf///UXrKQ+HOB3O41AuJurrBilFK9b2yniIbxqEQGZIoYLy76VScjvs/0mHA7YZ2FZAuJUSOhUB5LI5Ig4FBKFA0Zfq/IWKFVdWA4zQZymEVo028RA+uMZlMKxqMrnU7TWlbYzsVKIO11cWAvN9E+reaLWmkP98ULXXSesvlXbuo3ak63dw6wvGmPrhNVC5GfPdfOrrT1cvKajZE6JEx0NYe7+yAV87Y/O5LNvWV/1ekE4URABmSJBv70FknaYow7jMY5yC8SyLhwLCQMOrUxcXFhKKWJBf8VY28I0Qpc6EBhP27UwguguAhILVVSiD4ylaY4GK9rGF9NeHy6ZCVIoInSxQBa1RPH7VCEgPhjPMJrKOtaAWJyyoBGlYOvhIVLZHDuPjrDew7f8hc1RNq1o5d/v38VgPMO7zllSdY3F8vY6rjp9oWNBnyCciMi/5iniJCAZmypzi+JK9GKsKYUTTuPN5Ar1HnbEwv5KCyRdLYhe2dJda23bJbaYpmiQ4WSWXH7cKuuPp10rjcGIg/SOjguP1cbELQYS9PtY3BJlnxkQt9q0V7NAYqEAq+fV88yBAZ45MEA2r0sqxd34hzefyprOet5z3lLXtF9BmAtIEH2KhPw+23bpdt1vLQpZWJnyGIhlgTi7sLJ5TTaXL3GdxNM5R0sCjBumUwzEOYhe6cJKZHJkcrqqgIDh6rJEY2As7VgDYtFRH2Zv73hm1FHTAnHLpgKjMHC/mUVlzezwkrp60eoOvvv4AZa31RH0K84/yVsriXULG/nNRy/2dK0gzHbEApkiQb9DGu8kLJBkFQukMNa2zPXlFkQHoxo9MUkXVnEQ3YptuAmIXT+s/rHqFkhHY5jekRR503I5MpyktS5U1eWzor2OfcfHyOc1u46OEA74qlogAJet6ySdy/O9Jw5w3so2VwEWBMEeEZAp4tbKxC4DC5wr0QsWiGMQ3X6sbTKTc+3GGQv5GbOpA1HKvv0JjAtLsQVizQLxIiADRXGQgXh1C2RRc5R0Ll/IxOoaSDhWkxezbkEj8XSOfX1j7Dg6yqp59a6xFotzV7Ry+bpOGsKBivkUgiB4Q752TRGnGEgqmyfkIAROvbAKWViOabz2Y23d6kDAcGGVNzgcTWWpC9lXlBt78BPy+yZsgbTWGS6nPjOeobVmYCxT1QKxxKJrMMG8xghd/XHWLqjuiiqeK/H8wQHesH5B1TVgJBd8630bK9q6C4LgHbFApkjQpRuvowvLwQIpZGG5pPEa15Wui6edGzCCYYHEU5UuLKf4h4XRzmRceCwBabRp5W5RXhQ4nMySzuWr1lgsahlPyc3nNV0DCZZ4aFe9ep4xOva2xw8wnMzymlUTa4st4iEIk0cEZIoE/T6ytmm8OUIOabyFSYZOhYROWVgBexdWIl0lCysUsOmFlXPMwLIoH2vrxQJpL2tL0mtmU7m1JIFxC6R7MMGxkRTpXJ7FHmIZAb+Pi9d08MKhQZSC15wkmVGC8EohAjJFggEfaRsLxJg2aP/rVUoVxtoWM+7CqmaBjItBOpsnm9dVBMRv28qkWuC4oayhohcBiQT9NIQDHDddWFY9h91Mj9LXCtJaF2Jf71ihT5WXYDjAX1y6io6GMB++5KSqryMIwvQhMZApEnIpJHTrGRWyE5AqhYRWoLzYhTXeyt0ljTfsr7RAzBiIG/Xh0pbu1jApt2aKYAyI6jUtEGsoVLVBQwBr5zfw8tERth82KsRP8RADAWOg09OfuMzTtYIgTB9igUwRxzRelzoQMNxUFYWEVV1YlRZIPGPc4N0skLpQgFQ2X1LcN5qq7sJqiARLKtGHEhkaIoGqcYP2+hDHRywB8ebCAmPq3M4jI2ztHqK9PuxJdARBqB0iIFMk4GSB5JyzsMAIlFcUEmZyKGXfgBGK0niLsrcsy6KaC8u4dtyaGHOZRmjRXNZZd7hKFbpFR4NR0wGGCysc8NHgoc7ilAUNJDI5fvJctzQQFIQTABGQKWKk8Wq0Lo2DuNWBgOHCqigkzBrTCN1Sa8HBhVUljRcoqQVxm4du0VoXKrRiB6q2MbFY1BylezCB1pqe4STzmyKeZlkUz6+4SNqECMKrHhGQKWJlWpWn8hp1IM43zZDfR7q8DiTjPI0QxmMjxa1TrMcxl3hGfaEtybg1MepBQJpjQVLZfEGkBuJpx5bsxSxuiZHK5ukdTXnqjmuxsDnKxy5fw8VrOibUqFAQhNogQfQpEjStjGw+T6hIjzO56haIXRaWUxEhjFsgqUylC8utEr2h0FnXcGFlc3lS2XzVILpVPT4QTxMNRekbS1edtQEUuuF2DSQ40BfnzWd4K+4D+MtLV3u+VhCE2iIWyBSxBKR8LrqnILpNIaGbBRK1SeO1ely5xUCs+R1WFpXVWLFaIWGzKSCWG6tvNE1bXfVg+GKzAHDb4WGGEhmWtTqPihUE4cRFBGSKBB0aI7pVooMZA7GxQJwysMAQK79PlcRAvATRG8saI1ZrpGhhuasG4xmSmRyjqSxt9dVdWEtaoygFD718DIClbd5cWIIgnFiIgEyR8RjI+E09l9fk8pqQ3/mmbhdEN1xY7lZBJFA6EyTuIYhuDYCyUnItIXEbDAXQYjZG7I+nC1ZItZYkYMRjlrfVcb8pIOsWSEaVIMxGRECmSMBnurCKxMCyLKq5sMrTeA0XlvtHEgn6S9J4ExOIgVjCUehpVaUgsHieh9Ucsa3KfA6Lc5YbA5ra60NVJwQKgnBiIgIyRSwXlp2AOI20BQgF/LYWiFMVukUk6CeRtnNhObujokE/fp8qxECGPTRFBCOIHvApjg4nOW6Om/WShQXwvvOX01YX4q8uXe0phVcQhBMPycKaInZpvJYwuFkTIZs5IqlsjuYqbqVw0FdaSJjJEgr4XKvDlVI0FvW18tLTCsDnU3Q2RugZSnJ0yKwo99hrav2iJjZ/8jIRD0GYxYgFMkUKWVjFFkjOgwsr6KuYB5LMVLdAokF/SRpvtU68Fg2RYKE1uxULqebCAljYHKFnKEHXQAK/T7GgyXt7EREPQZjdiIBMEVsB8RADCfnt6kDc03jB6qxbGkSPVQm8g+GusupArMmCDS7NHi3mN0U5MpSkayDOgqZIySx2QRDmNnI3mCLBwnCoIheWJSAuWVh2dSCpjHshIRixjrF0qQXiFkC3aAiPWyBDiQz14YAnMVjQZLiwDvbHPQ14EgRh7lATAVFKtSql7lNK7TL/bnG47lrzml1KqWuLzj+klNqhlHre/DPvldt9KUGbNF7PWVjZfEkPrWqFhGAU/xVPF4yns64BdIuGSKBgeQwnM66t5otZ0hIllc3z7MFByaYSBKGEWlkgNwD3a61XA/ebxyUopVqBzwDnApuAz5QJzR9rrTeYf469Epu2wz4GYlgI7jEQay76+LpqhYRgtGYfKxKQRMabBdIYLbVAqtWAWKwtquE4eb63+RyCIMwNaiUgVwO3mY9vA95qc80VwH1a636t9QBwH3DlK7Q/z9gJSMpDGm95WxKttZnGW80CqXRheQuij8dAhuLeuuoCnLaoqfD4vJUTmzcuCMLsplZpvJ1a6x7z8RGg0+aaRcChouMu85zFrUqpHHAn8E+6vJ+6iVLqg8AHAZYuXTrVfVdgddwtTuO1HrtZE+XTBasNk7IoH08bT+dY1OLBAokEGU1lyeU1x0dTnOJx3kYk6Ofb79vIvuNjrC8SE0EQhBkTEKXUb4H5Nk99ovhAa62VUrY3fxf+WGvdrZRqwBCQ9wLftbtQa30TcBPAxo0bJ/o6VXHNwnIJokfKWrN7FZC6cIBMTheaNY56mG0OFHpY9Y+l6R1JcZHHinKAy9bZ6bsgCHOdGRMQrbXjkGql1FGl1AKtdY9SagFgF8PoBi4pOl4MPGT+7G7z7xGl1O0YMRJbAZlpJpvGa7mwrFYkVk1INRdW8XTBUCDEaDJLfbi6O8oqADzYH2cklaXDY0GgIAiCE7WKgdwNWFlV1wJ32VxzL/B6pVSLGTx/PXCvUiqglGoHUEoFgauAra/Anm0ppPGWVKKbYuAiIIXpgqZwWH2xIh6C6GC0ZM/nNaPpbGFglBvzGo0CwO2HhwCjR5UgCMJUqJWAfB64XCm1C7jMPEYptVEp9W0ArXU/8I/A0+afG81zYQwheRF4HsNS+dYr/xYMCmm8xdlUmeoWSEFAJmiBWFMEx1JZxtJZtMbTvHHLAtl2eBhALBBBEKZMTYLoWus+4FKb85uB64uObwFuKbtmDDh7pvfoFbcsLNcguuXCMmMgVjC9ahDdHAI1lsoymvJeUd5RJiDtE4iBCIIg2CGV6FNkfKRtZSW6mzVRmYVlCEm1eSCWCyuezjFqpuV6cWGFA35aYkG2dBsurAVNUhQoCMLUEAGZIpYLK11SEOghBhKwt0CqxUCsIPpoKluo6/CShQUU5pk3hAMSAxEEYcqIgEwRpRRBv6pwYSkFAZcW65GQ8au3Cgm9zPWA8RbsQ4nMhFxYAGvNSvJVnfXSKVcQhCkjAjINBP2+CgEJB3yuN+nySnSrOLBaW5Imc8zscCIz7sLykMYLcNXpC0v+FgRBmAoyUGoaCPhU6UCpbL5qU8RIWR2Il9G0APWhAD4Fg/FMwfLwaoFctKaDpz9xmWRgCYIwLYgFMg2EAr6ypog51xReMKyWgE8V6kCsWEi12R4+n6IpGmQokSl01/USRLcQ8RAEYboQC2QaCAf8pUH0TPWuulA63zzu0QIBIw4ymMhQH08T9CtPdSCCIAjTjdx5poFwwFcYYwve2rKDISAFCySdw6eq14EABQskGvTRVheWgLggCDVBBGQaCAV8JXPKUx5iIADRkK9QiR5P54iFAp7EoCkWYiiRIehThSaJgiAIrzQiINNAeBIxEDBqQazYRyKTq1pEaNEUDXKwbwyANqkoFwShRoiATAMVMRCPLqxoyF9I402ks54GQwE0R4MMxDNk85qV7XWT27QgCMIUEQGZBkIBX8mQp1Q272niXyQ4boHEPU4WBOhsDDOUyDCSzNBWJy4sQRBqg6TxTgPlQfT0BILoCbOFidfZ5gDzzT5WeS0uLEEQaocIyDRgBNEnHgOJBseD715nmwMsah5vhLisLTbB3QqCIEwPIiDTQEUQ3WMdSLTMhRUNTqynFcCqefUT3K0gCML0IDGQacA+iO4ljdc/3spkAi6slroQ5yxvYTSVY7UIiCAINUIEZBowWpmM14GkszlPFkhdKFDoqDuSzNA4gZYkP/zg+aRzeSkiFAShZoiATAOVdSDeXFgNkSDxdI5sLs9wIktDxFtXXTB6YkV83iwWQRCEmUBiINNAKOAruLC01p4FpDFq6Pfx0TTpXL5wLAiCcCIgd6xpIBzwk81rcnlNNl99nK2FZXF0D8YBaJyABSIIglBrRECmgXDQsDbS2fy4gHixQMyYR9dAAvA+10MQRO+5OQAACOFJREFUBOHVgNyxpoGQ3xCLVDZHNm8MlvJSB2JZIJaANHqoXhcEQXi1IAIyDVgWSCqbLwjIRGIg3YOmgIgFIgjCCYTcsaYBywJJZ/OF2ehe6kCsmEe3ZYFIDEQQhBMIEZBpwAqYF7uwvMVArCC6FQMRAREE4cRBBGQasMQimcmTm0AMxJplfqjfzMKSNF5BEE4gpA5kGrAGQSUzuUJBoRcXlt+naIkFSWXz1IcDxEIiIIIgnDiIgEwDVhfdRCZXGBAVCXr71XY2RgCY1yht2QVBOLEQAZkGoqYFEk/nCoOlvDZGXNBkCMh8U0gEQRBOFERApgFLLBLpHHGzu26dR3fU6k6jNftJHdJVVxCEE4uaCIhSqlUpdZ9Sapf5d4vDdb9WSg0qpX5Rdn6FUupJpdRupdQPlVI1netqubDiRQLidTjUG9bPZ9W8et6yYeGM7U8QBGEmqJUFcgNwv9Z6NXC/eWzHF4H32pz/AvBlrfUqYAC4bkZ26ZGYOQgqkZm4C+vMpS389mMXc87y1hnbnyAIwkxQKwG5GrjNfHwb8Fa7i7TW9wMjxeeUMQDjdcAd1da/Uoy7sLJFFohkVAmCMLuplYB0aq17zMdHgM4JrG0DBrXWWfO4C1jkdLFS6oNKqc1Kqc29vb2T220Vgn6F36eIp3Mk0sYwKb9PBj0JgjC7mbGvyUqp3wLzbZ76RPGB1lorpfRM7UNrfRNwE8DGjRtn5HWUUsTM+eaZXN5z/EMQBOFEZsYERGt9mdNzSqmjSqkFWusepdQC4NgEfnQf0KyUCphWyGKge4rbnTLWfPN0Li/uK0EQ5gS1cmHdDVxrPr4WuMvrQq21Bh4E3j6Z9TNFNOQvuLDEAhEEYS5QKwH5PHC5UmoXcJl5jFJqo1Lq29ZFSqlHgB8DlyqlupRSV5hP/R3wMaXUboyYyM2v6O5tiAYNARkTAREEYY5QE1+L1roPuNTm/Gbg+qLjCx3W7wU2zdgGJ0Es5Dd7YeXEhSUIwpxAKtGniVgoQDydZSwlFoggCHMDEZBpIhbyM5bKkcjkPBcRCoIgnMiIr2WaaIwGGUlmSOfyMhhKEIQ5gQjINNEYCTKUyJDK5mmOiYAIgjD7EQGZJpqiQcbMNibNUREQQRBmPxIDmSaaisbRNomACIIwBxABmSZa68cnCjbHatpdXhAE4RVBBGSaKJ4oaE0ZFARBmM2IgEwTJQLSLAIiCMLsRwRkmlhYJBrtdWGXKwVBEGYHkoU1TQT8Pr7yrg2EAj58MgtEEIQ5gAjINPLWMx3nWgmCIMw6xIUlCIIgTAoREEEQBGFSiIAIgiAIk0IERBAEQZgUIiCCIAjCpBABEQRBECaFCIggCIIwKURABEEQhEmhtNa13sMrhlKqFzgwyeXtwPFp3M6JgLznucFce89z7f3C1N/zMq11R/nJOSUgU0EptVlrvbHW+3glkfc8N5hr73muvV+YufcsLixBEARhUoiACIIgCJNCBMQ7N9V6AzVA3vPcYK6957n2fmGG3rPEQARBEIRJIRaIIAiCMClEQARBEIRJIQJSBaXUlUqpHUqp3UqpG2q9n5lAKbVEKfWgUmq7UmqbUuqvzPOtSqn7lFK7zL9bar3X6UYp5VdKPaeU+oV5vEIp9aT5ef9QKRWq9R6nE6VUs1LqDqXUy0qpl5RS58/2z1kp9VHz3/VWpdT3lVKR2fY5K6VuUUodU0ptLTpn+7kqg/8w3/uLSqmzJvu6IiAuKKX8wNeBNwDrgHcrpdbVdlczQhb4P1rrdcB5wJ+b7/MG4H6t9WrgfvN4tvFXwEtFx18Avqy1XgUMANfVZFczx78Dv9ZarwXOwHjvs/ZzVkotAv4S2Ki1Xg/4gWuYfZ/zd4Ary845fa5vAFabfz4I/OdkX1QExJ1NwG6t9V6tdRr4AXB1jfc07Wite7TWz5qPRzBuKosw3utt5mW3AW+tzQ5nBqXUYuBNwLfNYwW8DrjDvGRWvWelVBNwEXAzgNY6rbUeZJZ/zhiju6NKqQAQA3qYZZ+z1vphoL/stNPnejXwXW3wBNCslFowmdcVAXFnEXCo6LjLPDdrUUotB84EngQ6tdY95lNHgM4abWum+Arwt0DePG4DBrXWWfN4tn3eK4Be4FbTbfdtpVQds/hz1lp3A18CDmIIxxDwDLP7c7Zw+lyn7b4mAiIUUErVA3cCf621Hi5+Thv53rMm51spdRVwTGv9TK338goSAM4C/lNrfSYwRpm7ahZ+zi0Y37hXAAuBOipdPbOemfpcRUDc6QaWFB0vNs/NOpRSQQzx+B+t9U/M00ct09b8+1it9jcDvBZ4i1JqP4Zr8nUY8YFm09UBs+/z7gK6tNZPmsd3YAjKbP6cLwP2aa17tdYZ4CcYn/1s/pwtnD7XabuviYC48zSw2szYCGEE3+6u8Z6mHdP3fzPwktb634qeuhu41nx8LXDXK723mUJr/X+11ou11ssxPtcHtNZ/DDwIvN28bLa95yPAIaXUyeapS4HtzOLPGcN1dZ5SKmb+O7fe86z9nItw+lzvBt5nZmOdBwwVubomhFSiV0Ep9UYMX7kfuEVr/bkab2naUUpdADwCbGE8HvD3GHGQHwFLMdrgv1NrXR6oO+FRSl0C/I3W+iql1EoMi6QVeA54j9Y6Vcv9TSdKqQ0YSQMhYC/wAYwvkrP2c1ZKfRZ4F0a24XPA9Rg+/1nzOSulvg9cgtG2/SjwGeBn2HyuppB+DcOVFwc+oLXePKnXFQERBEEQJoO4sARBEIRJIQIiCIIgTAoREEEQBGFSiIAIgiAIk0IERBAEQZgUIiCCUEPM7rgfrvU+BGEyiIAIQm1pBkRAhBMSERBBqC2fB05SSj2vlPpirTcjCBNBCgkFoYaY3Y9/Yc6qEIQTCrFABEEQhEkhAiIIgiBMChEQQagtI0BDrTchCJNBBEQQaojWug94TCm1VYLowomGBNEFQRCESSEWiCAIgjApREAEQRCESSECIgiCIEwKERBBEARhUoiACIIgCJNCBEQQBEGYFCIggiAIwqT4/9tIjx3D44BTAAAAAElFTkSuQmCC\n" }, "metadata": { "needs_background": "light" } } ], "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 279 }, "id": "OKU-D2jAXfxO", "outputId": "3dbf8cef-1a8e-48c3-bc81-a2f5291059f1" }, "execution_count": 1 }, { "cell_type": "markdown", "source": [ "For a second example, suppose that you want to solve the following\n", "coupled, second-order differential equations, \n", "\n", "$$\n", "\\frac{d^2x}{dt^2} = ay\n", "$$\n", "and\n", "$$\n", "\\frac{d^2y}{dt^2} = b + c\\frac{dx}{dt}.\n", "$$ \n", "If we make the definitions,\n", "$$\n", "z[0]=x,\n", "$$ \n", "$$\n", "z[1] = \\frac{dx}{dt},\n", "$$ \n", "$$\n", "z[2]=y,\n", "$$\n", "and\n", "$$\n", "z[3] = \\frac{dy}{dt},\n", "$$ \n", "then the two second-order equations can be written as the four first-order equations, \n", "\n", "$$\n", "\\frac{dz[0]}{dt} = z[1],\n", "$$ \n", "$$\n", "\\frac{dz[1]}{dt} = az[2],\n", "$$ \n", "$$\n", "\\frac{dz[2]}{dt} = z[3],\n", "$$ \n", "and\n", "$$\n", "\\frac{dz[3]}{dt} = b + cz[1].\n", "$$ \n", "These equations are now in a form necessary for the derivative function, w\n", "hich would be an array with four elements. Notice that the index of the array is not the number of derivatives of a single function in this case.\n" ], "metadata": { "id": "aFtor63MXfxQ" } }, { "cell_type": "code", "source": [ "from scipy.integrate import odeint\n", "import pylab as pl\n", "\n", "def deriv(z,t): # return derivatives of the array y\n", " a = -0.1\n", " b = 0.2\n", " c = 0.1\n", " return pl.array([z[1], a*z[2], z[3], b+c*z[1]])\n", "\n", "time = pl.linspace(0.0,10.0,1000)\n", "zinit = pl.array([0,0.2,0,0.1]) # initial values\n", "z = odeint(deriv,zinit,time)\n", "\n", "pl.figure()\n", "pl.plot(time,z[:,0],label='x') # z[:,0] is the first column of z\n", "pl.plot(time,z[:,2],label='y') # z[:,2] is the third column of z\n", "pl.xlabel('t')\n", "pl.legend(loc='upper left')\n", "pl.show()" ], "outputs": [ { "output_type": "display_data", "data": { "text/plain": [ "
" ], "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEGCAYAAABmXi5tAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXhV5bn38e+dnXmAkIQZAiiIAiJIBEFsnacqqG2djrMetK219m3f08G+6ml7Wqseq61a63TQU4c61ok61gkHJCDzLIOEMQkSkpCQYd/vH2ujiGHMTlaS/ftcV669pux1b8TfXjzrWc9j7o6IiHR8SWEXICIirUOBLyKSIBT4IiIJQoEvIpIgFPgiIgkiOewCdqegoMD79+8fdhkiIu3GjBkzyty9a1P72nTg9+/fn+Li4rDLEBFpN8xs1a72qUlHRCRBKPBFRBKEAl9EJEHsUxu+mT0EnA5sdPdhsW15wN+B/sBK4Bx3/7yJ370E+FVs9bfu/vD+FFxfX09JSQm1tbX78+stLj09nT59+pCSkhJ2KSIiX7GvN20nA3cBj+yw7efAm+5+s5n9PLb+sx1/KfalcCNQBDgww8xeaOqLYU9KSkrIycmhf//+mNm+/nqLcnfKy8spKSlhwIABYZcjIvIV+9Sk4+7vApt22jwR2H61/jBwZhO/ejLwurtvioX868Ap+1grALW1teTn57e5sAcwM/Lz89vsvz5EJLHFow2/u7uviy2vB7o3cUxvYPUO6yWxbV9jZpPMrNjMiktLS5s8YVsM++3acm0iktjietPWg7GWmzXesrvf5+5F7l7UtWuTzw6IiHRM7rD0DZh6R4u8fTwCf4OZ9QSIvW5s4pg1QN8d1vvEtomISGMDzH0a7j0aHv02FD8I9fFvGo5H4L8AXBJbvgR4voljXgVOMrMuZtYFOCm2TUQkcdXXwPQH4M+HwzNXQOM2mHgPXDMDUtLjfrp9Cnwzexz4EBhsZiVmdgVwM3CimS0FToitY2ZFZvYAgLtvAn4DTI/9/Dq2rd2ZPn06w4cPp7a2lurqaoYOHcq8efPCLktE2pOazfDubXDHofDyTyCrK5z7KHx/Goz8N0hObZHT7lO3THc/fxe7jm/i2GLgyh3WHwIe2qfq9uA/X5zPgrVb4vmWDOnViRvPGLrL/UcccQQTJkzgV7/6FTU1NVx44YUMGzYsrjWISAe1ZR18dDcUT4a6Shh4Ioy/DvodBa3Q4aNND57WVt1www0cccQRpKen86c//SnsckSkrStbBh/cCbOfgGgDDD07CPoeh7ZqGe068Hd3Jd6SysvLqaqqor6+ntraWrKyskKpQ0TauPXz4L3bYP4/IDkNDr8Yxl4DeeE8mNmuAz8sV111Fb/5zW9YsWIFP/vZz7jrrrvCLklE2pI1M4I2+sVTIDUHxv8Yjvw+ZIfb1VyBv48eeeQRUlJSuOCCC2hsbGTcuHH861//4rjjjgu7NBEJ22cfwTu3wKdvQnouHHs9jP53yOgSdmWAAn+fXXzxxVx88cUARCIRpk2bFnJFIhIqd1jxTnBFv/I9yCyAE26CI66EtJywq/sKBb6IyP5wh6Wvw7u3QsnHkN0DTv49jLoUUjPDrq5JCnwRkX0RjcLil4OgXzcbOveFb/03jLiwRR6WiicFvojI3og2wvzn4L3/ho0LoMsAmHAXHHYeRNrH/BcKfBGR3WlsgLlPBd0ry5dBwWA4+/6gL32kfUVo+6pWRKS1NDbAvKeDXjebPoXuh8J3H4ZDJkBS+5wdVoEvIrKjaCPMewbe+UNwRd/90GCcm4O/1SrDH7QkBb6ICMSC/tlY0C+FbkPhnP+Fg09vt1f0O1Pgi0hi234z9p1boGwxdBsC5zwCB5/RYYJ+OwX+PrrhhhvIy8vjuuuuA+D666+nW7du/OhHPwq5MhHZJ9EoLPhHcEVfugi6HgLfnQyHTOxwQb9d+w78f/4c1s+N73v2OBROvXmXuy+//HLOPvtsrrvuOqLRKE888QQff/xxfGsQkZYTjcLC5+HtP0DpwqDXzXcegiFnddig3659B34I+vfvT35+Pp988gkbNmxg5MiR5Ofnh12WiOxJNAqLXgyCfuN8KDgIvv0gDD0LkiJhV9cq2nfg7+ZKvCVdeeWVTJ48mfXr13P55ZeHUoOI7KVoFBa9FDTdbJgH+YPg7Adg2NkJE/Tbte/AD8lZZ53FDTfcQH19PY899ljY5YhIU9xh0cvwzs1B02/+wOCBqWHfTrig367ZgW9mg4G/77DpAOAGd79jh2OOIZjcfEVs07Pu/uvmnjssqampHHvsseTm5hKJJOZfHJE2yx0W/xPe/j2snwN5B8BZf4Vh32l3T8bGW7M/vbsvBkYAmFkEWAM818Sh77n76c09X1sQjUb56KOPeOqpp8IuRUS2c4clrwRBv252MNbNmffCod9N+KDfLt5/CscDn7r7qji/b5uxYMECTj/9dM466ywGDRoUdjki4g5LXwuCfu0n0KU/TLwHhp+roN9JvP80zgMe38W+sWY2G1gL/NTd58f53K1iyJAhLF++POwyRGT7ePRv/x7WzoTcfjDx7ljQt4/RK1tb3ALfzFKBCcAvmtg9E+jn7lVmdhrwD6DJy2MzmwRMAigsLGzyXO6OtdExLdw97BJEOjZ3WPZmEPRriiG3ECb8GQ47X0G/B/F8yuBUYKa7b9h5h7tvcfeq2PIUIMXMCpp6E3e/z92L3L2oa9evT/ibnp5OeXl5mwxWd6e8vJz09LY9CYJIu7Q96B88ER79NlRtgDPuhGtmwOEXK+z3QjybdM5nF805ZtYD2ODubmajCb5oyvfnJH369KGkpITS0tL9r7QFpaen06dPn7DLEOk43GH528EV/epp0KkPnP7HYIap5NSwq2tX4hL4ZpYFnAhctcO2qwHc/V7gO8D3zKwBqAHO8/28RE9JSWHAgAHNL1pE2jZ3WPFuEPSffQidesO3boeRF0JyWtjVtUtxCXx3rwbyd9p27w7LdwF3xeNcIpIAVk6Ft34Hq96HnJ5w2m1Bs42CvlnUZ0lE2o5VHwRBv/I9yO4Bp94Ch1/S5icHby8U+CISvs8+CoJ+xTuQ3R1OuRlGXQopGWFX1qEo8EUkPKs/DoJ++VuQ1RVO/h2MugxSM8OurENS4ItI6yspDoL+0zchswBO+i0UXaGgb2EKfBFpPWtmwFu/h2WvQ0YenPCfMPrfITUr7MoSggJfRFre2k/g7ZuDwc0yusDxN8LoSZCWHXZlCUWBLyItZ93sIOgXT4H0XDjuVzD6KkjvFHZlCUmBLyLxt35uEPSLXoL0znDs9TDmqmBZQqPAF5H4WT8vmEpw4QuQ1hmO+QWMuRoycsOuTFDgi0g8rJsN79wSXNGn5sA3/gPGfj9or5c2Q4EvIvtvzQx451ZY8s/giv6bP4cjr1bQt1EKfBHZd6unB003y14PbsYe+ysYM0lt9G2cAl9E9t6qD4OgX/5W0I/++BvhiCvV66adUOCLyJ6teC8I+pXvBU/Gnvjr4MlY9aNvVxT4ItI092Aws3duCYYpzu6usW7aOQW+iHyVezDGzTu3BDNM5fSMDVN8sUavbOcU+CIScIelrwVBv6Y4mGHqtNtg5EUaj76DUOCLJLpoIyx4Ht67HTbMhc6FcPodMOICzTDVwSjwRRJVQx3MfRKm/hHKl0H+QJh4Nxx6jiYH76DiFvhmthKoBBqBBncv2mm/AXcCpwFbgUvdfWa8zi8ie6m+BmY+Au//CbaUQI9D4buT4ZAJkBQJuzppQfG+wj/W3ct2se9UYFDsZwzwl9iriLSG2gqY/iB8dA9Ul0LfI+GMO2DgCWAWdnXSClqzSWci8Ii7O/CRmeWaWU93X9eKNYgknuoy+Ogv8PH9sK0iCPijfwL9xoVdmbSyeAa+A6+ZmQN/dff7dtrfG1i9w3pJbNtXAt/MJgGTAAoLC+NYnkiCqVgDH94FMyYHzThDJsD4/wO9RoRdmYQknoE/3t3XmFk34HUzW+Tu7+7rm8S+KO4DKCoq8jjWJ5IYypbBB3fCrMfBozD8XBj/Y+h6UNiVScjiFvjuvib2utHMngNGAzsG/hqg7w7rfWLbRCQeVk8Pgn7hSxBJhVGXwrgfQpd+YVcmbURcAt/MsoAkd6+MLZ8E/Hqnw14ArjGzJwhu1lao/V6kmaJRWPpq0OPmsw+CkSu/8dNgvtjsbmFXJ21MvK7wuwPPBT0vSQYec/dXzOxqAHe/F5hC0CVzGUG3zMvidG6RxNOwDeY8CR/8GcoWQ+e+cMofYOSFGtBMdikuge/uy4HDmth+7w7LDvwgHucTSVi1FVD8P0Gvm6r1QR/6sx+AoWdCJCXs6qSN05O2Iu1BxRqY9hcongx1lXDAMXDWX+CAY9WHXvaaAl+kLduwIOhaOedJ8EYYejYcdS30/No/qEX2SIEv0tZEo8HUgR/dA8vfhpRMKLo8mBS8S/+wq5N2TIEv0lbUVcPsx+Gje6F8KeT0CqYQHHUpZOaFXZ10AAp8kbBVrIGP7wueiK3dDL1G6kastAgFvkhY1syAD++BBf8Inog9+HQY+wPoO0Y3YqVFKPBFWlNjAyx6KWifXz0NUnNg9FUwZpLa56XFKfBFWkN1WTAGffFDULEacvvBKTfDiH+D9E5hVycJQoEv0pJKZsD0+2HeM9BYBwO+Aaf8HgafpslGpNUp8EXirb4W5j8bjD+/diakZsPhl8ARV0K3g8OuThKYAl8kXjZ/FswoNfMRqNkEBQfBabcFwxOr2UbaAAW+SHNEo7Dibfj4AVjyz2Db4NOC0SoHfEO9baRNUeCL7I+qUpj1N5jxMHy+AjILgklGRl0GuX33/PsiIVDgi+ytaBRWvBM8ILXoZYjWQ7/xcOwvYchESE4Lu0KR3VLgi+xJ5QaY9SjMfBg+XwkZeTDmquBGrKYNlHZEgS/SlGgUlr8VXM0vngLRBuh/NBz3/4InYlPSw65QZJ8p8EV2tPmzYPLvWY/C5lWxq/mrgwHMCgaFXZ1IsyjwReq2wsIXg5Bf8S7gQQ+b42+AQ85Q27x0GM0OfDPrCzxCMK+tA/e5+507HXMM8DywIrbpWXffeZJzkdbjDqs/DnrazHsumEUqtx8c8ws47Dzo0i/sCkXiLh5X+A3AT9x9ppnlADPM7HV3X7DTce+5++lxOJ/I/tuyNhhzftZjUL4smFxkyJkw4gLodxQkJYVdoUiLaXbgu/s6YF1sudLMFgK9gZ0DXyQc2yqDbpRzngxuxHoUCsfBUdcFY86n5YRdoUiriGsbvpn1B0YC05rYPdbMZgNrgZ+6+/xdvMckYBJAYWFhPMuTRNJYD8vehLlPwqIp0FADnQth/P8JrubzDwy7QpFWF7fAN7Ns4BngOnffstPumUA/d68ys9OAfwBNdnlw9/uA+wCKioo8XvVJAnAPxpif8yTMfy4YzyajC4w4Hw49J5hYRE02ksDiEvhmlkIQ9o+6+7M779/xC8Ddp5jZPWZW4O5l8Ti/JLiNi4Ir+blPBd0qkzNg8Kkw/Bw48HhITg27QpE2IR69dAx4EFjo7rfv4pgewAZ3dzMbDSQB5c09tySw0sUw/x/BlXzpQrAkOOBYOPZ6OPhbapcXaUI8rvCPAi4C5prZrNi2XwKFAO5+L/Ad4Htm1gDUAOe5u5prZN9sXBTM/zr/H0HIY1A4Fk69BYaeBdndwq5QpE2LRy+dqcBux4B197uAu5p7LklAX4T8c1C6iC9D/tbgoahOPcOuUKTd0JO20ra4B7NELZoSTPa9PeT7jVPIizSTAl/C11AHK98NQn7xP6FybdAmXxgL+SETIKdH2FWKtHsKfAlHbQUsfT14IGrp68HQBimZcOBxcPANcNDJkJkXdpUiHYoCX1qHO5R/CstehyWvwsr3giGHs7rCsLNg8LfggG9CSkbYlYp0WAp8aTl1W2Hl1CDkl74WTB4CkD8Ixv4gCPk+RZAUCbVMkUShwJf4Kv80aKJZ+hqseh8aaoMHoQZ8A8ZeA4NOhC79w65SJCEp8KV5qsuDG67L34HlbwcTegPkDwwm9B50YjAKpWaIEgmdAl/2TV01rPoQVrwdhPz6uYBDag70PwqO/D4MOgHyDgi7UhHZiQJfdq9hG6yZGcwEteKdYNKQaD0kpQSDkR37SzjgGOg1EiIpYVcrIruhwJev2lYZjDi56kP47EMoKYbGbYBBz+Fw5PeC3jSFYyE1K+xqRWQfKPATXXVZEOyrPgh+1s8JJgixSBDwo/89CPd+49QvXqSdU+AnkoZtsH4erCkOrtxLpn95kzU5HfocAUf/FPqNhT6jIS073HpFJK4U+B2Ve9Dvfc2MINhLioOr98a6YH9OT+g9CkZdGlzB9xoByWlhViwiLUyB3xFEG4MJudfNgfWzY69zoObzYH9yRnBTdcxVwVV87yLo3DvcmkWk1Snw25u6rVC2+MtQXzcbNsyH+q3B/kgqdBsSjCrZc0TwJGu3IepBIyIK/DarvgbKlgTjwZcu/PL181VAbO6Y1BzocSgcfklwg7XHcOg6eL/DPRp1auob2VrXSG19I3WNUeoaYj+NUeobomzbcVtse11DlPrGKA1RJ+qOe/BeUSe2/uVy1MEJjjGD5CQjYkYkKYlIEl99NYhEkkiNGOkpkS9/kpPISN2+HCE9NYn0lAhZqclEknY7NYNIQlPgh8kdKtfDpuWw6dNgWIKypbFgXxn0lgFISg6eXO05Ag47H7oeDD0OZVunQiq3RamsbaCytp6qqga2lJUHy9sa2FrXSHXstaauka31jWyNrW+tb6SmLrZc18jWugZq66Mt8jHNIMmMJAPDMAu2RWNfDA3R+E1+lpOWTE56Mp0yUuiUnvLFck56Mp3SU+iSlUpBdir5WWkU5ASveVmp+qKQhKDAb2mNDcH47ptXx4I9CHcv/xQ2rcDqq784NGrJVGb2pTzzADb0Po6SlH6ssEJWeHc2bzOqNjRQuao+FvBLqGtctMfTJycZGakRMlMjZKYmx14jdM5IoWendDJTI2SkRshKSyYjJfLF/rSUCGnJSaRGkkiJJJGavMNPJCnYt8N6SnISyUkWC/Yg3JNse7jvOUy3B3/Ug9fG2E9DNEp9o1NbH3xpbWtopLY+GqzXf7lcW98Y+3NpYEttPVtqgj+n9VtqWbKxMtheU09T3y1mkJeZSn52KgXZafTsnEHv3HR65mbQKzeDXp2D5ew0/e8i7Vtc/gab2SnAnUAEeMDdb95pfxrwCDCKYPLyc919ZTzOHSp36qo3U1X6GbVlK6kv/wyvWE3SlhJSq9eSuXUt2XWlJPHllXMDEdbQjeWN3VnhR7PSu7PSe7DCe7DWC2isiUB5EELZacmxq9R6stOSKchOZUBBFjnpyeTErl6/+ElLITu23Ck9hey0ZLLSkklNTgrxD2jvJSUZqS18le3ubKlpoLRqG+VV2yivrqOsahtlVXWUV22jrGobpZXb+ODTMjZsqf3al0On9GR6d8mkf34m/fKzGFCw/TWLbjlpe/XFJhKmZge+mUWAu4ETgRJgupm94O4LdjjsCuBzdx9oZucBfwDObe6548Xdqa5rpKKmns1VNVRXlLO1ooxtleU0Vm6E6g2kVG8ktbaMjLpycurL6Nz4OXn+OelWz46PI9V7hPWex0oKWOuDKIuMpyK1O9XpPajK6kdDTm9yMjPIzUyhc0YKwzNTOToj5Yv17WGemRIhSc0McWVmdM5MoXNmCgO77f4Zg4bGKBsqt7Fucw1rK2pZu7mGdZtrKPm8hsUbKnlj4QbqG7/8RshIidAvP5MDu2UzuHsOg3vkMLh7Dn3zMtVcJG1GPK7wRwPL3H05gJk9AUwEdgz8icBNseWngbvMzNw9fo23Me5OydLZVFdVUFO1hdrqLdRt3UJdTSWNtVV4bSVeV43VV5FSX0lG4xayGyvpTBW5VkVPtpJkTZdVYTlURPKpSsljXXY/VqZ3pTGzazD9Xm4hKXmFZHTpRW52OodkpDI6PVmh3U4lR5LonZtB79ymJ2RpjDprN9ewoqyaVeXVrCjbysryauaUbOblOeu+OC49JYmDuudwUPccDu6Rw/A+uQzt1YksNQ9JCOLxt643sHqH9RJgzK6OcfcGM6sA8oGyOJz/a7o+egJ9rX6X++tIYVtSBrWRbGrTO1Of2o369IMoTc+lLDOf5Ow8UnPySe9UQGZud9K69MSyutE5OZXOLVGwtDuRJKNvXiZ98zKBrl/ZV72tgaUbq1i8fguL11exZEMlby8u5ekZJUDQXDewazaH9unM8N6dGd43lyE9O5GeoolgpGW1ucsMM5sETAIoLCzcn99nwdjbSE5JIz27E1lZncnqlEt2TmeS03MgNYvUSAqpQE6caxcByEpLZkTfXEb0zf3K9o2VtcxbU8GckgrmllTw7pIynp25BoCUiDGsd2dG98+jqH8eRf260CUrNYzypQOz5raqmNlY4CZ3Pzm2/gsAd//9Dse8GjvmQzNLBtYDXffUpFNUVOTFxcXNqk+krXJ31m+pZU5JBZ98tpnilZuYU1JBXWNwk39Qt2yOGJDHmAF5jDuwgK45GvpC9szMZrh7UVP74nGFPx0YZGYDgDXAecAFOx3zAnAJ8CHwHeBfLdF+L9KemBk9O2fQs3MGJw/tAUBtfSNzSiqYvnIT01du4sVZa3ls2mcAHNKzE0cPKmD8wAJGD8hTE5Dss2Zf4QOY2WnAHQTdMh9y9/8ys18Dxe7+gpmlA/8LjAQ2Aedtv8m7O7rCl0TXGHXmr63gvaVlTF1aRvGqTdQ3OqnJSYzun8cxg7tywiHd6V+guQkksLsr/LgEfktR4It81da6Bqat2MTUpWW8t7SUJRuqABjYLZsTh3TnhEO6M7JvrnqHJTAFvkgHtXrTVt5YuIE3Fm5g2vJNNESdguxUjj+4O98a3pNxB+aTHGkfD99JfCjwRRJARU097ywp5fUFG3hr0UaqtjWQn5XKKcN6cPrwXowekKeHwBKAAl8kwdTWN/LOklJenL2WNxdupKa+kW45aZx2aE/OGtmb4X06ayiIDkqBL5LAttY18ObCjbw0Zy1vLS6lriHKQd2z+e6ovpw5sre6e3YwCnwRAYJmn5fmrOWp4hJmrd5McpJxzOBunFPUh2MP7kaK2vvbPQW+iHzN0g2VPD2jhGdmrqGsahvdctI4f3QhF4wppHun9LDLk/2kwBeRXWpojPLW4lIenbaKtxeXkpxknDy0BxeN7ceYAXlq629nWvpJWxFpx5IjSZw4pDsnDunOqvJq/vbRKp4sLuHlues4qHs2l4zrz7cP76MnezsAXeGLyNfU1DXy4py1PPLhSuat2UJ+ViqXjOvPRUf206BubZyadERkv7g7Hy3fxH3vfspbi0vJSIlw7hF9uWL8gNjQ0NLWKPBFpNkWr6/kvneX88LsNTRGnTMO68UPjxu0x9nDpHUp8EUkbtZX1PLQ+yv420erqKlv5Izhvbj2+IEM7KYZJtoCBb6IxN2m6jruf285D3+wkpr6Rk4f3otrjxvIoO4K/jAp8EWkxWyqruOBWPBvrW9k4mG9+MlJg9XGHxIFvoi0uM+r67jvveX8z/sriEbhwiP7cc1xA8lTr55WpcAXkVazvqKWO95YwpPFq8lKTebqYw7ksqP6k5mqx35aw+4CXwNniEhc9eiczs3fHs5rP/4GRx6Yz62vLuaYW9/m6RklRKNt9wIzESjwRaRFDOyWw/0XF/H01WPp3SWDnz41m7P+8gGzVm8Ou7SEpcAXkRZV1D+PZ64ex+3nHMa6zTWceff7/PSp2WysrA27tITTrEY1M7sVOAOoAz4FLnP3r319m9lKoBJoBBp21b4kIh1TUpJx9uF9OGloD+761zIenLqcV+at59rjB3LZUQM0LHMrae6f8uvAMHcfDiwBfrGbY4919xEKe5HElZ2WzM9PPZjXfvxNxgzI43dTFnHGn6fyyWefh11aQmhW4Lv7a+7eEFv9COjT/JJEpKMbUJDFg5cewX0XjWLz1nrO/ssH3Pj8PCpr68MurUOL57+jLgf+uYt9DrxmZjPMbNLu3sTMJplZsZkVl5aWxrE8EWlrThragzd+8k0uGdufRz5axYm3v8sr89aHXVaHtcd++Gb2BtCjiV3Xu/vzsWOuB4qAs72JNzSz3u6+xsy6ETQD/dDd391TceqHL5I4Zq3ezC+encvCdVs4eWh3fnvmoZpvdz+06INXZnYpcBVwvLtv3YvjbwKq3P22PR2rwBdJLPWNUR6cuoLbX19CVmqE3555KN8a3jPsstqVFnvwysxOAf4DmLCrsDezLDPL2b4MnATMa855RaRjSokkcfU3D2TKteMpzMvkB4/N5JrHZvJ5dV3YpXUIzW3DvwvIAV43s1lmdi+AmfUysymxY7oDU81sNvAx8LK7v9LM84pIBzawWw7PfG8c//fkwbw6fz0n/vFdXl+wIeyy2j2NpSMibdrCdVv4yZOzWbBuC98d1YebJgwlK03j8uyKxtIRkXbrkJ6d+McPjuKaYwfyzMwSTv/zVOatqQi7rHZJgS8ibV5qchI/PXkwj/37kdTWN3LWPe/zwHvLNRjbPlLgi0i7ceQB+fzzR0dz3MHd+O3LC7l08nRKK7eFXVa7ocAXkXYlNzOVey8cxW/PHMa05eWceue7fLCsLOyy2gUFvoi0O2bGhUf244VrxtM5I4ULH5zGPW8vUxPPHijwRaTdGtwjh+evGc9ph/bkllcWM+l/Z1BRo/F4dkWBLyLtWnZaMn8+fyQ3njGEtxdvZMJdU5m/Vr14mqLAF5F2z8y47KgB/P2qoBfP2fd8wDMzSsIuq81R4ItIhzGqXx4vX3s0hxd24SdPzea/Xl5Ao9r1v6DAF5EOpSA7jUeuGM0lY/tx/3sruHzydLZonH1AgS8iHVBKJIn/nDiM3511KO8vK+PMu99nRVl12GWFToEvIh3WBWMK+duVY/i8uo6Jd03lvaWJPamSAl9EOrQjD8jnhWvG0ys3g0se+phHp60Ku6TQKPBFpMPrm5fJ098bxzcP6sr1z83jllcWJeRDWgp8EUkI2WnJ3H9xEReMKeSetz/lx0/OYltDY9hltSoNKi0iCZiStoQAAAqASURBVCM5ksR/nTmMPl0yuOWVxayvqOW+i4ronJkSdmmtQlf4IpJQzIzvHzOQO88bwSefbebb937A6k17nI67Q1Dgi0hCmjiiN49cMZqNW2r5zr0fsGRDZdgltbjmTmJ+k5mtic1nO8vMTtvFcaeY2WIzW2ZmP2/OOUVE4uXIA/J56upxuMM5f/2QWas3h11Si4rHFf4f3X1E7GfKzjvNLALcDZwKDAHON7MhcTiviEizDe6Rw9NXj6NTegoX3P8RU5d23LH1W6NJZzSwzN2Xu3sd8AQwsRXOKyKyVwrzM3n66rEU5mVy+eTpvDJvXdgltYh4BP41ZjbHzB4ysy5N7O8NrN5hvSS2rUlmNsnMis2suLQ0sZ+KE5HW061TOn+fNJZhvTvx/Udn8uT01Xv+pXZmj4FvZm+Y2bwmfiYCfwEOBEYA64D/bm5B7n6fuxe5e1HXrl2b+3YiInutc2YKf7tyDOMHdeU/npnDwx+sDLukuNpjP3x3P2Fv3sjM7gdeamLXGqDvDut9YttERNqczNRkHri4iGsem8mNL8ynIepcMX5A2GXFRXN76fTcYfUsYF4Th00HBpnZADNLBc4DXmjOeUVEWlJqchJ3/9vhnDqsB795aQH3vftp2CXFRXPb8G8xs7lmNgc4FvgxgJn1MrMpAO7eAFwDvAosBJ509/nNPK+ISItKiSTxp/NHcvrwnvxuyiLufmtZ2CU1W7OGVnD3i3axfS1w2g7rU4CvddkUEWnLUiJJ3HHuCJKTjFtfXUxDo/OjEwaFXdZ+01g6IiK7kRxJ4r/PGUFSkvHHN5bgONedcFDYZe0XBb6IyB5Ekoxbv3MYSWbc8cZS0pIjfO+YA8Mua58p8EVE9kIkyfjDt4ezrSHKH15ZREZKEpce1b567yjwRUT2UiTJuP2cw9hW38hNLy4gLSXC+aMLwy5rr2m0TBGRfZASSeLPF4zkmwd15ZfPzeW5T0rCLmmvKfBFRPZRWnKEv140iiMH5POTJ2czZW77GHtHgS8ish/SUyI8cEkRIwu7cO3jn/DOkrY/9pcCX0RkP2WlJfM/lx3BoO45fO9vM/jks8/DLmm3FPgiIs3QKT2Fhy8/goLsNC6fPJ1lG6vCLmmXFPgiIs3ULSed/71iNJEk4+IHp7GuoibskpqkwBcRiYN++VlMvmw0W2obuPjBj9m8tS7skr5GgS8iEifDenfmvotHsap8K5dPnk5NXWPYJX2FAl9EJI7GHVjAn84fwazVm/nh4zNpjHrYJX1BgS8iEmenDOvJf04YyhsLN/KblxaEXc4XNLSCiEgLuGhsf1aVb+WBqSuCydHbwKxZCnwRkRbyy9MOYfXnW/nNywvo0yWDk4b2CLUeNemIiLSQpCTjjnNHMrxPLtc+8QmzV28Ot55Qzy4i0sFlpEZ44OIiCrLTuOLhYlZv2hpaLc2dxPzvZjYr9rPSzGbt4riVsblvZ5lZcXPOKSLS3nTNSWPyZUdQ19DI5ZOns6W2PpQ6mhX47n6uu49w9xHAM8Czuzn82NixRc05p4hIezSwWw73XjSKFWXVXPfErFC6a8alScfMDDgHeDwe7yci0hGNO7CAGycM5V+LNnLrq4tb/fzxasM/Gtjg7kt3sd+B18xshplNitM5RUTanYuO7Me/jSnk3nc+bfXJU/bYLdPM3gCa6kt0vbs/H1s+n91f3Y939zVm1g143cwWufu7uzjfJGASQGFh+5k6TERkb900YSifllbxs2fmMqAgmxF9c1vlvObevHYkM0sG1gCj3H2PX1dmdhNQ5e637enYoqIiLy7WPV4R6Xg2Vdcx4a6p1DVEefGH4+neKT0u72tmM3Z1rzQeTTonAIt2FfZmlmVmOduXgZOAeXE4r4hIu5WXlcoDlxRRta2BSY8UU1vf8gOtxSPwz2On5hwz62VmU2Kr3YGpZjYb+Bh42d1ficN5RUTatYN7dOKP545gdkkFv3x2Ls1tcdmTZg+t4O6XNrFtLXBabHk5cFhzzyMi0hGdPLQHPz7hIP74xhJGFOZy8dj+LXYuPWkrIhKyHx43kOMP7savX1xA8cpNLXYeBb6ISMiSkozbzx1B7y4ZfO/RmWzcUtsy52mRdxURkX3SOSOFv140iqraBr7/6EzqGqJxP4cCX0SkjTi4Ryf+8J3hDOqegxP/G7gaD19EpA2ZcFgvJhzWq0XeW1f4IiIJQoEvIpIgFPgiIglCgS8ikiAU+CIiCUKBLyKSIBT4IiIJQoEvIpIgmj0BSksys1Jg1X7+egFQFsdy2gN95o4v0T4v6DPvq37u3rWpHW068JvDzIp3NetLR6XP3PEl2ucFfeZ4UpOOiEiCUOCLiCSIjhz494VdQAj0mTu+RPu8oM8cNx22DV9ERL6qI1/hi4jIDhT4IiIJosMFvpmdYmaLzWyZmf087Hpampn1NbO3zGyBmc03sx+FXVNrMbOImX1iZi+FXUtrMLNcM3vazBaZ2UIzGxt2TS3NzH4c+3s9z8weN7P0sGuKNzN7yMw2mtm8HbblmdnrZrY09tolHufqUIFvZhHgbuBUYAhwvpkNCbeqFtcA/MTdhwBHAj9IgM+83Y+AhWEX0YruBF5x94OBw+jgn93MegPXAkXuPgyIAOeFW1WLmAycstO2nwNvuvsg4M3YerN1qMAHRgPL3H25u9cBTwATQ66pRbn7OnefGVuuJAiB3uFW1fLMrA/wLeCBsGtpDWbWGfgG8CCAu9e5++Zwq2oVyUCGmSUDmcDakOuJO3d/F9i00+aJwMOx5YeBM+Nxro4W+L2B1Tusl5AA4bedmfUHRgLTwq2kVdwB/AcQDbuQVjIAKAX+J9aM9YCZZYVdVEty9zXAbcBnwDqgwt1fC7eqVtPd3dfFltcD3ePxph0t8BOWmWUDzwDXufuWsOtpSWZ2OrDR3WeEXUsrSgYOB/7i7iOBauL0z/y2KtZuPZHgy64XkGVmF4ZbVevzoO98XPrPd7TAXwP03WG9T2xbh2ZmKQRh/6i7Pxt2Pa3gKGCCma0kaLY7zsz+Fm5JLa4EKHH37f96e5rgC6AjOwFY4e6l7l4PPAuMC7mm1rLBzHoCxF43xuNNO1rgTwcGmdkAM0sluMHzQsg1tSgzM4J23YXufnvY9bQGd/+Fu/dx9/4E/43/5e4d+srP3dcDq81scGzT8cCCEEtqDZ8BR5pZZuzv+fF08BvVO3gBuCS2fAnwfDzeNDkeb9JWuHuDmV0DvEpwR/8hd58fclkt7SjgImCumc2Kbfulu08JsSZpGT8EHo1dzCwHLgu5nhbl7tPM7GlgJkFvtE/ogMMsmNnjwDFAgZmVADcCNwNPmtkVBEPEnxOXc2loBRGRxNDRmnRERGQXFPgiIglCgS8ikiAU+CIiCUKBLyKSIBT4IvsgNmLl98OuQ2R/KPBF9k0uoMCXdkmBL7JvbgYONLNZZnZr2MWI7As9eCWyD2Ijkr4UG59dpF3RFb6ISIJQ4IuIJAgFvsi+qQRywi5CZH8o8EX2gbuXA+/HJtXWTVtpV3TTVkQkQegKX0QkQSjwRUQShAJfRCRBKPBFRBKEAl9EJEEo8EVEEoQCX0QkQfx/7W3IGXsJoIUAAAAASUVORK5CYII=\n" }, "metadata": { "needs_background": "light" } } ], "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 279 }, "id": "ClguseKVXfxR", "outputId": "3a2847f2-eacf-4d42-f7cb-ce928bb367f7" }, "execution_count": 2 }, { "cell_type": "markdown", "source": [ "##Exercise" ], "metadata": { "id": "-Ol3efUMXfxS" } }, { "cell_type": "markdown", "source": [ "An example of a differential equation that exhibits chaotic behavior is \n", "$$\\frac{d^3x}{dt^3} = -2.017 \\frac{d^2x}{dt^2} + \\left(\\frac{dx}{dt}\\right)^2 - x.$$\n", "
    \n", "
  1. Write the differential equation as a set of three first-order\n", "differential equations.\n", "
  2. Modify the example program to solve the equations with the initial conditions of $x=0$, $dx/dt = 0$, and $d^2x/dt^2 = 1$.\n", "
  3. Plot the results for $x$ for $t$ from 0 to 100.\n", "
" ], "metadata": { "id": "K9Qovb3jXfxS" } }, { "cell_type": "markdown", "source": [ "##Additional Documentation" ], "metadata": { "id": "HE5fV0XDXfxU" } }, { "cell_type": "markdown", "source": [ "Further information is available at: \n", "http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html " ], "metadata": { "id": "cKF44YAWXfxU" } } ] }