w2/discrete_oscillator.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ada769df-0827-4614-a16a-5cb7ff34b0db",
   "metadata": {},
   "source": [
    "# Tutorial 3\n",
    "Before we start, let's import our libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "319505e0-e5f9-4ef9-8a76-3c683d64b89b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09c32f0c-1395-4678-b057-f7d13577b61a",
   "metadata": {},
   "source": [
    "## Finite differencing and the explicit Euler method\n",
    "\n",
    "Today we want to solve the oscillator problem (again). But this time, we do not rely on a fancy ODE solver like the `odeint`. Instead, we want to solve the problem with a [finite difference](https://en.wikipedia.org/wiki/Finite_difference_method) method.\n",
    "\n",
    "The finite difference method is a numerical method to solve differential equations. The broad idea behind the method is to approximate the differential with a discrete difference. For example, in our oscillator problem, we have the time derivative of $x$. This can be *discretised* as \n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{dx}{dt} \\approx \\frac{x(t+\\Delta t) - x(t)}{\\Delta t}\\tag{1}.\n",
    "\\end{equation}\n",
    "\n",
    "The process of approximating a continuous quantity with its discrete form is called *discretisation*.\n",
    "\n",
    "Say now we have an ODE\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{dx}{dt} = F(t,x(t))\\tag{2},\n",
    "\\label{eq:differential}\n",
    "\\end{equation}\n",
    "\n",
    "discretising this yields\n",
    "\\begin{equation}\n",
    "x(t+\\Delta t) = x(t) + \\Delta t \\, F(t,x(t)) \\tag{3} \\label{eq:discretised}\n",
    "\\end{equation}\n",
    "\n",
    "Generally, we write equation \\eqref{eq:discretised} as \n",
    "\\begin{equation}\n",
    "x^{n+1} = x^n + \\Delta t \\, F(t^n,x^n), \\tag{4}\\label{eq:euler_method}\n",
    "\\end{equation}\n",
    "\n",
    "where $n$ indexes the time-step at time $t$ and $n+1$ indexes the time-step at $t+\\Delta t$. Note that the notion of a \"time-step\" now makes sense, as we are integrating the problem discretely in time.\n",
    "\n",
    "Given a time-step size of $\\Delta t$, if we have the solution $(t^n,x^n)$ at time $n$, then we can evaluate the right-hand side of \\eqref{eq:euler_method}, and this gives us the solution of $x$ at the new time $n+1$.\n",
    "\n",
    "The method of numerically integrating an ODE with \\eqref{eq:euler_method} is also known as the [explicit Euler method](https://en.wikipedia.org/wiki/Euler_method).\n",
    "\n",
    "Now that we have some idea on a finite difference method, let's apply it..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f600fad-cdd7-4075-914e-24991adf7f37",
   "metadata": {},
   "source": [
    "## Setting up the harmonic oscillator problem\n",
    "\n",
    "We want to use the simple harmonic oscillator problem without damping and without forcing as an exercise to get an intuition of the finite difference numerical methods:\n",
    "$$m \\ddot{x} = -c x, \\qquad m,c > 0. \\tag{5}$$\n",
    "\n",
    "The analytical solution to this free oscillation scenario is on page 20, equation (75) of the lecture notes,\n",
    "\n",
    "$$x(t) = a \\cos(\\omega t) + b \\sin(\\omega t), \\tag{6}$$\n",
    "\n",
    "where $\\omega=\\sqrt{c/m}$.\n",
    "\n",
    "We can make our lives much easier by choosing the constants $c,m=1$, the initial position $x(t=0)=0$, and the initial velocity $\\dot{x}(t=0)=1$. This gives us the following analytical solutions:\n",
    "$$\\begin{align}\n",
    "x        &= \\sin(t), \\tag{7a}\\\\\n",
    "\\dot{x}  &= \\cos(t), \\tag{7b}\\\\\n",
    "\\ddot{x} &= -x. \\tag{7c}\n",
    "\\end{align}$$\n",
    "This is equivalent to\n",
    "$$\\begin{align}\n",
    "\\frac{dx}{dt}        &= \\dot{x}, \\tag{8a}\\\\\n",
    "\\frac{d\\dot{x}}{dt}  &= -x. \\tag{8b}\n",
    "\\end{align}$$\n",
    "\n",
    "Discretising this with the explicit Euler method,\n",
    "$$\\begin{align}\n",
    "x^{n+1}       &= x^n + \\Delta t \\, \\dot{x}^n, \\tag{9a}\\\\\n",
    "\\dot{x}^{n+1} &= \\dot{x}^n - \\Delta t \\, x^n. \\tag{9b}\n",
    "\\label{eq:ee_discretised}\n",
    "\\end{align}$$\n",
    "\n",
    "Please go through the working above on your own to verify that it makes sense!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "24bdebe0-9f71-4659-973e-2097c12b791a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Some helper functions for plots\n",
    "def xt_plot(t, x, xth=None, title=\"\", ylabel=\"x\"):\n",
    "    plt.plot(t, x, label='num. x')\n",
    "    if xth is not None:\n",
    "        plt.plot(t, xth, '--', label='theo. x')\n",
    "    plt.legend(loc='best')\n",
    "    plt.xlabel('t')\n",
    "    plt.ylabel(ylabel)\n",
    "    if len(title) > 0:\n",
    "        plt.title(title)\n",
    "    plt.grid()\n",
    "    plt.show()\n",
    "\n",
    "def vx_plot(x, v, xth, vth, title=\"\"):\n",
    "    plt.plot(x,v, label=\"num. sol.\")\n",
    "    plt.plot(xth,vth, '--', label=\"theo. sol.\")\n",
    "    plt.legend()\n",
    "    plt.xlabel(\"x\")\n",
    "    plt.ylabel(\"v\")\n",
    "    if len(title) > 0:\n",
    "        plt.title(title)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ec7d32a-b82e-4c32-800d-406b72f8cf07",
   "metadata": {},
   "source": [
    "## Exercise 1: The explicit Euler method \n",
    "1. Implement equation ([9](#mjx-eqn-eq:ee_discretised)). How would you do it?\n",
    "2. Compare the result for various step-sizes, e.g. `dt=0.1` and `dt=0.01`.\n",
    "3. Plot $\\dot{x}$ against $x$ for the numerical and the analytical solution. What do you see?\n",
    "4. We know that the potential energy in a spring is given by $PE=\\frac{1}{2}kx^2$ while the kinetic energy of the spring is $KE=\\frac{1}{2}mv^2$. Compute and compare the energy of the numerical solution with the theoretical energy. What happens with the energy?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40b00eb7-5e8c-41bb-b5e7-165eaa4a4365",
   "metadata": {},
   "source": [
    "### Hints:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95cad889-d746-4c5a-9590-af0763d51a75",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "source": [
    "1. We will need a *for-loop*. A short introduction to the features of the for-loop in Python is provided in the supplementary materials.\n",
    "2. What are we looping over?\n",
    "3. What do we have to do in each step of the loop?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "11b3fc8b-eeb5-4649-a55c-5ac7696d68a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here are some helper code to get you started.\n",
    "\n",
    "# Let's specify the number of time-steps\n",
    "Nt = 10000\n",
    "# and a time-step size\n",
    "dt = 0.01\n",
    "\n",
    "# What are numpy.arange and numpy.zeros?\n",
    "# Why do we need to initialise these arrays?\n",
    "t = np.arange(Nt+1)*dt\n",
    "x = np.zeros((Nt+1))\n",
    "v = np.zeros((Nt+1))\n",
    "\n",
    "# Since we know the analytical solution,\n",
    "# let's compare our numerical solution against them.\n",
    "xth = np.sin(t)\n",
    "vth = np.cos(t)\n",
    "\n",
    "# Fill in the rest of the code here:\n",
    "# ...\n",
    "\n",
    "# xt_plot(t,x,xth, title=\"explicit Euler\")\n",
    "# vx_plot(x,v,xth,vth, title=\"explicit Euler\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67d02d8f-e732-4e45-86c8-98e56d0bf80c",
   "metadata": {},
   "source": [
    "## Exercise 2: Euler-A and Euler-B methods\n",
    "Notice that our discretisation of equation ([2](#mjx-eqn-eq:differential)) is not unique, and there are actually many ways to discretise a differential. For example, we can write equation ([4](#mjx-eqn-eq:euler_method)) as\n",
    "\n",
    "$$x^{n+1} = x^n - \\Delta t \\, F(t^{n+1},x^{n+1}), \\tag{10} \\label{eq:implicit_euler} $$\n",
    "\n",
    "which would give us the [implicit Euler method](https://en.wikipedia.org/wiki/Backward_Euler_method).\n",
    "\n",
    "For the system of equations we are solving in ([9](#mjx-eqn-eq:ee_discretised)), we can apply the [semi-implicit Euler method](https://en.wikipedia.org/wiki/Semi-implicit_Euler_method). This is given by the following updates:\n",
    "\n",
    "a. For the Euler-A method:\n",
    "$$\\begin{align}\n",
    "x^{n+1}       &= x^n + \\Delta t \\, \\dot{x}^n, \\tag{11a}\\\\\n",
    "\\dot{x}^{n+1} &= \\dot{x}^n - \\Delta t \\, x^{n+1}. \\tag{11b}\n",
    "\\label{eq:euler_a}\n",
    "\\end{align}$$\n",
    "b. For the Euler-B method:\n",
    "$$\\begin{align}\n",
    "\\dot{x}^{n+1} &= \\dot{x}^n - \\Delta t \\, x^{n}, \\tag{12a}\\\\\n",
    "x^{n+1}       &= x^n + \\Delta t \\, \\dot{x}^{n+1}. \\tag{12b}\n",
    "\\label{eq:euler_b}\n",
    "\\end{align}$$\n",
    "\n",
    "Now on to the tasks:\n",
    "1. Implement equations ([11](#mjx-eqn-eq:euler_a)) and ([12](#mjx-eqn-eq:euler_b)).\n",
    "2. Again, try around with different number of time-steps and time-step sizes.\n",
    "3. Plot $\\dot{x}$ against $x$ for the numerical and analytical solutions. What do you observe?\n",
    "4. Can you explain your observation?\n",
    "5. Again, let's compare the energy of the numerical solution with the theoretical energy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "d65c9086-f241-48d0-9118-875b37727c90",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fill in the code here:\n",
    "# ...\n",
    "\n",
    "# xt_plot(t, xb, xth, title=\"Euler B\")\n",
    "# vx_plot(xb,vb,xth,vth, title=\"Euler B\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1a4f16a-f625-4359-a538-d32355fecd0d",
   "metadata": {},
   "source": [
    "## Exercise 3: The Störmer-Verlet method\n",
    "The explicit Euler and the semi-implicit Euler methods are [*first-order methods*](https://en.wikipedia.org/wiki/Truncation_error_(numerical_integration)#Local_truncation_error). We now turn our attention to the [Störmer-Verlet method](https://en.wikipedia.org/wiki/Verlet_integration), which is a second-order numerical integrator.\n",
    "\n",
    "The discretised equations are:\n",
    "$$\\begin{align}\n",
    "\\dot{x}^{n+1/2}    &= \\dot{x}^n       - \\frac{\\Delta t}{2} \\, x^n,     \\label{eq:sv1} \\tag{13a}\\\\\n",
    "           x^{n+1} &= x^n             + \\Delta t \\, \\dot{x}^{n+1/2},   \\label{eq:sv2} \\tag{13b}\\\\\n",
    "\\dot{x}^{n+1}      &= \\dot{x}^{n+1/2} - \\frac{\\Delta t}{2} \\, x^{n+1}. \\label{eq:sv3} \\tag{13c}\\\\\n",
    "\\end{align}$$\n",
    "\n",
    "Notice that ([13a](#mjx-eqn-eq:sv1)) is akin to an explicit Euler update, ([13b](#mjx-eqn-eq:sv2)) is a [midpoint method](https://en.wikipedia.org/wiki/Midpoint_method), and ([13c](#mjx-eqn-eq:sv3)) is an implicit update.\n",
    "\n",
    "Now, for the tasks:\n",
    "1. Implement equation ([13](#mjx-eqn-eq:sv3)).\n",
    "2. Plot $\\dot{x}$ against $x$ for the numerical and analytical solutions. What do you observe?\n",
    "3. Plot $v$ and $x$ and the energy. Are these results in line with your expectations?\n",
    "4. Plot the intermediate time steps in ([13a](#mjx-eqn-eq:sv1)). What do you observe?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "50f9f65b-56fe-43f8-abb5-9292accefb27",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fill in the code here:\n",
    "# ...\n",
    "    \n",
    "# xt_plot(t,x,xth, title=\"Störmer-Verlet\")\n",
    "# vx_plot(v,x,xth,vth, title=\"Stömer-Verlet\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f5a9ed8-3b39-4e15-8d9b-3f0776d14c02",
   "metadata": {},
   "source": [
    "## Further reading\n",
    "\n",
    "The implicit and explicit Euler methods are part of a family of numerical integrators called the [*Runge-Kutta*](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods) methods. The generalised Runge-Kutta method allows us to build higher-order integrators, and Wikipedia has such [a list](https://en.wikipedia.org/wiki/List_of_Runge%E2%80%93Kutta_methods). You will come across the Runge-Kutta methods frequently, especially as time integrators. We will go more in detail when we move on to numerical solution of PDEs later in this tutorial."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f7281e0-ac5a-408b-9ff1-e40bfc4f2201",
   "metadata": {},
   "source": [
    "## References\n",
    "[1] Benedict Leimkuhler and Sebastian Reich. Simulating Hamiltonian Dynamics. Cambridge Monographs on Applied and Computational Mathematics. Cambridge University Press, 2005.\n",
    "\n",
    "[2] Ernst Hairer, Christian Lubich, and Gerhard Wanner. Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Diﬀerential Equations; 2nd ed. Springer, Dordrecht, 2006."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a26ff85-eb87-4cb1-93a4-a754844859d5",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.7.3 64-bit ('anaconda3': virtualenv)",
   "language": "python",
   "name": "python37364bitanaconda3virtualenv7a28dc8db0264e168ad4f93d8f9f620c"
  },
  "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}