634 lines
31 KiB
Plaintext
634 lines
31 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "4ec40081b048ce2f34f3f4fedbb0be10",
|
|
"grade": false,
|
|
"grade_id": "cell-98f724ece1aacb67",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"# CDS: Numerical Methods Assignments\n",
|
|
"\n",
|
|
"- See lecture notes and documentation on Brightspace for Python and Jupyter basics. If you are stuck, try to google or get in touch via Discord.\n",
|
|
"\n",
|
|
"- Solutions must be submitted via the Jupyter Hub.\n",
|
|
"\n",
|
|
"- Make sure you fill in any place that says `YOUR CODE HERE` or \"YOUR ANSWER HERE\".\n",
|
|
"\n",
|
|
"## Submission\n",
|
|
"\n",
|
|
"1. Name all team members in the the cell below\n",
|
|
"2. make sure everything runs as expected\n",
|
|
"3. **restart the kernel** (in the menubar, select Kernel$\\rightarrow$Restart)\n",
|
|
"4. **run all cells** (in the menubar, select Cell$\\rightarrow$Run All)\n",
|
|
"5. Check all outputs (Out[\\*]) for errors and **resolve them if necessary**\n",
|
|
"6. submit your solutions **in time (before the deadline)**"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "raw",
|
|
"metadata": {},
|
|
"source": [
|
|
"team_members = \"Koen Vendrig, Kees van Kempen\""
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "7317367b26066eae942edeea4e7d4280",
|
|
"grade": false,
|
|
"grade_id": "cell-e9dc85dd9ad77a5e",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"## Ordinary Differential Equations - Initial Value Problems\n",
|
|
"\n",
|
|
"In the following you will implement your own easy-to-extend ordinary differential equation (ODE) solver, which can be used to solve first order ODE systems of the form\n",
|
|
"\n",
|
|
"\\begin{align*}\n",
|
|
" \\vec{y}'(x) = \\vec{f}(x, \\vec{y}),\n",
|
|
"\\end{align*}\n",
|
|
"\n",
|
|
"for $x = x_0, x_1, \\dots, x_n$ with $x_i = i h$, step size $h$, and an initial condition of the form $\\vec{y}(x=0)$.\n",
|
|
"The solver will be capable of using single-step as well as multi-step approaches. "
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 1,
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "910692c2f4b93f251a3a128caf2b0c37",
|
|
"grade": true,
|
|
"grade_id": "cell-acb87410eaef9fe1",
|
|
"locked": false,
|
|
"points": 0,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"import numpy as np\n",
|
|
"from matplotlib import pyplot as plt"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "79e708efbac889db28e537bcf01e26c6",
|
|
"grade": false,
|
|
"grade_id": "cell-a27c54a734bf2232",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"### Task 1\n",
|
|
"\n",
|
|
"Implement the Euler method in a general Python function $\\text{integrator(x, y0, f, phi)}$, which takes as input a one-dimensional array $x$ of size $n+1$, the initial condition $y_0$, the *callable* function $f(x, y)$, and the integration scheme $\\Phi(x, y, f, i)$. It should return the approximated function $\\tilde{y}(x)$. \n",
|
|
"\n",
|
|
"The integration scheme $\\text{phi}$ is supposed to be a *callable* function as returned from another Python function $\\text{phi_euler(x, y, f, i)}$, where $i$ is the step number. In this way we will be able to easily extend the ODE solver to different methods."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 2,
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "8396264309b023abf89a62200fa55379",
|
|
"grade": true,
|
|
"grade_id": "cell-1f9655333f2dc3cb",
|
|
"locked": false,
|
|
"points": 3,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"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 integration scheme provided by phi.\n",
|
|
"\n",
|
|
" Args:\n",
|
|
" x: size n + 1 numerical array\n",
|
|
" y0: an initial value to the function f\n",
|
|
" 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, 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",
|
|
" eta = np.zeros(len(x))\n",
|
|
" eta[0] = y0\n",
|
|
" \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 eta"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 3,
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "35d6f653b74f2c8e2ff509890754df72",
|
|
"grade": true,
|
|
"grade_id": "cell-04ac3a0537092706",
|
|
"locked": false,
|
|
"points": 3,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"def phi_euler(x, y, f, i):\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",
|
|
"\n",
|
|
" Args:\n",
|
|
" x: a size n + 1 numerical array\n",
|
|
" y: a size n + 1 numerical array with y[0] the initial value corresponding to x[0]\n",
|
|
" f: a callable function with signature (x, y), with x and y the current state\n",
|
|
" of the system\n",
|
|
" i: 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",
|
|
"\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",
|
|
" return f(x, y)\n",
|
|
" "
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "a300f75f62faab4732ebcac7bde5e18e",
|
|
"grade": false,
|
|
"grade_id": "cell-047f2dfab5489a88",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"### Task 2 \n",
|
|
"\n",
|
|
"Debug your implementation by applying it to the following ODE\n",
|
|
"\n",
|
|
"\\begin{align*}\n",
|
|
" \\vec{y}'(x) = y - x^2 + 1.0,\n",
|
|
"\\end{align*}\n",
|
|
" \n",
|
|
"with initial condition $y(x=0) = 0.5$. To this end, define the right-hand side of the ODE as a Python function $\\text{ODEF(x, y)}$, which in this case returns $f(x,y) = y - x^2 + 1.0$. You can then hand over the function $\\text{ODEF}$ to the argument $\\text{f}$ of your $\\text{integrator}$ function.\n",
|
|
"\n",
|
|
"Plot the solution you found with the Euler method together with the exact solution $y(x) = (x+1)^2 - 0.5 e^x$.\n",
|
|
"\n",
|
|
"Then implement a unit test for your $\\text{integrator}$ function using this system."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 4,
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "692bfda705027a7fa4452a7e03c8fa7d",
|
|
"grade": true,
|
|
"grade_id": "cell-8b8cb282dfe95a03",
|
|
"locked": false,
|
|
"points": 4,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"[<matplotlib.lines.Line2D at 0x1554d71a0610>]"
|
|
]
|
|
},
|
|
"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": [
|
|
"<Figure size 432x288 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"def ODEF(x, y):\n",
|
|
" return y - x**2 + 1.0\n",
|
|
"\n",
|
|
"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)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "aea0c6713e9bd3871c2e93e65e93d903",
|
|
"grade": true,
|
|
"grade_id": "cell-715c304f22415d45",
|
|
"locked": false,
|
|
"points": 3,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"def test_integrator():\n",
|
|
" # YOUR CODE HERE\n",
|
|
" raise NotImplementedError()\n",
|
|
" \n",
|
|
"test_integrator()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "8048460a3a4d5e94b9e9f086d90ff0d0",
|
|
"grade": false,
|
|
"grade_id": "cell-99399435c164c96d",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"### Task 3\n",
|
|
"\n",
|
|
"Extend the set of possible single-step integration schemes to the modified Euler (Collatz), Heun, and 4th-order Runge-Kutta approaches by implementing the functions $\\text{phi_euler_modified(x, y, f, i)}$, $\\text{phi_heun(x, y, f, i)}$, and $\\text{phi_rk4(x, y, f, i)}$. These can be used in your $\\text{integrator}$ function instead of $\\text{phi_euler}$."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "dc82ac823cbd262952a6990776e26d09",
|
|
"grade": true,
|
|
"grade_id": "cell-ea865d2efbc574d6",
|
|
"locked": false,
|
|
"points": 9,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"def phi_euler_modified(x, y, f, i):\n",
|
|
" # YOUR CODE HERE\n",
|
|
" raise NotImplementedError()\n",
|
|
" \n",
|
|
"def phi_heun(x, y, f, i):\n",
|
|
" # YOUR CODE HERE\n",
|
|
" raise NotImplementedError()\n",
|
|
" \n",
|
|
"def phi_rk4(x, y, f, i):\n",
|
|
" # YOUR CODE HERE\n",
|
|
" raise NotImplementedError()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "996fa4ae58be6c15aed465d8d7e5618e",
|
|
"grade": false,
|
|
"grade_id": "cell-a783908a0db32a24",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"### Task 4\n",
|
|
"\n",
|
|
"Add the possibility to also handle the following multi-step integration schemes:\n",
|
|
"\n",
|
|
"\\begin{align*}\n",
|
|
" \\Phi_{AB3}(x, y, f, i) = \n",
|
|
" \\frac{1}{12} \\left[ \n",
|
|
" 23 f(x_i, y_i) \n",
|
|
" - 16 f(x_{i-1}, y_{i-1})\n",
|
|
" + 5 f(x_{i-2}, y_{i-2})\n",
|
|
" \\right]\n",
|
|
"\\end{align*}\n",
|
|
"\n",
|
|
"and\n",
|
|
"\n",
|
|
"\\begin{align*}\n",
|
|
"\\Phi_{AB4}(x, y, f, i) = \n",
|
|
" \\frac{1}{24} \\left[ \n",
|
|
" 55 f(x_i, y_i) \n",
|
|
" - 59 f(x_{i-1}, y_{i-1})\n",
|
|
" + 37 f(x_{i-2}, y_{i-2})\n",
|
|
" - 9 f(x_{i-3}, y_{i-3})\n",
|
|
" \\right]\n",
|
|
"\\end{align*}\n",
|
|
" \n",
|
|
"In these cases the initial condition $\\text{y0}$ must be an array consisting of several initial values corresponding to $y_0$, $y_1$, $y_2$ (and $y_3$). Use the Runga-Kutta method to calculate all of the initial values before you start the AB3 and AB4 integrations."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "d301f6514b2bc50a77d929518f6928df",
|
|
"grade": true,
|
|
"grade_id": "cell-0bd05d93759bd652",
|
|
"locked": false,
|
|
"points": 6,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"def phi_ab3(x, y, f, i):\n",
|
|
" # YOUR CODE HERE\n",
|
|
" raise NotImplementedError()\n",
|
|
"\n",
|
|
"def phi_ab4(x, y, f, i):\n",
|
|
" # YOUR CODE HERE\n",
|
|
" raise NotImplementedError()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "c621385b0e71eb6fd2e4d6041268bd15",
|
|
"grade": false,
|
|
"grade_id": "cell-3fac52fc46f0a497",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"### Task 5\n",
|
|
"\n",
|
|
"Plot the absolute errors $\\delta(x) = |\\tilde{y}(x) - y(x)|$ with a $y$-log scale for all approaches with $0 \\leq x \\leq 2$ and a step size of $h=0.02$ for the ODE from task 2."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "22f375f78eeb46befb192c08db0f48e9",
|
|
"grade": true,
|
|
"grade_id": "cell-e61d969c23914a5f",
|
|
"locked": false,
|
|
"points": 4,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"# YOUR CODE HERE\n",
|
|
"raise NotImplementedError()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "9fd7c159424b1d536c11d53764cce54c",
|
|
"grade": false,
|
|
"grade_id": "cell-3c99ac2c49457d60",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"### Task 6\n",
|
|
"\n",
|
|
"Study the accuracies of all approaches as a function of the step size $h$. To this end use your implementations to solve\n",
|
|
"\n",
|
|
"\\begin{align*}\n",
|
|
" \\vec{y}'(x) = y\n",
|
|
"\\end{align*}\n",
|
|
" \n",
|
|
"with $y(x=0) = 1.0$ for $0 \\leq x \\leq 2$. \n",
|
|
"\n",
|
|
"Plot $\\delta(2) = |\\tilde{y}(2) - y(2)|$ as a function of $h$ for each integration scheme."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "0060c86da7630a17a6da4769bb0c728f",
|
|
"grade": true,
|
|
"grade_id": "cell-552ef9b066bc9b1c",
|
|
"locked": false,
|
|
"points": 4,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"# YOUR CODE HERE\n",
|
|
"raise NotImplementedError()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "0e3fcb92129c0b3ef4e2b8f99a38ddef",
|
|
"grade": false,
|
|
"grade_id": "cell-3070782f40a2599f",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"### Task 7\n",
|
|
"\n",
|
|
"Apply the Euler and the Runga-Kutta methods to solve the pendulum problem given by\n",
|
|
"\n",
|
|
"\\begin{align*}\n",
|
|
" \\vec{y}'(x) = \\vec{f}(x, \\vec{y})\n",
|
|
" \\quad \\leftrightarrow \\quad \n",
|
|
" \\begin{pmatrix}\n",
|
|
" y_0'(x) \\\\\n",
|
|
" y_1'(x)\n",
|
|
" \\end{pmatrix}\n",
|
|
" =\n",
|
|
" \\begin{pmatrix}\n",
|
|
" y_1(x) \\\\\n",
|
|
" -\\operatorname{sin}\\left[y_0(x) \\right]\n",
|
|
" \\end{pmatrix}\n",
|
|
"\\end{align*}\n",
|
|
" \n",
|
|
"To this end the quantities $\\text{x}$ and $\\text{y0}$ in your $\\text{integrator}$ function must become vectors. Depending on your implementation in task 1 you might have to define a new $\\text{integrator}$ function that can handle this type of input.\n",
|
|
"\n",
|
|
"Plot $y_0(x)$ for several oscillation periods. Does the Euler method behave physically?"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "9cb9ca1fb8e9c0722a7e6b5ec0ff0c2a",
|
|
"grade": true,
|
|
"grade_id": "cell-3a89944b55976666",
|
|
"locked": false,
|
|
"points": 4,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"# YOUR CODE HERE\n",
|
|
"raise NotImplementedError()"
|
|
]
|
|
}
|
|
],
|
|
"metadata": {
|
|
"kernelspec": {
|
|
"display_name": "Python 3",
|
|
"language": "python",
|
|
"name": "python3"
|
|
},
|
|
"language_info": {
|
|
"codemirror_mode": {
|
|
"name": "ipython",
|
|
"version": 3
|
|
},
|
|
"file_extension": ".py",
|
|
"mimetype": "text/x-python",
|
|
"name": "python",
|
|
"nbconvert_exporter": "python",
|
|
"pygments_lexer": "ipython3",
|
|
"version": "3.8.10"
|
|
}
|
|
},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 4
|
|
}
|