diff --git a/Week 5/9 Ordinary Differential Equations.ipynb b/Week 5/9 Ordinary Differential Equations.ipynb index 2365c9e..c9eee81 100644 --- a/Week 5/9 Ordinary Differential Equations.ipynb +++ b/Week 5/9 Ordinary Differential Equations.ipynb @@ -90,7 +90,8 @@ }, "outputs": [], "source": [ - "import numpy as np" + "import numpy as np\n", + "from matplotlib import pyplot as plt" ] }, { @@ -139,7 +140,7 @@ "def integrator(x, y0, f, phi):\n", " \"\"\"\n", " Numerically solves the initial value problem given by ordinary differential equation\n", - " f(x, y) = y' with initial value y0 using the Euler method.\n", + " f(x, y) = y' with initial value y0 using the integration scheme provided by phi.\n", "\n", " Args:\n", " x: size n + 1 numerical array\n", @@ -147,25 +148,27 @@ " f: a callable function with signature (x, y), with x and y the current state\n", " of the system\n", " phi: a callable function with signature (x, y, f, i), with x and y the current state\n", - " of the system, i the step number, and f as above\n", + " of the system, i the step number, and f as above, representing the integration\n", + " scheme to use\n", "\n", " Returns:\n", " An n + 1 numerical array representing an approximate solution to y' = f(x, y)\n", " given initial value y0 and steps from x\n", " \"\"\"\n", " \n", - " y = np.zeros(len(x))\n", - " y[0] = y0\n", + " eta = np.zeros(len(x))\n", + " eta[0] = y0\n", " \n", - " for i in range(len(y)):\n", - " y[i] = y[i - 1] + (x[i] - x[i - 1])*f(x[i - 1], y[i - 1])\n", + " for i in range(1, len(eta)):\n", + " h = x[i] - x[i - 1]\n", + " eta[i] = eta[i - 1] + h*phi(x[i - 1], eta[i - 1], f, i)\n", " \n", - " return y" + " return eta" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 3, "metadata": { "deletable": false, "nbgrader": { @@ -200,10 +203,7 @@ " given initial value y0 and steps from x\n", " \"\"\"\n", " \n", - " for i in range(len(y)):\n", - " y[i] = y[i - 1] + (x[i] - x[i - 1])*f(x[i - 1], y[i - 1])\n", - " \n", - " return y\n", + " return f(x, y)\n", " " ] }, @@ -241,7 +241,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 4, "metadata": { "deletable": false, "nbgrader": { @@ -256,12 +256,42 @@ "task": false } }, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAD4CAYAAAAdIcpQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAo0UlEQVR4nO3de3xV9Znv8c+TK5BAyI2ACUiQKAIqSkTUmbYqVWydoq1ttZ2WtrbMsTp12jkz1fGcdqadnrHTi1NbtcNUp9g6RYfakWlRCmpt6xQUvHATJNzvBEIC4ZLbfs4f+4fdhoSt7CRr7/B9v177lbWe9VtrPwsxD2v9fmv9zN0RERE5mayoExARkfSnYiEiIkmpWIiISFIqFiIikpSKhYiIJJUTdQK9payszEePHh11GiIiGWX58uX73L28c7zfFovRo0ezbNmyqNMQEckoZralq7huQ4mISFIqFiIikpSKhYiIJKViISIiSalYiIhIUhlTLMxsupmtM7M6M7sz6nxERE4nGVEszCwbuB+4FhgP3Gxm46PNSkTk9JEpz1lMAercfSOAmc0FZgBrevqLXnv2cY5s7uL5DLO3f5Au2hoGdH0Mt85twTu1tbDuCcdOmlFXOb8l1vURjke923NOdtwT2yZuPn5unffwLtp2e5Zdfp91s0vXbTsfwjGwLCwrO/7TsrHs7Ddj8Xg8lpVlYDkhHt+eFdpYVjaWnUNW1h/3y8rKJis7B8vJIyc3n5y8fHLzBpCTm0du3gDycnPjxxRJU5lSLCqBbQnr24FLOjcys1nALIBRo0ad0hcdW/MUl+574pT2FTlVHW4cI4c2cmizHNrJpZ0c2i2Hdsulw3JotzxaswbQnj2A9qyBdOQMxHMGEssdBLkDIXcQlldAVt4gsvMLyRpQSE7BUPILSygoKqVwSAlDCgeRm50RNxQkzWRKsXhb3H02MBugtrb2lGZ1uvjWHxHjR52P2833xbrKoYsYQBdx/PjGpG0TNnZa7LT/8fUkbd+y11vadpXniW2727/ThhODx/c/4Wu6P+6JOfmJf0Jd5RXinVtblw0BjxGLxYjFOoh1dOCx9rAcw2PtuMeIdbTjHR3EPMQ6Ynisg1isA2IdxLwD7+iIH6ujA/d43GMxiLXjHW14eyve0QodrQnLbVj4SUcrFmvDYq1kxdqwWBtZsTZyYi0MjB0jt62JvNhR8v0Y+d7CAG8h1zq6OKETHfZ89lHI4awCjmYV0pIzmGO5Q2nJL8ULyskqrCB3aAWDikcwuLSS4rIKSgoHkK2rntNephSLHcDIhPWqEOtxWdnZvXFYkV7l7a20HTvMscOHOHb0EK1Hm2k9cpD2w420HT5A+5EDxI404seayDrWSHbbQfJaDzK4vZ6Ctg0MbW4kd3/7Ccdt82x2U8zerHKa8oZztKAKH1JFbumZFFSMoaRyLKPKixiUlym/SuRUZcp/4ZeAGjOrJl4kbgI+Fm1KIunDcvLIK8wjr7CYIadyAHf86AEON+zi4L4dHN6/i9am3XQc2kN2806KDu9kZMtqihueJ7shBpvju7V5Nlu8gh05VTQVjKGjpIb84edQUn0BZ1VWUD44vwfPUqKUEcXC3dvN7HZgIZANPOzuqyNOS6T/MMMGlVA4qITCqgndt+toJ9a0g6ZdG2javYG2PevJ2v8G5xzaSNmh5eQc6oAt0LHE2OBn8FLOWJqKxpNVeSFlYydz7uhKzhg6sO/OS3qMdXc/PtPV1ta63jor0oc62vCGTTRuXcWhzS/DrtcoalxNUfv++GY3Vnk1a3In0Dz8EorOeTeTzqlmbHmhRoKlETNb7u61J8RVLESkVx3aTcvWl2l44wVsyx8obVpBrrcBsDp2JkuyJ3No5BWcOendvOvs4ZQW6tZVlFQsRCQ9tB3Ddyyn8fXf0Lb+WUobXiGbDg54Ic/HLmB18VWUTrqW910wmlGlg6LO9rSjYiEi6eloI7G6Z2lc8UsGblrMwPYmDvognu64mFUl0xg5eTrXTRrJiCL1dfQFFQsRSX8dbbDxeQ6/PJfcNxaQ13GYHV7K4x1XsPXMDzH9sou4ctwwPVjYi1QsRCSztB2DN57m6JKHGLjtt7STxeKOyTyedwMXXnY1fz71TIoL8qLOst9RsRCRzLV/A7FlP6Z9+SPktTayNDaOh/16hk++js++6yxGlqhvo6eoWIhI5mtphpcfoe3395F7eBerfTT/3H4zlbXv5wtX1jC8aEDUGWY8FQsR6T/aW2Hlf9L+3D+Rc3AbL8Qm8p3Yx5h86RV8/j1jdXsqBd0VC/USiUjmycmDCz9OzheWwzX/xNRBO3gi5++oWXIXN3x7Pj9dsoWOWP/8h3BUVCxEJHPl5MOlnyf7r16Dy77Ah3N/z3z7Iq/99w+44Qe/47VtjVFn2G+oWIhI5htQBFd/HfuL3zG4agLfyp3N3Qfu5tYH5vOthWtpaX97r3CX7qlYiEj/UTEe+9QCuO5fmJKzgcUD72Lr8z/hA99/gVU7mqLOLqOpWIhI/5KVBbWfxv7X7xh0xrl8P+8H3H7wu9z8wG/46ZIt3U5mJienYiEi/VPpWfDpp+HdX+bP/Dn+e9DX+OGTz3LH3Fdpbjlxoic5ORULEem/snPgir+Djz3Omdn1LC74CgdWLuSG+19gW8ORqLPLKCkVCzP7sJmtNrOYmdV22naXmdWZ2TozuyYhPj3E6szszoR4tZktDfHHzCwvxPPDel3YPjqVnEXkNHT2Ndis3zCgZCRz8r/FlIMLueGBF3hl64GoM8sYqV5ZrAI+CPw2MWhm44lPfToBmA48YGbZZpYN3A9cC4wHbg5tAb4J3OvuY4EDwC0hfgtwIMTvDe1ERN6ZkjHwmafJGn053+B+/sJ+wU2z/8CClbuiziwjpFQs3P11d1/XxaYZwFx3b3H3TUAdMCV86tx9o7u3AnOBGWZmwJXAvLD/HOD6hGPNCcvzgKtCexGRd2bAEPj4PDj/o3yu7VHuHfwf3PYfy3l82baoM0t7vTUHdyWwJGF9e4gBbOsUvwQoBRrdvb2L9pXH9wlzcTeF9vt6J3UR6ddy8uCGf4WCct73hx/wb6Xw2XnQ0tbBJy4dHXV2aStpsTCzxcDwLjbd7e5P9nxKp87MZgGzAEaNGhVxNiKStszg6n+ErGymvfA9/r0cPv3kTbS0x/jsn46JOru0lLRYuPu0UzjuDmBkwnpViNFNfD8w1MxywtVFYvvjx9puZjlAUWjfVa6zgdkQf5HgKeQtIqcLM5j2DwBc8cL3+PeKHD79qxspyM/h5in6x2ZnvTV0dj5wUxjJVA3UAC8CLwE1YeRTHvFO8Pkef0rmOeDGsP9M4MmEY80MyzcCz7qeqhGRnnC8YFx6O1c0/YJvVTzD3b9YyVPq9D5BSn0WZnYD8H2gHPiVmb3q7te4+2ozexxYA7QDt7l7R9jndmAhkA087O6rw+G+DMw1s38EXgEeCvGHgJ+YWR3QQLzAiIj0DDN479eheS8fXvkwu8uKuGNuFkMG5nL52LKos0sbms9CRATic2T87KP4xuf58qC/Z0HzOTzx+cs4u2Jw1Jn1Kc1nISJyMjl58JFHsLIa7un4LmNz6/nsnGU0HG6NOrO0oGIhInJc/mC4+Wdk4fys8HscOniAW3+6nNb2WNSZRU7FQkQkUckY+PCPGdi0gf8eNZelm/Zzz1Nro84qcioWIiKdnXUFXPV/qdq5kO/XvMbDL2zi6VW7o84qUioWIiJduewOOOsqrtv5PT4wvIG/mffaaf2mWhULEZGuZGXBDf+KDSzmO3YvA2jh9p+9QnvH6dl/oWIhItKdwnL44L+Se2ADj41dzGvbGnngNxuizioSKhYiIicz5j1w8ecYU/cIX6yp575n1p+W83mrWIiIJDPt76H4TG4/dC+VBTG++NirHGvriDqrPqViISKSTH4hXP8g2Y1b+OmYxazf28wDz9VFnVWfUrEQEXk7zrwMJn+KkW88wq3jjvLg8xuo29scdVZ9RsVCROTtuuorMLCYL7U8yKBc4+5frKS/vl+vMxULEZG3a1AJXP11cnctY/bEtSzd1MDPX96RfL9+QMVCROSduOBmGHUZUzbcx59U5XLPU2tpbmlPvl+GU7EQEXknzODae7CjB/jOGc+yr7nltOjsVrEQEXmnRlwAF9xExeqH+cx440e/39TvXwWSUrEws2+Z2VozW2FmvzCzoQnb7jKzOjNbZ2bXJMSnh1idmd2ZEK82s6Uh/liYdpUwNetjIb7UzEankrOISI+48v+AGf879z/JMvjm0/37zbSpXlksAia6+/nAG8BdAGY2nvj0pxOA6cADZpZtZtnA/cC1wHjg5tAW4JvAve4+FjgA3BLitwAHQvze0E5EJFpFVXDp7Qxa9wT/98Jj/HLFLl7b1hh1Vr0mpWLh7r929+M9O0uAqrA8A5jr7i3uvgmoA6aET527b3T3VmAuMMPMDLgSmBf2nwNcn3CsOWF5HnBVaC8iEq3L74CBxXz08KMUD8rlO4veiDqjXtOTfRafAZ4Ky5XAtoRt20Osu3gp0JhQeI7H33KssL0ptD+Bmc0ys2Vmtqy+vj7lExIROakBQ+CyvyRnwyK+MukIv32jnpc2N0SdVa9IWizMbLGZreriMyOhzd1AO/BobyabjLvPdvdad68tLy+PMhUROV1MmQUDS/hA4yOUFebz7YXr+uWDekmLhbtPc/eJXXyeBDCzTwHXAR/3P/4J7QBGJhymKsS6i+8HhppZTqf4W44VtheF9iIi0csfDJd/geyNz/APFzazdFMDL9T1v19RqY6Gmg78LfABd08cNzYfuCmMZKoGaoAXgZeAmjDyKY94J/j8UGSeA24M+88Enkw41sywfCPwrPfHsi0imeviz8GgUq7dP4cRRQP4/rPro86ox6XaZ/EDYDCwyMxeNbMfArj7auBxYA3wNHCbu3eEPofbgYXA68DjoS3Al4EvmVkd8T6Jh0L8IaA0xL8EvDncVkQkLeQXwtTPk7XxWf5mUjtLNzXw8tYDUWfVo6y//iO9trbWly1bFnUaInK6OHoAvjuB9rPfx+Q1H+GS6hJmf7I26qzeMTNb7u4nJK4nuEVEesLAYpg8k5zXf8Ftkwfw6zV7qNt7KOqseoyKhYhIT5l6K7jzCRYwIDeLHz6/MeqMeoyKhYhITxk6CiZ+iIErfsKnLizmv17Zwd6Dx6LOqkeoWIiI9KTLbofWZmYNWUKHO48u3Rp1Rj1CxUJEpCeNuACqplCy5hGuqCnl0aVbaW2PRZ1VylQsRER62sWfhYYN3DFmJ/uaW1iwclfUGaVMxUJEpKdNuB4GlXL+rnlUlxXw4//ZHHVGKVOxEBHpaTn5cNEnsTee4tZJeby6rTHjX1+uYiEi0htqPwPuzIgtYlBeNv+R4R3dKhYiIr1h6CioeS/5K+dy3cRh/HLFTg63tCffL02pWIiI9JZJH4dDO/ls5RYOt3bwqxWZ29GtYiEi0lvOuRYGFlOzcz5jygt4bNm25PukKRULEZHekpMP530YW/srPnFBEcu3HKBub3PUWZ0SFQsRkd406ePQ0cKN+UvJyTL+M0OvLlQsRER604gLoGIig19/jCvHDePnL++gvSPznuhOdaa8r5vZijDx0a/N7IwQNzO7z8zqwvaLEvaZaWbrw2dmQnyyma0M+9xnZhbiJWa2KLRfZGbFqeQsItKnzGDSx2Dny3yipoV9zS38YWPmTbua6pXFt9z9fHefBPwS+EqIX0t8KtUaYBbwIMR/8QNfBS4BpgBfTfjl/yDwuYT9pof4ncAz7l4DPINmyhORTDPhg4Ax9ehvGZyfw5Ov7ow6o3cspWLh7gcTVguA49PuzQAe8bglwFAzGwFcAyxy9wZ3PwAsAqaHbUPcfUmYX/sR4PqEY80Jy3MS4iIimWHICDjzcnLXPMH0CRU8vWo3x9o6os7qHUm5z8LMvmFm24CP88cri0ogsRdne4idLL69izhAhbsfH5y8G6g4SS6zzGyZmS2rr68/xTMSEekFEz8I+97g5tHNNLe089zavVFn9I4kLRZmttjMVnXxmQHg7ne7+0jgUeD23kw2XHV0O2m4u89291p3ry0vL+/NVERE3pnxM8CymdT0DGWF+Rl3KyppsXD3ae4+sYvPk52aPgp8KCzvAEYmbKsKsZPFq7qIA+wJt6kIPzOrHIuIABSUwZj3kLX65/zZ+cN5dt1emo62RZ3V25bqaKiahNUZwNqwPB/4ZBgVNRVoCreSFgJXm1lx6Ni+GlgYth00s6lhFNQngScTjnV81NTMhLiISGaZ+CFo3MLNVftobY/x69W7o87obUu1z+KecEtqBfFf/HeE+AJgI1AH/BvweQB3bwC+DrwUPl8LMUKbH4V9NgBPHf8O4L1mth6YFtZFRDLPuPdDdh41e3/NGUUDWLh6T9QZvW05qezs7h/qJu7Abd1sexh4uIv4MmBiF/H9wFWp5CkikhYGDoXqd2HrFnDNhI/x6IvbONzSTkF+Sr+K+4Se4BYR6Uvj3g8HNnFD5UFa22P8Zl1mjNxUsRAR6UvnvA+Aic0vUFqQx9MZ0m+hYiEi0pcGD4eqi8la9yveO76C59bupaU9/R/QU7EQEelr494PO1/hA9VOc0s7/1OX/u+KUrEQEelr464D4OLWJQzOz+GpVek/g56KhYhIXyurgdIact9YwHvGDePZtXuJxbp9OUVaULEQEYnCuPfB5t9z9VmD2NfcyqqdTVFndFIqFiIiURj7Xoi18+7cNZjBc2vTewitioWISBRGXgJ5gxmy/XkmjRzKs+vS+7V3KhYiIlHIyYMx74a6xVxxdjkrtjeyr7kl6qy6pWIhIhKVsdOgaRvXDm/CHZ5P46e5VSxERKIydhoAZzUtoXxwPs+l8a0oFQsRkagMHQnl48jasJj3nF3Ob9+op70jFnVWXVKxEBGJ0thpsOV/eO/YQg4ea+eVbY1RZ9QlFQsRkSiNnQYdrVye/Tpm8ELdvqgz6lKPFAsz+2szczMrC+tmZveZWZ2ZrTCzixLazjSz9eEzMyE+2cxWhn3uCzPmYWYlZrYotF8UZtgTEekfRl0KOQMo2PkC51UW9d9iYWYjic+StzUhfC1QEz6zgAdD2xLgq8AlwBTgqwm//B8EPpew3/QQvxN4xt1rgGfCuohI/5A7AEZOgU2/4/KxZbyytZHDLe1RZ3WCnriyuBf4WyDxxSYzgEc8bgkw1MxGANcAi9y9wd0PAIuA6WHbEHdfEmbZewS4PuFYc8LynIS4iEj/UP0u2LOS91Rl0R5zXtzUkHyfPpZSsTCzGcAOd3+t06ZKYFvC+vYQO1l8exdxgAp3P/5Kxt1ARSo5i4iknep3A3BhbBV5OVn8Pg1vRSWd+NXMFgPDu9h0N/B3xG9B9Ql3dzPr9tWMZjaL+G0vRo0a1VdpiYik5owLIa+QvK2/5+LRH0nLfoukVxbuPs3dJ3b+ABuBauA1M9sMVAEvm9lwYAcwMuEwVSF2snhVF3GAPeE2FeFnt0+tuPtsd69199ry8vJkpyYikh6yc+HMy2DTb7l8bBlrdx+i/lB6vfrjlG9DuftKdx/m7qPdfTTxW0cXuftuYD7wyTAqairQFG4lLQSuNrPi0LF9NbAwbDtoZlPDKKhPAk+Gr5oPHB81NTMhLiLSf1S/C/av593D453b/7Mhva4ueus5iwXErzzqgH8DPg/g7g3A14GXwudrIUZo86OwzwbgqRC/B3ivma0HpoV1EZH+pfpdAJx77FWGDMhJu1tRSfss3q5wdXF82YHbumn3MPBwF/FlwMQu4vuBq3oqTxGRtFRxHgwYStbm3zGl+tNpNyJKT3CLiKSDrCwY/Sew+XdMqS5m8/4j7D14LOqs3qRiISKSLs68DBq3cNmweL/Fi5vT5+pCxUJEJF2MmgrAuLbVDMrL5qU0uhWlYiEiki6Gnw+5g8jZvpSLRhWzVMVCREROkJ0LlZNh6x+YUl3Cuj2HaDrSFnVWgIqFiEh6GXUp7F7JJZV5uMOyLelxdaFiISKSTkZNBY9xYVYdudmWNp3cKhYiIumk6mKwLPJ2vMj5VUPT5nkLFQsRkXQyYAhUTHiz32Ll9iaOtnZEnZWKhYhI2hl1KWxfxsUjB9Mec1buaIo6IxULEZG0M2oqtB3movz4y7df2Xog4oRULERE0k/VxQAMbXiNM0sH8crWxmjzQcVCRCT9FI2EgmGwYzkXjhzKy1sPEH8/a3RULERE0o1Z/OG8Hcu5cFQxew+1sKsp2pcKqliIiKSjysmw7w0mVxhA5LeiVCxERNJR1WQAzunYQH5OVuSd3CkVCzP7ezPbYWavhs/7ErbdZWZ1ZrbOzK5JiE8PsTozuzMhXm1mS0P8MTPLC/H8sF4Xto9OJWcRkYxwxoUA5O5+mfMqi3hlW2Ok6fTElcW97j4pfBYAmNl44CZgAjAdeMDMss0sG7gfuBYYD9wc2gJ8MxxrLHAAuCXEbwEOhPi9oZ2ISP82sBhKx4Z+i6Gs3NFEa3sssnR66zbUDGCuu7e4+ybi82pPCZ86d9/o7q3AXGCGmRlwJTAv7D8HuD7hWHPC8jzgqtBeRKR/q6yF7cu4cORQWttjvL7rYGSp9ESxuN3MVpjZw2ZWHGKVwLaENttDrLt4KdDo7u2d4m85VtjeFNqfwMxmmdkyM1tWX1+f+pmJiESpcjIc3svkkiMAvBxhv0XSYmFmi81sVRefGcCDwFnAJGAX8J3eTffk3H22u9e6e215eXmUqYiIpK4y3sld0bSKYYPzWbk9utd+5CRr4O7T3s6BzOzfgF+G1R3AyITNVSFGN/H9wFAzywlXD4ntjx9ru5nlAEWhvYhI/zZ8ImTnwY7lnFf5gUjfEZXqaKgRCas3AKvC8nzgpjCSqRqoAV4EXgJqwsinPOKd4PM9/mjic8CNYf+ZwJMJx5oZlm8EnvWoH2UUEekLOflQMRF2vsLEyiI21DdzpLU9+X69INU+i382s5VmtgK4AvgigLuvBh4H1gBPA7e5e0e4argdWAi8Djwe2gJ8GfiSmdUR75N4KMQfAkpD/EvAm8NtRUT6vREXwK4VnHfGEGIOa3ZG08md9DbUybj7J06y7RvAN7qILwAWdBHfSHy0VOf4MeDDqeQpIpKxRpwPy/+dSYMbAVi5o4na0SV9noae4BYRSWcjLgCg9NBaygrzI+u3ULEQEUlnwyaAZWO7V3Je5RBWqViIiMgJcgdA+Tmw6zXOqyyibm9zJNOsqliIiKS7ERfA7hVMrCyKd3JH8CS3ioWISLobfj407+GCofE5LaK4FaViISKS7kacD8Cww+soK8yLpJNbxUJEJN0NPw8AC7eidGUhIiInGlAExdWwawUTzyhi/d5mjrX1bSe3ioWISCYYcQHseo1zRwyhI+bU7W3u069XsRARyQQjzofGLYwviU+AtHb3oT79ehULEZFMUBHvtxjVtpkBuVl9PhGSioWISCaoiM9AnV2/hnMqBrN2t4qFiIh0NqQS8otg7xrGDR/C67sO0ZezNahYiIhkArP41cWeNYwbMZiGw63UN7f02derWIiIZIph42Hv65w7fDAAr+/qu07ulIuFmf2lma01s9Vm9s8J8bvMrM7M1pnZNQnx6SFWZ2Z3JsSrzWxpiD8WZtIjzLb3WIgvNbPRqeYsIpKRhp0LLU2cWxAvEmv7sJM71WlVrwBmABe4+wTg2yE+nviUqROA6cADZpZtZtnA/cC1wHjg5tAW4JvAve4+FjgA3BLitwAHQvze0E5E5PRTMQGAooPrGVE0oE+Hz6Z6ZXErcI+7twC4+94QnwHMdfcWd98E1BGfBW8KUOfuG929FZgLzDAzA64E5oX95wDXJxxrTlieB1wV2ouInF6GnRv/uXc144YP7tPhs6kWi7OBPw23h543s4tDvBLYltBue4h1Fy8FGsMc3YnxtxwrbG8K7U9gZrPMbJmZLauvr0/x1ERE0szA4vioqD1rOHfEEDbUN9PaHuuTr046B7eZLQaGd7Hp7rB/CTAVuBh43MzG9GiG74C7zwZmA9TW1vbdmDIRkb4SOrnHXTaEtg5nQ30z544Y0utfm7RYuPu07raZ2a3AEx4f7PuimcWAMmAHMDKhaVWI0U18PzDUzHLC1UNi++PH2m5mOUBRaC8icvoZdi5sep5zywcA8Pqug31SLFK9DfVfwBUAZnY2kAfsA+YDN4WRTNVADfAi8BJQE0Y+5RHvBJ8fis1zwI3huDOBJ8Py/LBO2P6s9+WTKCIi6aRiAnS0MjprN7nZxvo+eqFg0iuLJB4GHjazVUArMDP8Il9tZo8Da4B24DZ37wAws9uBhUA28LC7rw7H+jIw18z+EXgFeCjEHwJ+YmZ1QAPxAiMicnoaFh9AmrvvdarLyli/JwOKRRjR9OfdbPsG8I0u4guABV3ENxIfLdU5fgz4cCp5ioj0G+XngGXDnjXUDHs/q3f2zURIeoJbRCST5ORDyRjYt46xwwrZ2nCkTyZCUrEQEck0ZWdD/RvUVBQSc9hYf7jXv1LFQkQk05SfDQ0bqCmNj4hav7f3n+RWsRARyTRl50CsnersPWRnWZ9MsapiISKSacrPBiDvwAZGlw7qkxFRKhYiIpmmtCb+c986aoYN1m0oERHpwoAhMPiMNzu5N+8/Qkt7746IUrEQEclE5WfDvjcYO6yQjpized+RXv06FQsRkUxUdjbsW09NeSHQ+yOiVCxERDJR2dnQeogxA5rIMnq9k1vFQkQkE5WfA8CAxg2MKhnU68NnVSxERDJRWXz4bLzfovdHRKlYiIhkosIKyC+C+nWcVV7A5v1H6Ij13uwNKhYiIpnI7M0RUdVlBbS2x9jZeLTXvk7FQkQkU5XFi8WYMCJq477ee6FgSsXCzB4zs1fDZ7OZvZqw7S4zqzOzdWZ2TUJ8eojVmdmdCfFqM1sa4o+FmfQIs+09FuJLzWx0KjmLiPQbZTXQvIcxg+MP5G2q771O7pSKhbt/1N0nufsk4OfAEwBmNp74jHYTgOnAA2aWbWbZwP3AtcB44ObQFuCbwL3uPhY4ANwS4rcAB0L83tBORERKx8Z/tGxjcH4Om9L1yuI4MzPgI8DPQmgGMNfdW9x9E1BHfBa8KUCdu28Ms+zNBWaE/a8E5oX95wDXJxxrTlieB1wV2ouInN5KzgLAGjZSXV6QvrehEvwpsMfd14f1SmBbwvbtIdZdvBRodPf2TvG3HCtsbwrtRURObyXV8Z/7NzCmrKBXJ0FKWizMbLGZreriMyOh2c388aoiMmY2y8yWmdmy+vr6qNMREelduQNhSBU0bKC6rJCdTUd7bYrVnGQN3H3aybabWQ7wQWByQngHMDJhvSrE6Ca+HxhqZjnh6iGx/fFjbQ/fVRTad5XrbGA2QG1tbe8NOBYRSRelY2D/BqrPKsAdtuw/wjnDB/f41/TEbahpwFp3354Qmw/cFEYyVQM1wIvAS0BNGPmUR7wTfL67O/AccGPYfybwZMKxZoblG4FnQ3sRESk5Cxrit6EANu3rnRFRSa8s3oab6HQLyt1Xm9njwBqgHbjN3TsAzOx2YCGQDTzs7qvDbl8G5prZPwKvAA+F+EPAT8ysDmgI3yciIgClZ8HRA1QXtAKwoZf6LVIuFu7+qW7i3wC+0UV8AbCgi/hG4qOlOsePAR9ONU8RkX4pjIgqaN5CxZD8Xhs+qye4RUQyWWm8WLB/A9VlBSoWIiLSheLRYFnQsJEx5YUqFiIi0oWcfCiqerOTu+FwK41HWnv8a1QsREQyXclZb96Ggt55oaCKhYhIpiuND5+tKS/kT2vKyOqFNyL1xNBZERGJUslZcKyJUQOP8ZNbLumVr9CVhYhIpjs+IqphQ699hYqFiEimK/nj8NneomIhIpLpis98c/hsb1GxEBHJdNm58bfPNm7pta9QsRAR6Q+Kz4QDm3vt8CoWIiL9QfFoFQsREUmieDQ074HWI71yeBULEZH+oHh0/Gcv9VuoWIiI9AfFYT7uXroVpWIhItIfHL+ySMdiYWaTzGyJmb1qZsvMbEqIm5ndZ2Z1ZrbCzC5K2Gemma0Pn5kJ8clmtjLsc59Z/OUmZlZiZotC+0VmVpxKziIi/dKgEsgbnJ7FAvhn4B/cfRLwlbAOcC3xebdrgFnAgxD/xQ98FbiE+Kx4X0345f8g8LmE/aaH+J3AM+5eAzwT1kVEJJFZr46ISrVYODAkLBcBO8PyDOARj1sCDDWzEcA1wCJ3b3D3A8AiYHrYNsTdl7i7A48A1ycca05YnpMQFxGRRL34rEWqb539K2ChmX2beOG5LMQrgW0J7baH2Mni27uIA1S4+66wvBuo6C4ZM5tF/EqGUaNGvfOzERHJZMWjoW4xuMevNHpQ0isLM1tsZqu6+MwAbgW+6O4jgS8CD/Vodp2Eqw4/yfbZ7l7r7rXl5eW9mYqISPopHg3tx+LPW/SwpFcW7j6tu21m9ghwR1j9T+BHYXkHMDKhaVWI7QDe0yn+mxCv6qI9wB4zG+Huu8Ltqr3JchYROS0lDp8dPLxHD51qn8VO4N1h+UpgfVieD3wyjIqaCjSFW0kLgavNrDh0bF8NLAzbDprZ1DAK6pPAkwnHOj5qamZCXEREEpXVwLl/Fp+Xu4el2mfxOeB7ZpYDHCP0FwALgPcBdcAR4NMA7t5gZl8HXgrtvubuDWH588CPgYHAU+EDcA/wuJndAmwBPpJiziIi/VPxmfDRn/bKoS3eDdD/1NbW+rJly6JOQ0Qko5jZcnev7RzXE9wiIpKUioWIiCSlYiEiIkmpWIiISFIqFiIikpSKhYiIJKViISIiSfXb5yzMrJ74Q3ynogzY14PpREnnkn76y3mAziVdpXIuZ7r7CS/X67fFIhVmtqyrh1Iykc4l/fSX8wCdS7rqjXPRbSgREUlKxUJERJJSseja7KgT6EE6l/TTX84DdC7pqsfPRX0WIiKSlK4sREQkKRULERFJSsWiEzObbmbrzKzOzO6MOp9TYWYjzew5M1tjZqvN7I7ke6U3M8s2s1fM7JdR55IKMxtqZvPMbK2ZvW5ml0ad06kysy+Gv1+rzOxnZjYg6pzeLjN72Mz2mtmqhFiJmS0ys/XhZ3GUOb4d3ZzHt8LfrxVm9gszG9oT36VikcDMsoH7gWuB8cDNZjY+2qxOSTvw1+4+HpgK3Jah55HoDuD1qJPoAd8Dnnb3ccAFZOg5mVkl8AWg1t0nAtnATdFm9Y78GJjeKXYn8Iy71wDPhPV092NOPI9FwER3Px94A7irJ75IxeKtpgB17r7R3VuBucCMiHN6x9x9l7u/HJYPEf+FVBltVqfOzKqA9wM/ijqXVJhZEfAu4CEAd29198ZIk0pNDjAwTKs8CNgZcT5vm7v/FmjoFJ4BzAnLc4Dr+zKnU9HVebj7r929PawuAap64rtULN6qEtiWsL6dDP4lC2Bmo4ELgaURp5KKfwH+FohFnEeqqoF64N/DLbUfmVlB1EmdCnffAXwb2ArsAprc/dfRZpWyCnffFZZ3AxVRJtNDPgM81RMHUrHox8ysEPg58FfufjDqfE6FmV0H7HX35VHn0gNygIuAB939QuAwmXGr4wThfv4M4gXwDKDAzP482qx6jsefKcjo5wrM7G7it6Qf7YnjqVi81Q5gZMJ6VYhlHDPLJV4oHnX3J6LOJwWXAx8ws83EbwteaWY/jTalU7Yd2O7ux6/y5hEvHploGrDJ3evdvQ14Args4pxStcfMRgCEn3sjzueUmdmngOuAj3sPPUynYvFWLwE1ZlZtZnnEO+zmR5zTO2ZmRvy++Ovu/t2o80mFu9/l7lXuPpr4f49n3T0j/wXr7ruBbWZ2TghdBayJMKVUbAWmmtmg8PftKjK0sz7BfGBmWJ4JPBlhLqfMzKYTv237AXc/0lPHVbFIEDqFbgcWEv+L/7i7r442q1NyOfAJ4v8KfzV83hd1UgLAXwKPmtkKYBLw/6JN59SEq6N5wMvASuK/SzLmdRlm9jPgD8A5ZrbdzG4B7gHea2briV853RNljm9HN+fxA2AwsCj8v//DHvkuve5DRESS0ZWFiIgkpWIhIiJJqViIiEhSKhYiIpKUioWIiCSlYiEiIkmpWIiISFL/H183XZDeEOkrAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], "source": [ - "# Do your plotting and your own testing here ...\n", + "def ODEF(x, y):\n", + " return y - x**2 + 1.0\n", "\n", - "# YOUR CODE HERE\n", - "raise NotImplementedError()" + "x = np.linspace(0, 12, 500)\n", + "y0 = .5\n", + "\n", + "eta = integrator(x, y0, ODEF, phi_euler)\n", + "y = (x + 1)**2 - 0.5*np.exp(x)\n", + "\n", + "plt.plot(x, eta)\n", + "plt.plot(x, y)" ] }, {