{ "metadata": { "name": "", "signature": "sha256:9150ebd94535c4ffe58299306d649d13cb2cca95f73fd2d87600f857292c4097" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "heading", "level": 1, "metadata": {}, "source": [ "Solving Ordinary Differential Equations (ODEs) using Python" ] }, { "cell_type": "markdown", "metadata": {}, "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", "\\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", "\\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", "\\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." ] }, { "cell_type": "code", "collapsed": false, "input": [ "from scipy.integrate import odeint\n", "from pylab import * # for plotting commands\n", "\n", "def deriv(x,t): # return derivatives of the array y\n", " a = -2.0\n", " b = -0.1\n", " return array([ x[1], a*x[0]+b*x[1]])\n", "\n", "time = linspace(0.0,100.0,1000)\n", "xinit = array([0.0005,0.2]) # initial values\n", "x = odeint(deriv,xinit,time)\n", "\n", "figure()\n", "plot(time,x[:,0]) # x[:,0] is the first column of x\n", "xlabel('t')\n", "ylabel('x')\n", "show()" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEPCAYAAACKplkeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4VOXdPvB7MEF2CUsCWTQsgbAbCCBabSwGRCVuvIqi\nUqUWUbSvSwu/vlpDWwXq60KlrUutYhfE11ZBxQjWBiuIqQKKBCQgkZCNJWxhSSB5fn98Pcxkcs6c\n58xMMsmc+3NdXGSWJ3Nyksyd77Mdj1JKgYiIKAzaRPoAiIgoejBUiIgobBgqREQUNgwVIiIKG4YK\nERGFDUOFiIjCJqKhkpeXh/T0dKSlpWHhwoWNHt+2bRvGjRuHdu3a4cknn2zwWGpqKoYPH46MjAyM\nGTOmuQ6ZiIgCiInUC9fV1WH27Nn44IMPkJSUhNGjRyMnJweDBg0685zu3bvj2WefxVtvvdWovcfj\nQX5+Prp169ach01ERAFErFIpKChA//79kZqaitjYWEydOhXLly9v8JyePXsiMzMTsbGxpp+D6zaJ\niFqWiIVKaWkpUlJSztxOTk5GaWmpdnuPx4PLLrsMmZmZePHFF5viEImIyKGIdX95PJ6Q2q9duxa9\ne/fGvn37kJ2djfT0dFx88cVhOjoiIgpGxEIlKSkJJSUlZ26XlJQgOTlZu33v3r0BSBfZtddei4KC\ngkahEmpwERG5VbDDCxHr/srMzERRURGKi4tRW1uLZcuWIScnx/S5/l/c8ePHcfToUQDAsWPHsGrV\nKgwbNsyyLf8pPProoxE/hpbyj+eC54LnIvC/UESsUomJicHixYsxceJE1NXVYcaMGRg0aBCef/55\nAMDMmTNRUVGB0aNH48iRI2jTpg0WLVqEwsJC7N27F9dddx0A4PTp05g2bRomTJgQqS+FiIi+E7FQ\nAYBJkyZh0qRJDe6bOXPmmY979erVoIvM0KlTJ2zatKnJj4+IiJzhinqXyMrKivQhtBg8F148F148\nF+HhUaF2oLVgHo8n5P5BIiK3CeW9k5UKERGFDUOFiIjChqFCRERhw1AhIqKwYagQEVHYuCJUlAL6\n9QP++c9IHwkRUXRzRaiUlwPffAO8/36kj4SIKLq5IlSKi+X/nTsjehhERFHPNaGSnAyUlUX6SIiI\nopsrQqWiAhg1iqFCRNTUXBEqVVXAkCFAZWWkj4SIKLq5IlQOHACSkuTjEycieyxERNHMFaFSVQV0\n6wbExQEHD0b6aIiIoperQqVrV+DQoUgfDRFR9HJFqFRXA507s1IhImpqrgmVjh1ZqRARNTVXhMqx\nY0CnThIsx45F+miIiKKXK0LFqFQYKkRETcsVoWJUKp06MVSIiJpS1IdKXR1w8iTQvj0rFSKiphb1\noXL8uARKmzYSKtXVkT4iIqLoFfWhYnR9AaxUiIiaWtSHijFIDzBUiIiaWtSHCisVIqLmE/WhEo5K\npb4+vMdERBStXBEqRqXSqZPzgfrqauCss4ANG8J/bERE0SaioZKXl4f09HSkpaVh4cKFjR7ftm0b\nxo0bh3bt2uHJJ5901NYQavfXxx/L/7y+PRGRvYiFSl1dHWbPno28vDwUFhZi6dKl2Lp1a4PndO/e\nHc8++yweeughx20NoXZ/FRUBHToAW7Y4a0dE5EYRC5WCggL0798fqampiI2NxdSpU7F8+fIGz+nZ\nsycyMzMRGxvruK3h2DEJBSC4UPn2W+D73wd273bWjojIjSIWKqWlpUhJSTlzOzk5GaWlpWFva6ym\nB4ILlT17gNGj5Tr3REQUWMRCxePxNEvbkyeBdu3k42BW1O/fD4wYAZSXO2tHRORGMZF64aSkJJSU\nlJy5XVJSguTk5LC3zcvLhccD5OYCF12UhZMnsxwd5/79QJ8+wKlTsuWL0ZVGRBQt8vPzkZ+fH5bP\nFbFQyczMRFFREYqLi5GYmIhly5Zh6dKlps9VSgXd9oILchEXB8ydCygla05OnwZiNL/y/fuBHj3k\nAl+HDzNUiCj6ZGVlISsr68ztefPmBf25IhYqMTExWLx4MSZOnIi6ujrMmDEDgwYNwvPPPw8AmDlz\nJioqKjB69GgcOXIEbdq0waJFi1BYWIhOnTqZtjXjO6bi8UhXWE2Ns1Dp3h045xwJld69w/HVExFF\nJ4/yLwOiiMfjwY9/rDByJDBzptzXvTuwfbv8b6emRq5tX1sLjB0LLFoEXHBB0x4zEVGkeTyeRj1E\nuqJ+Rb3vQD0gH584odf26FEJFcBbqRARkTVXhsrJk3ptGSpERM64LlTatw8uVLp2BQ4dCv/xERFF\nE9eFCisVIqKmw1AJ4MgRoEsX+ZihQkRkj6ESgH+lEkz3Fy8KRkRuEvWhcuJE4zGVYGZ/GYsfndi0\nSbbdLyx01o6IqLWK+lCJ5JjK6tXyP6/FQkRuwVAJwDdUgtnh+KuvgFGjgK+/dtaOiKi1YqgEcORI\naKHy7bdAdrZc6IuIyA1cESrG3l+A83UqxuyvYEPlggu4bT4RuYcrQiUc3V8dOsjW907s3QsMHw7s\n2+esHRFRaxX1oXL6NOB7NWIne3/5Xj/FaaVy8qS89rnnylTkujr9tkRErVXUh0q7drLlve9t3Uol\nlFA5cEB2Qj7rLJk5duCAflsiotYq6kPl7LMb3nYSKidONAyV48flQl86jFABgJ492QVGRO7gulBx\nMlB//Lh3kP+ss6QbraZGr+2BA0C3bvJxfDxDhYjcIepDpW3bhreDHVMB5GPdLjD/SmXvXr12RESt\nWdSHin+l0ratXMlRx4kTDacjOxlXYfcXEblR1IeKf6Vy9tn6XVj+lUqwoRIX53wzyhMngBUrnLUh\nIoo0hkoAvgP1gHewXkdVlTdUunSR1flOPPcccPXVsiklEVFrwVAJwHegHnBWqRw+LDsbA8FtRvnR\nR9J+zRpn7YiIIsl1oaI7pnLqlKxv8V046WSg3nc1fjCVypdfAj/8IbBli7N2RESRFPWh4j9Qr1up\n+FcpgLPur6NH5VoqgPNQqasD9uwBxo8HduzQb0dEFGlRHyrBdn/5D9IDzrq//K/F4iRUSkuBHj2A\n9HRg1y79dkREkea6UNHt/vKfTgwEf9VIp5VKeTmQlAT06gVUVuqv4iciijTXhUoolUqwOxw7DZW9\ne2VtS6dOMq5TXa3flogokhgqFvynEwOhVSpOZn/t2ydbuwBAQoJUK0RErUHUh4rZivpgB+pDrVR0\nu7GMSgVgqBBR6xL1oWJWqeiOqQRbqZw+LcFltG/bFoiJ0Q8k/0rF6b5h9fUMIiKKjIiGSl5eHtLT\n05GWloaFCxeaPue+++5DWloaRowYgY0bN565PzU1FcOHD0dGRgbGjBlj+RqhjKkEW6lUV8tMMd/r\nuDjpAtu3L7RK5dFHZZC/qspZOyKiUEUsVOrq6jB79mzk5eWhsLAQS5cuxdatWxs8Z+XKldixYweK\niorwwgsvYNasWWce83g8yM/Px8aNG1FQUGD5Omazv3RCxf8yxIB+peLb9WXo1El/wH3v3tDGVN55\nR9qtWuWsHRFRqCIWKgUFBejfvz9SU1MRGxuLqVOnYvny5Q2es2LFCkyfPh0AMHbsWBw6dAiVPu+w\nSmOQwj9UYmLkf7vL+5qFipNKxT9UnKxx8a1Uund3VnHU1QFffw386EfAZ5/ptyMiCoeIhUppaSlS\nUlLO3E5OTkZpaan2czweDy677DJkZmbixRdftHwd/4F64z67aqWmxvwCX8FWKk5C5eBB7wW+unaV\n27q++UaqlFGjgKIi/XZEROEQE6kX9vgOOARgVY18/PHHSExMxL59+5CdnY309HRcfPHFjZ730Ue5\nyM2Vj7OyspCVlXUmVPwH4n2FWql07Njwvk6d9EPl0CFZhQ9IqDjZNn/LFmDIEOC884Bvv9VvR0Tu\nlZ+fj/z8/LB8roiFSlJSEkpKSs7cLikpQXJycsDn7NmzB0lJSQCAxMREAEDPnj1x7bXXoqCgwDRU\nrrgiFw880PA+nXGVmprgx1ROnGgcKrqVilIy/bhLF7kdF+esUtm5E+jfH0hNBYqL9dsRkXsZf3Ab\n5s2bF/Tnilj3V2ZmJoqKilBcXIza2losW7YMOTk5DZ6Tk5ODV199FQCwfv16dO3aFQkJCTh+/DiO\nHj0KADh27BhWrVqFYcOGmb6O/5gKoDet+OTJxt1fupWK1WaUOqFy7Ji8rrE7stNKpaICSEyUMKqr\nc35xMCKiUESsUomJicHixYsxceJE1NXVYcaMGRg0aBCef/55AMDMmTNxxRVXYOXKlejfvz86duyI\nl19+GQBQUVGB6667DgBw+vRpTJs2DRMmTDB9HatQsatUTp70dkEZnFQqwW5Gefhww9d1WqmUlwPD\nhsl05tRU6QIzruuiQykJXLOxKCIiOxELFQCYNGkSJk2a1OC+mTNnNri9ePHiRu369u2LTZqXRDR7\ncwy2+yvUSkVnSrHveAoQXKXSu7d8nJIClJQAI0bot7/rLuD//k+mMfteS4aISIfrVtQDwXd/NVel\n4ltZdOokx3LqlH1bQCqVXr3k4/h4mZ6sq7YWeO01Gc/597/12xERGVwbKpGoVILp/vJ4nFUrvpVK\nfLyzLV6++EK6zKZOBdau1W9HRGRwZajodH9FslLxH8vRXatSUyNrZIw1Lj17OguVL78Ezj8fGDkS\n2LBBvx0RkcGVoaI7UB+JSsV/TAXQr1T27ZMrRrb57rvqtPvLmI6cliYfExE5FfWhYrWi3m5Mxar7\nq6bGfgv7UC5FbFapdO4sFYidqirZ1sXgtPvrm2+Avn3l3zff8IqTRORc1IdKKJWKfyB5PHpdZ2aX\nItYNlaNHvQsfDU5Cxej6ApyHys6dQL9+Emrt2jnfcj8/H1i50lkbIoouEZ1S3ByCHVMxq1QA77iK\n2WOGUCqV6mrpwvIVbKj07Om8+6tfP/m4Xz/vPmI6jh8HJk2Sa8lUVDSsmIjIPVxbqQQzpRjQG1ex\nqlR01qlUV8s0Yl9dugQXKj166IdKdbUctxFoffoAu3bptQWAf/4TGDcOuPJKIC9Pvx0RRZeoD5Vg\ndyk2G6gH9GaAhVqp+IdKsJVK+/YyLqIzY82Yimzs85mYCJSV2bczrF0LfP/7wEUXAQEub0NEUS7q\nQyXc3V/t2tm/SZtVKrq7FIczVDweua0zHbmiwrtoEpCAKS+3b2f46itZuT9ihKx3ISJ3cmWohNL9\n1b69ffdXS6lUAP29w3wXTQLOQ2XbNiA93RsqTmeOff21HAMRtW6uDZXmrlSchIr/tvmdO8t2+HYO\nHAg+VHy3dwGchUpNDbBnjwzuJyTI+XXSdWZsgnnRRTLQT0StlytDJdgV9UDwlUr79vKadpcxDnWg\nPi6u4X1OKhXfUHEyprJjh1wUzNiA0ljnomvpUuDWW2WRJ/ccI2rdoj5Uwj1QH2yl4vFI0Bw/Hrht\nKN1f/ptRAs3T/WWsxDc4nTn20UfAxInAhAnAhx/qtyOilifqQ8Vs+3a7MZX6eumGMatygq1UAL0u\nsFBC5ciRxqvxdUNl715ZLGk45xw5R3YhCMj2+ikp3ttOQkUp4JNPgAsuAH7wA+Bf/9Jr56u62r4C\nJKLmEfWhctZZje+z6/6qqZHnGNNrfdlVKvX11uMxdqGilPWYim6l4r8aXzdUDhxouGDR45HFk/v3\n27cNJVSKi+V7lJICjB4NbNok51DX11/L13j99fptiKjpRH2omLHr/rIKBcC+UjHGYtqYnFm7UKmp\nkTdY/wqpOSoVs5ljPXrohcqePcGHSmEhMHy4d4v/uDi5WqWuJ58Efv5zCaONG/XbEVHTcG2oBOr+\nshqkB+wrFbPxFIPdmIpZ1xcg1Yfd7C8j6PzDMNRQ0VmRX1ICJCd7bycnA6Wl9u0AqTQGDvTeHjIE\n2LJFr61SwDvvALfdBlx3nXxMRJHl2lAJVKlYDdID9pWK1XgKYF+pWIWKUakEWvtx5Ejjri9AL1SU\nkuf4h0qw3V/GzDGdtSpffw0MGOC9PXSoLKTUUVQExMTIbLOJE4H339dr53vcOTly+WQiCg9XhorO\nmIpVqESiUmnbVrrTAh2z2Zb5gIRKVZV1O0ACqX37xpMadCoVpSRAEhO993XqJJ9L5xow27c3rFTS\n0mSKso41a2RrGI8HuPhiubCY3aw+Xw89JOdn5ky9ao6I7LkyVHQqFavur0hUKoD9uIrZID2gt02L\n2aJJQK9SOXJEAsR/coHuOhf/6chO1rhs2gRkZsrHHTrI59GtcqqqZOPL3/4WuOwy4LXX9Nr50pkZ\nR+Q2rg2VQGMqLa1SAexDxWyQHtDr/vK/uJdBp1KprDTfHj8pyX5cpa5O1sckJXnv69tXf5D/yy9l\nkN8wciTw+ed6bd9+Gxg/Xs7ZrbcCr7+u186waJEE6U9/6qwdUbRzZajYdX+1xErFblW9VaWiGypm\nlYrO7C+rUNGpVMrL5TV8u93OPVfanToVuK1SEirDhnnvGzlSfwbYmjVSoQBAVhbwn//o7eYMyPnK\nzZVK6dVX9asjIjdwZaiEMlDf2iqVdu3kDThQEIbS/RVKpeI/wA9IwCQmArt3B267Z4+cT98Lmg0Z\nIhtb6vjoI+CSS+Tjzp2l4lm3Tq/tkiVy3ZgRI4DZs4HnntNrZ1iwQGbI/fnPztoRtQauDZVgu7+a\nslI5dixwqASaVmxVqXg89tVKU3R/6VQqZqEC6I2rFBU1nDUGyC7JW7cGbgdIUO7bBwwe7L3ve98D\n1q+3bwsAf/+7dJkBwM03y+wx3Y0wCwpkHGfJEuCBB+TrcOL4cZly7XQXaKLm4spQCaX7K5RKpWPH\n5q9UAFlUGGgmllX3VyiVik6o+C+aNPTpYx8qO3Y0HOA3XvPYMfvuvs2bpdvMd4FqZibw2WeB2wFy\nrr78UmadAbIzc69e0n2mY/584JFHZDzn7ruBJ57QawdIEJ5/vmxnM2sWg4VaJleGSlOuqA9UqXTo\n0PyzvwD7SsV/ixZDt27yJhpo25RQu798F00adCoVs1DxeKRa+frrwG2NUPGlGypr1khV4/vzMX68\nXE7Zzt69QH4+cMstcvvuu2WCgO4ssvvvByZPlhlzH30EvPmmXjtA/ui4804Jw2D2VyPSxVAx0RIr\nFbuB+kCVSlxccJVKbKyEWaBAsgoVnV2Orbq/UlPtt2rZsUPWtPhLT7cfV/GfNQZIdXTsmP2Fwtav\nBy68sOF9uqHy3nsyOaBzZ7mdkCD7na1cad92505g1SqZINCpE/DUU/KxTrWilHTX1dQA99wD3HCD\nnANdlZUSgLfcIoFMFEhEQyUvLw/p6elIS0vDwoULTZ9z3333IS0tDSNGjMBGn6k9Om2ttG1rv01L\noIH61lapdO1qX6mYhQogXWCBxlWsQqVXL3ksUJUTKFSKi63bAeaVCqAXKmaVisejV62sXw+MHdvw\nvksuke4vu4rj/fdl5b+vG24Ali0L3A6Qa87ceKM3kCZOlHGc/Hz7tmvWyFjTH/8or/fYYzLBQCeQ\nTpyQr69dO2DUKAnQTZvs2xnee0/aZ2UBq1frtzPs3au37x21HBELlbq6OsyePRt5eXkoLCzE0qVL\nsdVvlHXlypXYsWMHioqK8MILL2DWrFnabQMJtfsrUmMqgQbqQ61UzLq/APtpxVahcvbZEnKB2gYb\nKvX18pd7v36NH7MLlfp6GegeOrTxY6NGyap8K6dPy+NjxjS8v3NnCalPPw38uqtXyzVjfF17rVQg\ngf5QUQr429+Am27y3ufxAPfdBzz7rHU7w1NPyaQAY6PSGTPkj5Dly+3bPvmkVHVPPSXdb08/DUyb\nZn/5BwB44w15rfvvB+69V/Zo0wlQQMbbLrtMdltITJQKS3e3BKWk+vvpT4H//V8JJifq6qSN3SXH\nyZxtqBQWFja6L1/nzyMbBQUF6N+/P1JTUxEbG4upU6diud9P+YoVKzB9+nQAwNixY3Ho0CFUVFRo\ntQ0k1O6vllqpBBqot5v9Fe5KBQjcBXbqlASO74XBDL16SQhaBXB5uZwP4692X3ahsmuXfK3+FzMD\nZJ1LoFD56itZR2N2ni+5JPBVKzdulIA+99yG9/foIVOTA3Wfffml/LEyblzD+2++WS5qVllp3Xb7\ndrlezW23ee876yzgV78C5s0LXK2Ul0uI+HYE3HyznOPHHrNuB8istrvukjf3a6+VSxO8/75USHaV\nzsGDQHa2nNN9+2R6eWUlcPnl8jsSyOHDsp/bnDnyx9S2bTLL79VXA7cDJEyeeELG+QYPlp+T6dPt\np7cDch7z8oA77pAwnD4deOst/ev8fPutzAh86ikZZ7PbWslXfb3s9v3pp/p77vm3D2eA2obKDTfc\ngIULF0IphePHj+Pee+/F3LlzQ37h0tJSpPj8mZqcnIxSv5Fdq+eUlZXZtg0kJkZOpNU3PJSBertK\nxS5U/Lc7MYQ6UB+oUgnU/RVoWnF1tfwAWwVhoFApK5MwMrveTZs28uZrNa5iNZ4CSJdYcbH1L4nZ\neIrBLlTWr5eLiZm55BIZPLdi1vVluPpqYMUK67Z/+xswdWrj6/t06SJv2IHWuyxaJHub+f+hM3my\nfO/eftu67cMPS6XRt6/3Po9HpkT/4Q/WEyLq66Xdww/LbDXD8OHAM8/I12L1e1BbC0yZIgHyi1/I\n72pcnLzR9usHXHWVdduqKnlDT06WnRV+/nPp8vvXvyQEf/IT66nfFRXymu+8I0G9f7+ESWqqdIv+\n5S/Wb9abN8uMvIceAjIypEK64AKZ6TdqlHzvrdquWydfU2amhNLu3fL97NNHfi7eecf8fUop+Vn9\n2c/kuZMnSzU3fLj8bjz0kPy8mnU/19bKOZkzR74/xjZLPXrIOXjySfNj1WUbKp9++ilKSkowbtw4\njBkzBr1798Y63VViAXjMroBlQoU4bzI3N/fMP6PC8ngCVyuhDNTbVSqRGKgPVKnU10vgBKpUrLqw\njCrF6lsZKFSsur4MgbrAduww7/oC5PuWkmI9e8xsPMXQt6+cR6sQNRtPMVx0kfylaBVmgUIlJ0dC\nxewNoL5e9iXz7fryNWMG8NJL5m9aVVUSSPfc0/gxj0fetH/5S/O2mzYB774L/M//NH4sKUnutxqX\n+d3v5I3w3nsbPzZtmlRcd93VuK1Scn+nTtJt5atNG+CFF+TnIien8e9RRYWM21x6KfD73ze8JtGw\nYfK927YNuOKKxr8Lq1bJm/+FF0rFOGiQ3N+tm1Rz778vAXHttQ3/0Dl4UAJk/Hjgv/5Lztm998r3\nedYsec1HH5Uwu/RSCYgjR+T3bflyCaJp0yRUSkpk3OyZZyToS0slVH79a/maH31Uvh9//7sExsCB\n8pqxsfJ5d+yQ8cB9+6TbsUMHqZrOO09+RubNkwC6/HL5nZ47V35Xfv97eR9ctSofP/xhLs45JxdL\nl+Y2/sY5EGP7hJgYtG/fHidOnMDJkyfRt29ftDG7ApVDSUlJKCkpOXO7pKQEyX7zS/2fs2fPHiQn\nJ+PUqVO2bQ25ubmm9xtrVcwCIJSt71tqpWIVKkeOyGvGWPwk9Owp/dtmAnV9AYHXqoQaKlaVCuDt\nAktPb/zY5s3y5mDG4/FWK2YB8OmnMjZhpmtXqZI2bGhczRw9KvcbK/j9paXJ9+izzxqP16xb5x2z\nMXPRRfJmvG6dfOzrhRfkDdisixEArrlGZpCtXCk7BBiUkq/z0Uet/1C5917glVekgrjxRu/9X38t\nb2Br15pXoYCEzrhxMh50333e+x99FPjiC5lYYFXBvvSSdC1dfTXwj3/Iufn8c6l+pk+XsDP7Iycu\nTt6Uf/YzqQoefFC6WV9/Xc7dq69KOJjJyJDv329+I3/ZDxkib8gbNsjP0ldfNbwUt8HjkccnT5ZK\nZ+FCOVdt2sgx/PCH0p1odsnzTp0kFO64Q6rrV14BFi+W1x05EvjrX+Vz+H+tHo8c4/nnyx8MW7ZI\nt2xpqZyDu++Wtv5jqJdemoVLL83y+TzzzE+GDmVj+PDh6uGHH1a1tbWqrKxMTZ48WU2ZMsWuma1T\np06pvn37ql27dqmamho1YsQIVVhY2OA57777rpo0aZJSSqlPPvlEjR07VrutkhLH8vV79lSqosL8\nsR//WKnnnjN/7PRppTweperrzR+/4gql3n7b/LG9e5Xq0cPykFRqqlI7d5o/9tlnSmVkmD9WV6dU\nmzZybGY++ECpSy81f2zHDnldK6+8otS0aeaP/eMfSuXkWLddtEipe+4xf2zhQqUefNC67WOPKfWz\nn5k/NmWKUkuXWrd96CGl5s83f2zgQKU2b7Zu++CDSj3+eOP7q6qU6tzZ+hwrpdR998nX5W/5cqXG\nj7dup5RSc+Yo9T//0/j+u+9W6te/Dtz2N79R6vbbG95XU6NUUpJSGzcGbvv660qNGdPw5/mNN5Qa\nPFipU6cCt/34Y6USE5WqrJTb1dVKjRyp1O9+F7idUvJzl5ys1COPKLVpk3ydgwdb/076OnVKqZkz\nlYqLU+r88+V3+a9/tW9n+OADpW65RamrrlLqqaeUOnxYv211tVIffqhUXp5S+/bptzPU11u/d7Qk\nGtFg3dbuCQUFBY3uW7JkSdAv6GvlypVqwIABql+/furx736Tn3vuOfWczzv6Pffco/r166eGDx+u\nPv/884Bt/QU6McnJSn37rfljt92m1MsvWx93bKxSJ0+aP3bppfJDa6a6Wqn27a0/b48e3l9Qf19/\nrVT//uaPHT6sVMeO1p/388/ll89MQYFSo0ZZt333XaUmTDB/7A9/UOrOO63bvv66UtddZ/7Y7NlK\nPfOMddu//lWpG24wfywjQ6n//Me67R//KN9Df8ePK9WunVK1tYFf1+xvprw8pbKyrNspJW/GV17Z\n+P5Zs+SNP5C1a5UaOrThfSdPys+E1R8ahooKpbp2bfjm+NxzSk2cGLidUvIHyeDBSv3973J7/36l\neveWwNDxyCNKpafLz8KYMRJuum+apaVK3XijBP1ddyl14IBeO8OePUp9+qlSJ044a0f2mjRUWrNA\nJ6ZfP6W2bzd/7MYblfrb36w/b5cuSh08aP7Y2LFKrVtn/lhdnVQ5dXXmj7dvL8FjprxcqYQE88dK\nSuQvRitpq0EGAAAZPklEQVQ7d1pXI++9p1R2tnXbggL569NMbq5SDz9s3fbjj5W64ALzx665Rt6E\nraxdK29S/urrpWKoqgr8umZtP/tMqeHDrdsppdTWrUr16dP4/kceUernPw/ctrJS3tx9q5n6eqXO\nPVepLVsCtz19Wqn4eKWKirz3vf66dYXpb+pUpX71K/n4+HGlzjvP+ufQ37p1El4LF8ofH//v/+m1\nU0q+vmXLlJo+XamXXrL+2abWJZRQceWKesB+oN5qTAUIPK040JhKmzbWA/11dfI5rdoGGlMJNJ4C\nBB5TCTSdGAg8pdhuTKUpBur375c+6Lg467bGmIr/QHCgmV+GAQPk6/U/Xx9/3HjMwl98vAxi++4D\n9tVXMj5gDP5aOessWfX+wgve+154QQZZdTz+uMz0evNNaXPRRY2nIFsZN07GVb75RsZK7KYL+/J4\nZEHlK69I/38YhluplXPtj0CgUAk0pRgIPFgfaPYXYL0A8tgxeczql7JDBzkus+mFgWZ+AfJYdbX5\n7KJACx8Bb6iYzfLRCZWKCvO2dqHSq5eEpf+5CjTzy9C9u3x//bdc0QmVNm1kkNN3avGpUxIUOm/S\nV1/dcFHhihUyCK4z2fHuu4E//UkmN/zrX/K1Tpli3w6QaaWvvSazhc46y/l2/KNHS5s77tA7ViIr\nDBUTgaYUA4GnFQeqVADrBZBGqFjxeORxs2rFrlJp00Zmkxw+3PixQGtUAHlNj8f8mO1CpX17+ef/\nV/+JE3IsgdparVWx2p7FX3q6LAjzpRMqgMzeWrvWe/uLL6RyClQdGa6/Xt7c6+okxF95RaaN6ujb\nVyqF8eNlCrEx20fX+PEyE+rPfzZfGErUHBgqJkLp/gq2Ugm0RsVg1QVmV6kA1gsg7bq/AOsFkHah\nAki14j+teM8eWaBm11Vi1gVmtT2LvxEjJAwMSsltnVDx3yDy/fdlTYGOkSOlynrjDemK6tjRem2L\nmV/8QqaevvFGw2m+RK2Fa0Ml0KaSNTWB/0IMtP9XsJWKTqhYLYAMtEWLwWoBpF33F2C9AFI3VPzH\nVXbvDtz1ZTALFd1KxX8fr5ISWYvTq5d924svlm1VjBB++21Za6DriSdkweFdd8nqcyfdSR6PrC35\n3vf02xC1JK4NlaaoVOrr7cdj7MZUArGqVOy6vwDrwXq77i/AfLD+xAkJZbswS0w0DxX/PbDMmF2s\nq6hIL1RGjpSuIMO6dbJiWucNvmNHWfy4bJns7rtrl/XCRTPf+55sHvn++87aEUUD2xX10cpuoN6u\nUjELFWMsJlC3jtWqet3uL7OdinUrFbPur/37G17n3YxZqFRWymwnuzfpxMTGF+vSDZVBg2TvJoOx\ncZ7vZYADtS0r8359RqjomjlTttpYvlwG0H23/dCRkeHs+UTRgpWKCZ1Kxaz7y248BQit+ytQpaIz\npmJVqeh0f5mFil3XFyB7D/kPtuuGyuDBDQfbd++Wr1NnwDw2VqqEf/5TxlPef1/2htJ12WWyjUZC\nguzvRER6GComgp1SbDeeAgTu/gplTMVsK3dfVpXKgQPBVyo6oWI2LqIbKn37SteZEcJffWV+HRQr\nEybIlN4tW+R7M2qUfltA9pF6+WX7PxSIyIuhYiLYKcWhViqhjKkEU6mcOiWva9fWbKDeSaXiHyrf\nfqsXKjExshjRqFY2bbLeWNHMLbfIduL33y/Tern+gqjpMVRM2FUqVgP1oVQqTT2mYjaluKpK7reb\n2tuzZ+Or5zkJld27vQsgT5+W23362LcFZMGhcaWFtWudjYvExcnsq/h42eqbiJqea0PF2PreX12d\n/LPaCh6wnlIcSqWi0/1lVakcOhTclGKdri/AfFqwbqh06iRfs3F1wuJi+XyBQtvX974nW6TU1ckV\nDJ1OtZ02Tbb6tjs/RBQerg2Vs882X6diVCmBukqaqlJp6u4v/0pl/377QXpA9rPyX8BYWiozu3QM\nGOC9SuD27XJb1yWXAPn5clXFlBSpmoio5XJ1qJhVKnbjKYD1QH1Tz/4K9+JH3UolIUEG6n0vxbpn\nj94CRkAuamSMi2zfHvgCW/7OPVdWwV9zjexLRUQtG0PFj914ChDaQH0os7/MKhWl9Ldp8Q8V3Uol\nJkbCx+jCAmSFusXFNhvxnRq8ebPeOhNfr7wiW5fMnu2sHRE1P4aKn6auVAItftTp/vIfqD9+XNZk\n2C3OM5tSrLNGxZCU5F3EWFsrba0uU+tvyBCZDgzIJXn9L7drJyVFtjwxu+wqEbUsDBU/dgsfAetK\n5cQJve6vcG4oqdP1BchWLFVVDbeh1+3+Ahpeb76sTPbQsroGub8xY+T66xUVMlDvZFowEbUuDBU/\ndlu0ANYD9aFWKsGMqeiGytlnS4Xlu/29bvcX0HAbeiddX4BUSSNHAg8/LFUKKw6i6OXavb+sphTr\nVCqBphTbzf6yqlR0NpTs0qXxNVF0QwWQAffKSu/qe519vwxpabKZI6C/qaOv6dPlioRvvOGsHRG1\nLq6uVKymFLfUSsVsXMRJqMTHN1zEWF6uPy7iGyrbtslFsJy44w7p/rr+emftiKh1cXWohFKpWK1T\nacoxlY4dJQh9w1Bn3y+Df6hUVOhdXwRoHCp211w3o7NYkohaN4aKn+aYUhzM5YQBWZDpX63orKY3\nGN1fgGwjv3evfqikpsoA/cmTwVUqROQODBU/OlOKA3V/6Yyp+IfK6dNSfdi1BRqvNzl0KLhKZf9+\nGaPRvU5IbKxUJ/n50m3mdEyFiNyBoeJHp1IJde+v48cbTu01qhSdXXT9KxUna018KxUn4ymGrCzg\ngQdkk0fO4CIiMwwVP6FWKnahEhMj/3xfW6fry+BfqThdwLhnj3wcTKj86Edyed377nPWjojcw7VT\nikOtVIIdqAe8W7UYr6MzSG/w38Orqsr+GvOG1FS53jog4ZKUpNfOMGSI7BZst1U+EbmXa98e2rY1\nn1KsW6kE2/0FNB6sP3pUP1T8dxt2UqkYV2FUSmZyOdnY0cBAIaJAXPsW0RSVis5APSABUl3tva2z\nIaTBv1JxEirnnCNhun9/8KFCRBRIREKlqqoK2dnZGDBgACZMmIBDZhdPB5CXl4f09HSkpaVh4cKF\nZ+7Pzc1FcnIyMjIykJGRgby8PMfHEMqYSkyMTMn13Qoe0K9UunRpuDHkkSNyn464OOnyMjgJFUCu\nuLhrl/PrmhAR6YhIqCxYsADZ2dnYvn07xo8fjwULFjR6Tl1dHWbPno28vDwUFhZi6dKl2Lp1KwDA\n4/HggQcewMaNG7Fx40Zcfvnljo8hlMWPHo95F1iwoeJ0Vfy+ffKxUhIqumMqgPea7998w2nBRBR+\nEQmVFStWYPr06QCA6dOn46233mr0nIKCAvTv3x+pqamIjY3F1KlTsXz58jOPK985uUEIZUNJQMZF\n/ENFd6A+lEolPt47Lbi6Wqb26l6aFwBGjwaef16uHa8744yISFdEQqWyshIJ3+3ZkZCQgErfqz99\np7S0FCk+lxZMTk5GqXFBDwDPPvssRowYgRkzZlh2nwXSti1w6pR0Y/nSqVSAxosYldIfUzGrVHRD\nxXetyb59+htCGiZPBtavB3JynLUjItLRZFOKs7OzUVFR0ej+xx57rMFtj8cDj8mqP7P7DLNmzcIv\nfvELAMAjjzyCBx98EC+99JLpc3Nzc898nJWVhaysrO8+v3cGmG+IOKlUfPfwqq2V64vEaJxRs0ol\nLs6+HdAwVEpLnU8LTkuTa5twmxUiMuTn5yM/Pz8sn6vJQmX16tWWjyUkJKCiogK9evVCeXk54uPj\nGz0nKSkJJSUlZ26XlJQg+buLePg+/0c/+hEmT55s+Vq+oeLP2P7eN1SCrVR0x1MA81A57zy9tsaY\nSn297MXlNFQAYNQo522IKHr5/sENAPPmzQv6c0Wk+ysnJwdLliwBACxZsgTXXHNNo+dkZmaiqKgI\nxcXFqK2txbJly5DzXZ9NeXn5mee9+eabGBbkpQTNtr/XmVIMNN5tWHc8BQhtoP7ss2VK8sGDwVUq\nRERNKSKhMnfuXKxevRoDBgzAhx9+iLlz5wIAysrKcOWVVwIAYmJisHjxYkycOBGDBw/GjTfeiEHf\n7bc+Z84cDB8+HCNGjMCaNWvw9NNPB3UcZoP1OlOKgcbdX6FWKrpjKoC3C4yhQkQtTUS2aenWrRs+\n+OCDRvcnJibi3XffPXN70qRJmDRpUqPnvfrqq2E5DrNQcVKp+Hd/6QzSA6EN1AMSJCUlEirsyiKi\nlsS1K+qByFYqvpcFdrKiHgD69pUFjHv2OLtWPBFRU2OotIBKxcmFtgCgXz/ZZmXrVmDgQP12RERN\njaESpkqluhro3Fnvdf3379q/39l6k/PPB1askM0dTSbOERFFjKtDxZhS7CvY2V9Otq/v0UO2VwFk\n1tjp0/ptAeCCC4AdO4CLL9a7sBcRUXNx7fVUAPMpxbqVSocODasNJ5VK9+4SKsbeXT16OAuHrl2B\nVauAoUP12xARNQfXh4r/Fva6ix87dpTZVwYn10Rp21ZC6fBh6fpyssuwITvbeRsioqbm6u4vs+ui\n6G7TEkr3FyDVyf793kqFiCgauDpUzK4176RS8Z395aT7C/CGitNBeiKilszVodK+fcPt60+flnEO\nnU0h/SsVJ91fQMNQCab7i4ioJWKo+ITKiRNyn86geadOEiSGYCsVp1duJCJqyVwfKr7dX043hfQP\nFaeVyr59stNwYqJ+OyKilszVoeJ/SWAnq+LPOafhVitOu79SUmT/rt27gXPP1W9HRNSSuTpUrLq/\ndPhvteK0+ys1FSgullDxucAlEVGrxlDxC5Vgt6932v3Vrx+wZQvwzTdyNUYiomjg+lDxHVNx0v3V\nrp3MFjNW5B896qxSGTRIAqVDB/0gIyJq6Vy9oj6USsXj8Q7Wd+/ufEwlJga4/37Zxp6IKFq4OlT8\nB+qdjKkA3uuidO0qoeJk+3oAeOopZ88nImrpXN/9FezsL8A7rmJcufGss8J/jERErYnrQyXYdSqA\nN1SqqoBu3cJ/fERErY2ru7/CVam0a8dQISICGCpBD9QD3lCJiWGoEBEBLg+VUAfq4+K8V3BkqBAR\nuTxUzNapOFlr0rs3UFEhHzNUiIg4UB9S91evXkB5OQfqiYgMDJUQur9692aoEBH5cn2oHD8uF+YC\n5GMnlYrR/VVRAcTHN80xEhG1Jq4eU4mNlZlbNTXeQftgKpX27YHzzmu64yQiai1cHSqA7NdVXS2h\n4nRTyPh47+wvbl9PRBSh7q+qqipkZ2djwIABmDBhAg4dOmT6vDvuuAMJCQkYNmxYUO11GKECOA+V\nmBjg7LOBykogOTnoQyAiihoRCZUFCxYgOzsb27dvx/jx47FgwQLT591+++3Iy8sLur0O32vNOw0V\nAJg2Dbj6aqCNq0eniIiERyljmLr5pKenY82aNUhISEBFRQWysrKwbds20+cWFxdj8uTJ2Lx5s+P2\nHo8Hdl/euHGyW/C4cXI1xvx8+Z+IyK103jutROTv68rKSiQkJAAAEhISUFlZ2aztfflWKkeOOK9U\niIjIq8kG6rOzs1FhLDf38dhjjzW47fF44PF4gn6dUNsbYypKBdf9RUREXk0WKqtXr7Z8zOi26tWr\nF8rLyxHvcJGHk/a5ublnPs7KykJWVlaDxzt3llCpqZFxkbZtHR0KEVGrl5+fj/z8/LB8rohMKc7J\nycGSJUswZ84cLFmyBNdcc02TtfcNFTNG9xerFCJyK/8/uOfNmxf054rImMrcuXOxevVqDBgwAB9+\n+CHmzp0LACgrK8OVV1555nk33XQTLrzwQmzfvh0pKSl4+eWXA7YPhtH9xVAhIgpdRGZ/NRedGQy/\n+hVQWwtMmQLceivw5ZfNdHBERC1Uq5v91ZKw+4uIKHwYKuz+IiIKG9eHSpcuwOHDDBUionBwfah0\n7y6bQh44wGuiEBGFiqHyXajs3w/07BnpoyEiat0YKj6h0qNHpI+GiKh1c32o9OghgbJvH0OFiChU\nrg+VDh3k3xdfAImJkT4aIqLWzfWhAshW94WFQN++kT4SIqLWjaECb7cXr95IRBQa11+jHgDuvBM4\n5xy5PDAREQXP9Xt/ERFRQ9z7i4iIWgSGChERhQ1DhYiIwoahQkREYcNQISKisGGoEBFR2DBUiIgo\nbBgqREQUNgwVIiIKG4YKERGFDUOFiIjChqFCRERhw1AhIqKwYagQEVHYMFSIiChsGCpERBQ2EQmV\nqqoqZGdnY8CAAZgwYQIOHTpk+rw77rgDCQkJGDZsWIP7c3NzkZycjIyMDGRkZCAvL685DpuIiGxE\nJFQWLFiA7OxsbN++HePHj8eCBQtMn3f77bebBobH48EDDzyAjRs3YuPGjbj88sub+pBbvfz8/Egf\nQovBc+HFc+HFcxEeEQmVFStWYPr06QCA6dOn46233jJ93sUXX4y4uDjTx3iZYGf4C+PFc+HFc+HF\ncxEeEQmVyspKJCQkAAASEhJQWVnp+HM8++yzGDFiBGbMmGHZfUZERM2ryUIlOzsbw4YNa/RvxYoV\nDZ7n8Xjg8Xgcfe5Zs2Zh165d2LRpE3r37o0HH3wwnIdORETBUhEwcOBAVV5erpRSqqysTA0cONDy\nubt27VJDhw4N6nEA/Md//Md//BfEv2DFIAJycnKwZMkSzJkzB0uWLME111zjqH15eTl69+4NAHjz\nzTcbzQ4zKI67EBE1K4+KwDtvVVUVbrjhBuzevRupqal4/fXX0bVrV5SVleHOO+/Eu+++CwC46aab\nsGbNGhw4cADx8fH45S9/idtvvx233XYbNm3aBI/Hgz59+uD5558/M0ZDRESRE5FQISKi6BS1K+rz\n8vKQnp6OtLQ0LFy4MNKH02xKSkpw6aWXYsiQIRg6dCh++9vfAtBfcBqN6urqkJGRgcmTJwNw77k4\ndOgQpkyZgkGDBmHw4MH49NNPXXsu5s+fjyFDhmDYsGG4+eabUVNT45pzYbaoPNDXPn/+fKSlpSE9\nPR2rVq2y/fxRGSp1dXWYPXs28vLyUFhYiKVLl2Lr1q2RPqxmERsbi6effhpbtmzB+vXr8bvf/Q5b\nt27VXnAajRYtWoTBgwefmWXo1nPxk5/8BFdccQW2bt2KL7/8Eunp6a48F8XFxXjxxRexYcMGbN68\nGXV1dXjttddccy7MFpVbfe2FhYVYtmwZCgsLkZeXh7vvvhv19fWBXyDoIf4WbN26dWrixIlnbs+f\nP1/Nnz8/gkcUOVdffbVavXq1GjhwoKqoqFBKKVVeXh5wxl00KSkpUePHj1cffvihuuqqq5RSypXn\n4tChQ6pPnz6N7nfjuThw4IAaMGCAqqqqUqdOnVJXXXWVWrVqlavOhf+sWauv/fHHH1cLFiw487yJ\nEyeqTz75JODnjspKpbS0FCkpKWduJycno7S0NIJHFBnFxcXYuHEjxo4dG5YFp63R/fffjyeeeAJt\n2nh/1N14Lnbt2oWePXvi9ttvx8iRI3HnnXfi2LFjrjwX3bp1w4MPPohzzz0XiYmJ6Nq1K7Kzs115\nLgxWX3tZWRmSk5PPPE/nvTQqQ8XpYspoVF1djeuvvx6LFi1C586dGzwWzILT1uidd95BfHw8MjIy\nLKeXu+VcnD59Ghs2bMDdd9+NDRs2oGPHjo26d9xyLnbu3IlnnnkGxcXFKCsrQ3V1Nf7yl780eI5b\nzoUZu6/d7rxEZagkJSWhpKTkzO2SkpIGaRvtTp06heuvvx633nrrmTVACQkJqKioACDrfOLj4yN5\niM1i3bp1WLFiBfr06YObbroJH374IW699VZXnovk5GQkJydj9OjRAIApU6Zgw4YN6NWrl+vOxWef\nfYYLL7wQ3bt3R0xMDK677jp88sknrjwXBqvfCf/30j179iApKSng54rKUMnMzERRURGKi4tRW1uL\nZcuWIScnJ9KH1SyUUpgxYwYGDx6M//7v/z5zv7HgFEBQC05bo8cffxwlJSXYtWsXXnvtNfzgBz/A\nn//8Z1eei169eiElJQXbt28HAHzwwQcYMmQIJk+e7LpzkZ6ejvXr1+PEiRNQSuGDDz7A4MGDXXku\nDFa/Ezk5OXjttddQW1uLXbt2oaioCGPGjAn8ycI9ANRSrFy5Ug0YMED169dPPf7445E+nGbz73//\nW3k8HjVixAh1/vnnq/PPP1+999576sCBA2r8+PEqLS1NZWdnq4MHD0b6UJtVfn6+mjx5slJKufZc\nbNq0SWVmZqrhw4era6+9Vh06dMi152LhwoVq8ODBaujQoeq2225TtbW1rjkXU6dOVb1791axsbEq\nOTlZ/elPfwr4tT/22GOqX79+auDAgSovL8/283PxIxERhU1Udn8REVFkMFSIiChsGCpERBQ2DBUi\nIgobhgoREYUNQ4WIiMKGoULUzA4fPow//OEPkT4MoibBUCFqZgcPHsTvf//7SB8GUZNgqBA1s7lz\n52Lnzp3IyMjAnDlzIn04RGHFFfVEzezbb7/FVVddhc2bN0f6UIjCjpUKUTPj33EUzRgqREQUNgwV\nombWuXNnHD16NNKHQdQkGCpEzax79+646KKLMGzYMA7UU9ThQD0REYUNKxUiIgobhgoREYUNQ4WI\niMKGoUJERGHDUCEiorBhqBARUdgwVIiIKGwYKkREFDb/Hxu/IkiMyOjQAAAAAElFTkSuQmCC\n", "text": [ "" ] } ], "prompt_number": 1 }, { "cell_type": "markdown", "metadata": {}, "source": [ "For a second example, suppose that you want to solve the following\n", "coupled, second-order differential equations,\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", "\\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" ] }, { "cell_type": "code", "collapsed": false, "input": [ "from scipy.integrate import odeint\n", "from pylab import * # for plotting commands\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 array([z[1], a*z[2], z[3], b+c*z[1]])\n", "\n", "time = linspace(0.0,10.0,1000)\n", "zinit = array([0,0.2,0,0.1]) # initial values\n", "z = odeint(deriv,zinit,time)\n", "\n", "figure()\n", "plot(time,z[:,0],label='x') # z[:,0] is the first column of z\n", "plot(time,z[:,2],label='y') # z[:,2] is the third column of z\n", "xlabel('t')\n", "legend(loc='upper left')\n", "show()" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEPCAYAAACneLThAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlcVXX+x/EXIBIqIqKgAoqmsihuuG9RhmXmUlojlpqp\nLWal+Zhstkf1+/1KGafUNnVKHW3RMjOaJsnRQs0E3ENFNBEXENxA2RS5nN8fZ3IZzWQ9l8v7+Xic\nx12895zPveqbL9/zPd+vk2EYBiIi4pCcrS5AREQqj0JeRMSBKeRFRByYQl5ExIEp5EVEHJhCXkTE\ngZU75B9//HF8fX0JCwu7/NzZs2eJjIykbdu2DBw4kJycnPIeRkREyqDcIT9+/HhiY2OveW7WrFlE\nRkZy4MABBgwYwKxZs8p7GBERKQOnirgYKi0tjSFDhpCUlARAcHAwGzZswNfXl8zMTCIiIti/f3+5\nixURkdKplD75rKwsfH19AfD19SUrK6syDiMiIr+h0k+8Ojk54eTkVNmHERGRG6hVGTv9pZumSZMm\nnDhxAh8fn+teo+AXESmb0vSyV0pLfujQoSxduhSApUuXMnz48Bu+zjAMbYbByy+/bHkN9rLpu9B3\noe/ixlv8sXhGrxpd6jwud8hHRUXRu3dvUlJSCAgIYMmSJbz00kv8+9//pm3btnz33Xe89NJL5T2M\niEiNU2Qr4pOkT+jxQQ+iVkXRtWnXUu+j3N01y5cvv+Hz69atK++uRURqpJP5J1m4bSELti8guFEw\nf+r3Jwa3GYyLswsv8EKp9lUpffJSOhEREVaXYDf0XVyh7+KKmvJd7Dyxk3kJ84hJieGh0IeIfSSW\nMN+w337jTVTIOPkyHdjJCYsOLSJiN2wlNmJSYpgbP5fDOYd5ptszTOoyCe863jd8fWmz0+5a8o40\n6kY/xETk1+QV5bFk5xLmJszFt64v03pO44GQB6jlXLGxbHchD44Rjo70w0pEKk76+XTeTnybD3Z8\nQERgBB898BG9AnpV2vHsMuRFRBzNrsxdvLnlTb4+8DVjOowhcVIirbxaVfpx7a5P3lH66h3lc4hI\n2ZUYJcT+HMsbW94g5XQKz3Z/lifCn8DL3avM+6z2ffIiItXdheILfLj7Q+bEz8GtlhvTe03n4XYP\nU9uldpXXopAXEakgJ/NPMn/rfOZvm094s3Deue8d7gy809JzdAp5EZFySjmdwhtb3mDlvpU8FPoQ\n34/7npDGIVaXBSjkRUTKLDE9kVk/zOKHoz8wudtkUqak4FP3+gkZraQ1Xkvh0KFDeHt7s3PnTgAy\nMjJo3LgxGzdutLgyEakqhmHw7c/fctfSu3ho5UPcGXgnh58/zCsRr9hdwING15TaBx98wJw5c9i2\nbRvDhw+nY8eO/PWvf73udfb+OUSkdIpLivl83+dEb47mku0SM/rMYFT7Ubi6uFZpHaXNlmoZ8hV1\nDqOsn3zYsGGkpqbi4uLC1q1bcXW9/i9ZIS/iGAovFbJ091Jm/zibZh7NmNFnBve1uQ9nJ2s6QmrE\nEEqrs3PixIkMGzaM999//4YBLyLVX86FHOZvnc9biW/RtVlXlg5fSt/mfa0uq9SqZUveSnl5eXTs\n2JEBAwbwzTffkJSUhJfX9Rc22PvnEJEby8jNYG78XBbtXMTgNoN5sc+LtPdpb3VZl9WI7horTZgw\ngYKCApYvX86TTz5JTk4On3766XWvs/fPISLXOnDmALM3z2ZV8irGdBjDC71eoEWDFlaXdZ0a0V1j\nlZiYGNauXUtSUhIAb775Jp06dWL58uVERUVZXJ2IlMXW9K1Eb45mw5ENPNPtGQ48e4BGdRpZXVaF\nUUu+kjjK5xBxRIZhsC51HdGbozlw5gDTe01nQpcJ1Ktdz+rSfpNa8iIiv8JWYmNV8iqiN0dzofgC\nL/Z+kaiwKEvmlKkqCnkRcXgXii+wdNdS/rblbzSu05iX73iZ+9veb9kwyKqkkBcRh3XuwjkWbFvA\n3IS5hDcNZ/HQxfRt3rdGLeqjkBcRh3Mi9wTzEubx/o73GdR6EN8++i0dfDtYXZYlFPIi4jB+Pvsz\nszfPZuW+lTwS9gjbn9hOYINAq8uylEJeRKq97Rnbid4czfdp3/N016dJmZJC47qNrS7LLijkRaRa\nMgyD7w5/x6zNs9h/ej8v9HyBxcMWV4thkFVJIS8i1YqtxMbq/auJ3hxNXlEeM/rMYHTYaIceBlke\nCnkRqRYuFl9k2e5lzP5xNt51vPlzvz8zJGhIjRgGWR4KeRGxa+cvnmfhtoXMTZhLR9+OfDD0A/o1\n71ejhkGWh0JeROxSVl4W8xLm8fftf+ee1vfwzehv6Niko9VlVTv6PacUZs+ezciRI6957rnnnmPq\n1KkWVSTieA6dPcTTXz9NyLshnL94nq2TtvLxgx8r4MtIE5SVQmZmJq1btyY9PR1PT0+Ki4vx8/Mj\nNjaWzp07X/Nae/4cIvZo54mdRG+OZl3qOp7u+jTP9njWLtdMtVqNmKDM6dWK6YszXi5dCDdp0oR+\n/fqxcuVKJk6cSGxsLI0bN74u4EXk1vwyDDJ6czT7Tu1jWs9pvD/kfTzcPKwuzWGoJV9KK1asYMGC\nBcTFxTFq1Cg6d+7MjBkzrnudvX8OESvZSmx8kfwFf/3xr+QV5fFi7xd5pMMjGgZ5C7QyVCUrLCzE\nz8+PjRs30qtXL5KTk/H397/udfb+OUSscKH4wuVhkI3qNGJGnxkMDRqqYZCloJCvApMmTSIhIQEf\nHx/WrVt3w9dUh88hUlXOXTjH/G3zmZcwjy5NuzCjzwwNgyyj0maLfnyWwbhx49izZw9jxoyxuhQR\nu5aRm8GL/36RVm+1Yu+pvax9dC3/Gv0v+rfor4CvItXyxKvVWrRogbu7OyNGjLC6FBG7dPWi2I92\neFSzQVpIIV9KJSUlvPHGG0RFRVGvniZCErlaYnoi0Zuj2XRkE5O7TXa4RbGrI4V8KeTn5+Pr60vL\nli2JjY21uhwRu2AYBt8e+pbozdGkZqcyvdd0lg1fRt3ada0uTdCJ10rjKJ9D5NcUlxSzcu9KojdH\nYzNsvNj7RUa1H4Wri6vVpTm0GnExlIhYp+BSAUt2LuGNLW/gX9+f1+56jfva3KcTqXZKIS8ityQr\nL4t3t77Lgm0L6BXQi48e/IjeAb2tLkt+g0JeRG4q5XQKb2x5g5X7VvK7dr9j0/hNBDUKsrosuUV2\nGfL6tU/EWoZh8MPRH/jblr+x5dgWJnebTMqUFE0YVg3Z3YlXEbFOcUkxq5NX87ctf+Ns4Vmm95rO\n2I5jqeNax+rS5D904lVESi2vKI8lO5cwJ34OTT2a8lKflxgaNBQXZxerS5NyUsiL1GCZeZm8nfA2\nf9/xd/q36K+TqQ6oUkM+MDCQ+vXr4+LigqurK4mJiZV5OBG5RftO7ePNLW+yKnkVo9uPZsuELbRu\n2NrqsqQSVGrIOzk5ERcXR8OGDSvzMCJyCwzD4Pu073lzy5tsy9jGM92e4eCzBzXtgIOr9O4anVwV\nsdaF4gt8kvQJc+PnYjNsPN/jeVY+tBJ3V3erS5MqUKmja1q1aoWnpycuLi48+eSTTJo06cqBNbpG\npFJl5mUyf+t8FmxfQHjTcKb2nEpkq0gNUa7m7Gp0zebNm2natCmnTp0iMjKS4OBg+vXrd/nPX3nl\nlcv3IyIiiIiIqMxyRGqEnSd2MjdhLl+lfMWodqOIGxdHSOMQq8uSMoqLiyMuLq7M76+ycfKvvvoq\n9erVY/r06eaB1ZIXqTC2Ehv/PPBP5sbP5eezPzOl+xQmdZmEdx1vq0uTCmY3LfmCggJsNhseHh7k\n5+ezdu1aXn755co6nEiNdP7ieZbsXMJbiW/RqE4jpvWcxoiQEZoJUi6rtJDPysrigQceAKC4uJhH\nHnmEgQMHVtbhRGqUw9mHeSvhLZbuXkrk7ZF89MBH9PTvqf52uY6mNRCpJgzDYP3h9byT+A4/HP2B\nxzs/zpTuU2ju2dzq0qQKlTY7FfIidu78xfMs272Md7e+Sy3nWkzpNoVHOjxCvdpafrImsps+eREp\nn/2n9/NO4jt8kvQJA1oNYMHgBfRv0V9dMlIqCnkRO2IrsfH1ga95Z+s7JGUlManLJH56+if86/tb\nXZpUUwp5ETtwuuA0i3YsYv62+TT1aMqUblMYGToSt1puVpcm1ZxCXsRC2zO2887Wd1idvJrhwcP5\n/OHP6dqsq9VliQNRyItUsfyifFbsWcHC7QvJys/i6a5Pc/DZgzSu29jq0sQBaXSNSBVJykpi4faF\nLN+znD4BfXgy/EnubX2vFuaQUtHoGhE7UnipkM/3fc7C7Qs5nHOYiZ0nsuvJXQR4BlhdmtQQasmL\nVIKU0yks3L6QD3/6kPCm4TzV9Snub3s/tZzVrpLyUUtexCJFtiJWJ69mwfYFJJ9KZnyn8SRMTKCV\nVyurS5MaTCEvUk7Jp5JZvHMxH/70Ie182vF016cZHjyc2i61rS5NRCEvUha5F3P5bO9nLNq5iMM5\nhxnXcRwbx2+krXdbq0sTuYb65EVukWEY/HjsRxbtXMTq/au5o8UdTOg8gUFtBqmvXaqMJigTqWCZ\neZl8uPtDFu9ajGEYTOg8gTEdx9CkXhOrS5MaSCEvUgGKS4pZc3ANi3YuYsORDTwY/CATukygl38v\nTRAmllLIi5SRYRjsztrNst3L+CTpE1p5tWJC5wk83O5hPNw8rC5PBNAQSpFSy8jN4OOfPubDnz4k\ntyiXMR3G6CSqOAy15KVGyi/K58v9X7Lsp2VsTd/KgyEPMrbjWPo274uzk7PV5Yn8KnXXiPyKEqOE\nuLQ4lu1eRkxKDL0DejO2w1iGBg3F3dXd6vJEbolCXuQqhmGwK3MXy/csZ8WeFTSq04gxHcYQFRal\n0TFSLSnkRTDnjvkl2ItsRYxqP4qo9lGE+YZZXZpIuSjkpcY6eu4on+75lOV7lpOZl8nD7R4mqn0U\n3f26a9ijOAyFvNQoJ/NPsnLvSpbvWc7+0/t5MORBotpH0b9Ff83TLg5JIS8OLysviy/3f8nnyZ+z\nNX0r97e9n6j2UUTeHqlJwcThKeTFIaWfT+eL5C9YlbyKXZm7uK/NfYwIGcGgNoOo41rH6vJEqoxC\nXhzGkZwjrEpexef7Pmf/6f0MCRrCyJCRRN4eyW21brO6PBFLKOSlWjt45iBfJH/B58mfk5aTxvCg\n4YwIHcFdLe9SV4wICnmpZmwlNhLSE/gq5StiUmI4d+Ecw4KGMTJ0JHcE3qEpfEX+i0Je7F5+UT7r\nUtcRkxLD1we+pkm9JgwLGsbQoKGENwvXtAIiN6GQF7uUmZfJ1we+JiYlhg1pG+ju152hQUMZ0nYI\nLb1aWl2eSLWhkBe7YCuxsS1jG2t+XsOan9dw4MwB7m19L0PbDuXe1vfi5e5ldYki1ZJCXixzMv8k\n3/78LWt+XsPaQ2tp6tGUQa0HcW/re+nbvK9OnIpUAIW8VBlbiY3E9MTLrfWDZw5yV8u7Lgd7gGeA\n1SWKOByFvFQawzBIzU5l/eH15pa6Hr/6fpdDvXdAb7XWRSqZQl4qVFZeFt8d/o71h9ezLnUdRbYi\nBrQawN0t7+buVnfjV9/P6hJFahSFvJRL7sVcNhzZwPrU9aw7vI5j544RERjBgJYDuLvV3QQ3CtaM\njiIWUshLqWQXZvPD0R/YeGQjG49uZO/JvXT363451MObheuCJBE7opCXmzqRe4JNRzex8chGNh3d\nRGp2Kj39e9K/eX/6t+hPd7/uWgpPxI4p5OWyEqOEA2cOEH88/nJr/XTBafo270v/Fv3p17wfXZp2\nwdXF1epSReQWKeRrsJwLOSQcTyD+eDzx6fEkHE/A8zZPevr3pLd/b+4IvIP2Pu01bYBINaaQryGK\nS4pJPpVM/PF4thzfQvzxeI6eO0p4s3B6+feip39Pevr31GLVIg5GIe+ALtkuse/UPraf2M6OEzvY\nfmI7P2X9hH99f3r49aCnf096+fcizDdMJ0lFHJxCvpq7WHzxcqBvz9jOjswd7Dm5h+aezQlvGk54\n03C6NO1C56adqe9W3+pyRaSKKeSriRKjhCM5R0g6mURSVpJ5ezKJ1OxUbve6nfBm4XRp0oXwZuF0\natKJerXrWV2yiNgBhbydMQyDzLxM9p/ez56Tey6H+d6Te6nvVp8w3zDCfP6z+YYR0igEt1puVpct\nInZKIW+RC8UXOHjmIClnUth/ej8pZ1JIOW3ed6vlRpB3EO0at6ODbwfCfMNo79Oehu4Ny3y8oiIo\nKID8/Gtvf7lfWAgXL5qv+7XbX3uupOTKZrNd+/hGzzk5gYvLlc3Z+drHVz/n6gq33XZrm4cH1K9v\nblffd3MzjylSEynkK9H5i+dJzU7lcPZhUrNTSc1O5VD2IQ6cOUBGbgYtvVoS3CiYIO8ggryDzPuN\ngvCs3ZDcXDh/3tyuvv9rz90swAsKwDCgbl1zq1Pn+lt3dzMM3dygdu0rt1ffv9Gtq+uVUP5l++/H\nVz/n5GTWYrNdCf9f7l+9/fL8pUvmD5QLF26+FRRAXt71301urrmfXwLf0xO8vaFRoyu3V9//5bZJ\nE/MHh0h1ZzchHxsby9SpU7HZbEycOJEZM2Zce2A7C3nDMMi5kMPx88c5fv44x84fIy0njUNnUzl4\nOpW0c6lcKL5AE7dWeDu3xNNoRZ2LrXArbIlbblvIbkl+rusNw7uwEOrVu3Gr9EbPeXjcPMDr1DFD\nuaYqKroS/Dk5cOaMuZ0+bW7/ff/UKcjKMv8OmjW78ebnB4GB0LixfksQ+2YXIW+z2QgKCmLdunX4\n+fnRrVs3li9fTkhISJkLLSvDgDPnLvDziSwOnzzJ0bNZpGUfJ/38cTILjnP60nFybMfJdToOhgu1\nC/1xzvfHyPHn0qmWXDrVCo9LrfCiFd7ujfFq4ESDBuDlBQ0amC1JT8+bh3fdumarV6xTUgJnz0J6\nOmRkXL8dOwZHjpi/RQQGQsuW5nb1/datzb9PESuVNjsrZVB1YmIirVu3JjAwEIBRo0YRExNzTciX\nxqVLcO4cZGfD6bOXOHY6h4yz2WTmZHMyN5vT+TlkF57lbNFJzttOkm9kUeiSRVHtLEpuOwm1LuB8\nwYfaRb64l/jgYfjTwNmfxm79aePuT0D9AJp7+dHMu/41Ae7lZbb+FNDVn7Pzla6cjh1//XXnzkFa\nmrkdPmzebtwIqalw6BA0bAhBQRAcbG6/3Pf3178TsU+VEvLp6ekEBFxZFcjf35+EhITrXvfUu5+Q\nk59PdkEeuRfyyb2QT/6lPAqK87lgy+eikUeRUx4ltXNwrpONcVs2hkshrjZPapd44e7kRV0XL+rf\n5oVnfS+C6/rQpF47/L3uIqChD618fLm9iS8+9T01Pa7cEk9P84fAjX4QlJTA0aOQkgL798O+fbB6\ntXn/3DkIDb3y3o4doUMHs7EgYqVKCflbDdQVH7+Oq3NtatdypUmb1gS1b4dnnSZ4e9TD26MujT3r\n4uNVj2ZeDWhYxwuv27zwcPPQ3CtiCWdns/smMBDuuefaPzt3Dvbsgd27ze3jj83Hv/zm0KkTdOsG\n3buDj48V1Ut1FRcXR1xcXJnfXyl98vHx8bzyyivExsYCMHPmTJydna85+WpvJ15FKlpJidnFs3s3\n7NoFiYmwdavZuu/eHXr0MG+7dDFPpovcCrs48VpcXExQUBDr16+nWbNmdO/e3bITryL2pKQEDh6E\nhAQz9BMSYO9eCAmB/v3hjjugXz9z6KfIjdhFyAOsWbPm8hDKCRMm8Ic//OHaAyvkRQBzRM+2beYJ\n3o0b4ccfoUULM/D79ze3JppMVP7DbkL+Nw+skBe5oeJi2LkTNmwwQ3/TJggIMM8DDBxotvR1YVfN\npZAXcTDFxWZL/9tvzS0pCfr2NQP/nnvMrh4NHqs5FPIiDi4nB9avh7VrITYWatWCYcPMrU8f87E4\nLoW8SA1iGObonZgYczt6FAYPNgP/nnvMq63FsSjkRWqwI0fgq6/MwE9MhMhIGDXKDH4N03QMCnkR\nAcy5elavhk8/NQN/8GAz8AcONGcclepJIS8i18nKglWrYMUKc1z+Aw/AY4+Zffg6aVu9KORF5KaO\nHYPly+Ef/zAn/3vsMRg71hymKfZPIS8it8QwzG6cJUvgs8/MuXUee8xs5Wscvv1SyItIqRUWwpdf\nmoG/c6cZ9k89BbffbnVl8t9Km52azlFEcHeHqChz7H18vNlP37MnDBoE//ynueSiVE9qyYvIDRUW\nwsqV8N57kJlptuyfeMJcOEWso5a8iFQId3fzhGx8vDkyZ/9+cwnEZ581p1CW6kEhLyK/KTzcHI2z\nZ4+5JGaPHjBypPkDQOybumtEpNTy8mDxYpg7F5o2hT/8wbzYSmPuK59G14hIlSkuhi++gNdeAxcX\n+POfYfhwLWpemRTyIlLlSkrMUTj/+79w8aIZ9iNHmsEvFUshLyKWMQxYs8YM+5wc+MtfzPly1LKv\nOAp5EbGcYZhz3v/lL1BQAK+/Dvfdpz77iqCQFxG7YRjm1Md//CN4ecHMmebyhVJ2CnkRsTs2G3z8\nMbz8MoSGmi37jh2trqp60sVQImJ3XFzMC6v274d77zW3xx+HEyesrszxKeRFpMq4uZlXzKakQOPG\nEBZmtuoLC62uzHEp5EWkytWvD9HRkJAA27dDSIi5gpV6cCue+uRFxHJxcTB1qjllwjvvQKdOVldk\nv9QnLyLVTkSE2aIfM8Zcg3baNMjNtboqx6CQFxG74OICTz5prkGbk2N24Xz2mbpwykvdNSJilzZt\ngsmToVkzswunTRurK7IP6q4REYfQrx/s2AGRkdCrl3kh1aVLVldV/aglLyJ2Ly3NXJXq9GlziuOa\nfGJWLXkRcTiBgfDtt+YY+4EDzTlxLl60uqrqQSEvItWCkxOMHw+7dkFSEnTpYo6zl5tTd42IVDuG\nYY68ef55c3qEV16B2rWtrqpqqLtGRByekxP87newe7fZqu/ZE/bts7oq+6SQF5Fqy9fXnMr46afh\njjvgrbfMVarkCnXXiIhD+PlnePRRc16cJUvAz8/qiiqHumtEpEZq3Rp++AH69jVPyn75pdUV2Qe1\n5EXE4cTHm2vLDh9uznbp5mZ1RRVHLXkRqfF69oSdO+HIEejTBw4dsroi6yjkRcQheXnBF1/AuHFm\n6H/2mdUVWUPdNSLi8LZvN4dcRkbCnDlw221WV1R26q4REfkv4eFm0J8+bU58dvSo1RVVHYW8iNQI\nnp5ml83DD0OPHvD991ZXVDXUXSMiNc769eaY+t//3lyFysnJ6opuXWmzUyEvIjXSkSMwYoQ5vn7R\nIqhb1+qKbo365EVEbkGLFubqU+7u5uib1FSrK6ocCnkRqbHc3c1FSJ56Cnr3NkPf0SjkRaRGc3KC\nZ56BZcvM7pt//MPqiipWpYT8K6+8gr+/P507d6Zz587ExsZWxmFERCrMwIGwYQP83//BjBmOM5tl\npZx4ffXVV/Hw8OCFF1749QPrxKuI2KEzZ+DBB80rZj/6COrVs7qia9nNiVcFuIhUR97e8O9/m7d9\n+8Lx41ZXVD6VFvJvv/02HTt2ZMKECeTk5FTWYUREKlzt2vDBBzB6tHlCds8eqysquzJ310RGRpKZ\nmXnd86+99ho9e/akcePGAPzlL3/hxIkTLFq06NoDOznx8ssvX34cERFBREREWUoREak0n3xiXjD1\n2Wfm6lNVLS4ujri4uMuPX331Vfu6GCotLY0hQ4aQlJR07YHVJy8i1cT69RAVBe++Cw89ZG0tdtEn\nf+LEicv3V69eTVhYWGUcRkSkSgwYAGvXmi36efOsrqZ0KqUlP3bsWHbt2oWTkxMtW7Zk4cKF+Pr6\nXntgteRFpJo5cgTuvRfuv99cccrZgiuNNHeNiEglOnsWhgy5MudNrVpVe3y76K4REXFUDRuaQyyz\nssxpiy9etLqim1PIi4iUUp06EBNjdtcMGQL5+VZX9OsU8iIiZeDmBitWgJ+fOSWCvV4OpJAXESmj\nWrXMfvmuXeHOO+HkSasrup5CXkSkHJydYe5cs9umf3/7mwahis8Li4g4Hicn+J//AQ8PiIgw148N\nCLC6KpNCXkSkgvz+92bL/s474bvvoHlzqytSyIuIVKjp068E/fffWx/0CnkRkQo2bZp5+0vXTYsW\n1tWikBcRqQTTpl3bdRMYaE0dCnkRkUry/PPmSdk774S4OGta9Ap5EZFK9Nxz5u2AAbBxIzRrVrXH\nV8iLiFSy554zpz6IjDRb9P9ZU6lKaBZKEZEq8qc/wZo1Zh99gwZl24emGhYRsVOGAVOnwtat5iIk\n9eqVfh8KeRERO1ZSAk88Aamp8K9/gbt76d6vkBcRsXM2Gzz6KOTmwhdfQO3at/5eLRoiImLnXFxg\n2TLz9rHHzNZ9ZVHIi4hYwNXVnI/++HFzKoTK6thQyIuIWMTd3Vxhat06mD27co6hcfIiIhby8oLY\nWOjTB3x8zO6biqSQFxGxmJ+fGfQREeaFUoMHV9y+1V0jImIHgoPNrpvx4yE+vuL2q5AXEbETPXrA\n0qUwfDgkJ1fMPhXyIiJ2ZNAgmDXL7LKpiIXBFfIiInbmscfMi6WGDYPCwvLtS1e8iojYIcMwg/7S\nJXM8vfN/muS64lVExAE4OcHixXDiBPzxj2Xfj0JeRMROubnBl1+a89u8/37Z9qFx8iIidszb25yt\nsl+/sq0Tqz55EZFqYNMmGDECTp1Sn7yIiMPp1w/mzCn9+9SSFxGpRjS6RkRELlPIi4g4MIW8iIgD\nU8iLiDgwhbyIiANTyIuIODCFvIiIA1PIi4g4MIW8iIgDU8iLiDgwhbyIiANTyIuIODCFvIiIAytz\nyK9cuZJ27drh4uLCjh07rvmzmTNn0qZNG4KDg1m7dm25ixQRkbIpc8iHhYWxevVq+vfvf83z+/bt\n49NPP2Xfvn3ExsYyefJkSkpKyl2oI4uLi7O6BLuh7+IKfRdX6LsouzKHfHBwMG3btr3u+ZiYGKKi\nonB1dSUwMJDWrVuTmJhYriIdnf4BX6Hv4gp9F1fouyi7Cu+Tz8jIwN/f//Jjf39/0tPTK/owIiJy\nC266kHdFJ1X+AAAFe0lEQVRkZCSZmZnXPf/6668zZMiQWz6Ik5NT6SsTEZHyM8opIiLC2L59++XH\nM2fONGbOnHn58T333GPEx8df9z5AmzZt2rSVYSuNm7bkb5Vx1XqDQ4cOZfTo0bzwwgukp6dz8OBB\nunfvftP3iIhI5Shzn/zq1asJCAggPj6ewYMHM2jQIABCQ0N5+OGHCQ0NZdCgQbz33nvqrhERsYiT\noSa1iIjDsuSK19jYWIKDg2nTpg3R0dFWlGAXjh07xp133km7du1o3749b731ltUlWc5ms9G5c+dS\nndh3RDk5OYwcOZKQkBBCQ0OJj4+3uiTLzJw5k3bt2hEWFsbo0aO5ePGi1SVVmccffxxfX1/CwsIu\nP3f27FkiIyNp27YtAwcOJCcn56b7qPKQt9lsTJkyhdjYWPbt28fy5ctJTk6u6jLsgqurK3PmzGHv\n3r3Ex8fz7rvv1tjv4hfz5s0jNDS0xnfxPf/889x3330kJyfz008/ERISYnVJlkhLS+P9999nx44d\nJCUlYbPZWLFihdVlVZnx48cTGxt7zXOzZs0iMjKSAwcOMGDAAGbNmnXTfVR5yCcmJtK6dWsCAwNx\ndXVl1KhRxMTEVHUZdqFJkyZ06tQJgHr16hESEkJGRobFVVnn+PHjfPPNN0ycOLFGn5g/d+4cmzZt\n4vHHHwegVq1aeHp6WlyVNerXr4+rqysFBQUUFxdTUFCAn5+f1WVVmX79+uHl5XXNc1999RXjxo0D\nYNy4cXz55Zc33UeVh3x6ejoBAQGXH+tiKVNaWho7d+6kR48eVpdimWnTpjF79mycnWv2vHmHDx+m\ncePGjB8/ni5dujBp0iQKCgqsLssSDRs2ZPr06TRv3pxmzZrRoEED7r77bqvLslRWVha+vr4A+Pr6\nkpWVddPXV/n/ppr+a/iN5OXlMXLkSObNm0e9evWsLscSX3/9NT4+PnTu3LlGt+IBiouL2bFjB5Mn\nT2bHjh3UrVv3N38ld1SHDh1i7ty5pKWlkZGRQV5eHh9//LHVZdkNJyen38zUKg95Pz8/jh07dvnx\nsWPHrpkGoaa5dOkSI0aM4NFHH2X48OFWl2OZH3/8ka+++oqWLVsSFRXFd999x9ixY60uyxL+/v74\n+/vTrVs3AEaOHHndTK81xbZt2+jduzfe3t7UqlWLBx98kB9//NHqsizl6+t7eSaCEydO4OPjc9PX\nV3nId+3alYMHD5KWlkZRURGffvopQ4cOreoy7IJhGEyYMIHQ0FCmTp1qdTmWev311zl27BiHDx9m\nxYoV3HXXXSxbtszqsizRpEkTAgICOHDgAADr1q2jXbt2FldljeDgYOLj4yksLMQwDNatW0doaKjV\nZVlq6NChLF26FIClS5f+duOw9BMZlN8333xjtG3b1rj99tuN119/3YoS7MKmTZsMJycno2PHjkan\nTp2MTp06GWvWrLG6LMvFxcUZQ4YMsboMS+3atcvo2rWr0aFDB+OBBx4wcnJyrC7JMtHR0UZoaKjR\nvn17Y+zYsUZRUZHVJVWZUaNGGU2bNjVcXV0Nf39/Y/HixcaZM2eMAQMGGG3atDEiIyON7Ozsm+5D\nF0OJiDiwmj2MQUTEwSnkRUQcmEJeRMSBKeRFRByYQl5ExIEp5EVEHJhCXuQq586dY/78+VaXIVJh\nFPIiV8nOzua9996zugyRCqOQF7nKSy+9xKFDh+jcuTMzZsywuhyRctMVryJXOXLkCPfffz9JSUlW\nlyJSIdSSF7mK2jziaBTyIiIOTCEvchUPDw9yc3OtLkOkwijkRa7i7e1Nnz59CAsL04lXcQg68Soi\n4sDUkhcRcWAKeRERB6aQFxFxYAp5EREHppAXEXFgCnkREQemkBcRcWAKeRERB/b/OrZJOw6jpQgA\nAAAASUVORK5CYII=\n", "text": [ "" ] } ], "prompt_number": 2 }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Exercise" ] }, { "cell_type": "markdown", "metadata": {}, "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", "
" ] }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Additional Documentation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Further information is available at: \n", "http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html \n", "http://docs.scipy.org/doc/scipy/reference/tutorial/integrate.html" ] } ], "metadata": {} } ] }