From 4fd1c63aa0d29bc25a3f8a748f988adf5e4e7c19 Mon Sep 17 00:00:00 2001 From: Kees van Kempen Date: Tue, 8 Mar 2022 10:26:52 +0100 Subject: [PATCH] Fetch week 5 --- .../9 Ordinary Differential Equations.ipynb | 564 ++++++++++++++++++ 1 file changed, 564 insertions(+) create mode 100644 Week 5/9 Ordinary Differential Equations.ipynb diff --git a/Week 5/9 Ordinary Differential Equations.ipynb b/Week 5/9 Ordinary Differential Equations.ipynb new file mode 100644 index 0000000..9b681f3 --- /dev/null +++ b/Week 5/9 Ordinary Differential Equations.ipynb @@ -0,0 +1,564 @@ +{ + "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 = \"\"" + ] + }, + { + "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": null, + "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 packages here ...\n", + "\n", + "# YOUR CODE HERE\n", + "raise NotImplementedError()" + ] + }, + { + "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": null, + "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", + " # YOUR CODE HERE\n", + " raise NotImplementedError()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "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", + " # YOUR CODE HERE\n", + " raise NotImplementedError()" + ] + }, + { + "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": null, + "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": [], + "source": [ + "# Do your plotting and your own testing here ...\n", + "\n", + "# YOUR CODE HERE\n", + "raise NotImplementedError()" + ] + }, + { + "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 (ipykernel)", + "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 +}