{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "\n", "\n", "\n", "# Numerical solution of ordinary differential equations: Euler's and Heun's method\n", "\n", " \n", "**Anne Kværnø, André Massing)**\n", "\n", "Date: **Apr 2, 2020**\n", "\n", "If you want to have a nicer theme for your jupyter notebook,\n", "download the [cascade stylesheet file tma4125.css](https://www.math.ntnu.no/emner/TMA4125/2020v/part_II/notebooks/tma4125.css)\n", "and execute the next cell:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [ { "data": { "text/html": [ " \n", "\n" ], "text/plain": [ "" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "from IPython.core.display import HTML\n", "from numpy import pi\n", "from numpy.linalg import norm, solve\n", "\n", "\n", "def css_styling():\n", " try:\n", " with open(\"tma4125.css\") as f:\n", " styles = f.read()\n", " return HTML(styles)\n", " except FileNotFoundError:\n", " pass # Do nothing\n", "\n", "\n", "# Comment out next line and execute this cell to restore the default notebook style\n", "css_styling()" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "As always, we start by calling the necessary modules:\n", "And of course we want to import the required modules." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "%matplotlib inline\n", "\n", "\n", "# Use a funny plotting style\n", "plt.xkcd()\n", "newparams = {\n", " \"figure.figsize\": (6.0, 6.0),\n", " \"axes.grid\": True,\n", " \"lines.markersize\": 8,\n", " \"lines.linewidth\": 2,\n", " \"font.size\": 14,\n", "}\n", "plt.rcParams.update(newparams)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Introduction: Whetting your appetite\n", "The topic of this note is the numerical solution of systems of\n", "ordinary differential equations (ODEs). This has been discussed in\n", "previous courses, see for instance the web page\n", "[Differensialligninger](https://wiki.math.ntnu.no/tma4100/tema/differentialequations)\n", "from Mathematics 1, as well as in Part 1 of this course, where the\n", "Laplace transform was introduced as a tool to solve ODEs analytically.\n", "\n", "Before we present the first numerical methods to solve ODEs, we quickly\n", "look at a number of examples which hopefully will will serve as test examples\n", "throughout this topic." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Scalar first order ODEs\n", "A scalar, first-order ODE is an equation on the form\n", "\n", "$$\n", "y'(t) = f(t,y(t)), \\qquad y(t_0)=y_0,\n", "$$\n", "\n", "where $y'(t)=\\frac{dy}{dx}$.\n", "The *inital condition* $y(t_0)=y_0$ is required for a unique\n", "solution. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "**Notice.**\n", "\n", "It is common to use the term *initial value problem (IVP)* for an ODE\n", "for which the inital value $y(t_0)=y_0$ is given, and we only are\n", "interested in the solution for $x>t_0$. In these lecture notes, only\n", "initial value problems are considered." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Example 1: Population growth and decay processes\n", "\n", "
\n", "One of the simplest possible IVP is given by\n", "\n", "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "y'(t) = \\lambda y(t), \\quad y(t_0)=y_0.\n", "\\label{ode:eq:exponential} \\tag{1}\n", "\\end{equation}\n", "$$\n", "\n", "For $\\lambda > 0$ this equation can present a simple model for the growth of\n", "some population, e.g., cells, humans, animals, with unlimited resources\n", "(food, space etc.). The constant $\\lambda$ then corresponds to the\n", "*growth rate* of the population.\n", "\n", "Negative $\\lambda < 0$\n", "typically appear in decaying processes, e.g., the decay of a radioactive\n", "isotopes, where $\\lambda$ is then simply called the *decay constant*.\n", "\n", "The analytical solution to ([1](#ode:eq:exponential)) is\n", "\n", "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "y(t) = y_0 e^{\\lambda(t-t_0)}\n", "\\label{ode:eq:exponential_sol} \\tag{2}\n", "\\end{equation}\n", "$$\n", "\n", "and will serve us at several occasions as reference solution to assess\n", "the accuracy of the numerical methods to be introduced." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Example 2: Time-dependent coefficients\n", "\n", "Given the ODE\n", "\n", "$$\n", "y'(t) = -2ty(t), \\quad y(0) = y_0.\n", "$$\n", "\n", "for some given initial value $y_0$.\n", "The general solution of the ODE is\n", "\n", "$$\n", "y(t) = C e^{-t^2},\n", "$$\n", "\n", "where $C$ is a constant. To determine the constant $C$,\n", "we use the initial condition $y(0) = y_0$\n", "yielding the solution\n", "\n", "$$\n", "y(t) = y_0 e^{-t^2}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Systems of ODEs\n", "A system of $m$ ODEs are given by\n", "\n", "$$\n", "\\begin{align*}\n", "y_1' &= f_1(t,y_1,y_2,\\dotsc,y_m), & y_1(t_0) &= y_{1,0} \\\\ \n", "y_2' &= f_2(t,y_1,y_2,\\dotsc,y_m), & y_2(t_0) &= y_{2,0} \\\\ \n", " & \\vdots & &\\vdots \\\\ \n", "y_m' &= f_m(t,y_1,y_2,\\dotsc,y_m), & y_m(t_0) &= y_{m,0} \\\\ \n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Or more compactly by\n", "\n", "$$\n", "\\mathbf{y}'(t) = \\mathbf{f}(t, \\mathbf{y}(t)), \\qquad \\mathbf{y}(t_0) = \\mathbf{y}_0\n", "$$\n", "\n", "where we use boldface to denote vectors in $\\mathbb{R}^m$,\n", "\n", "$$\n", "\\mathbf{y}(t) =\n", "\\left(\n", "\\begin{array}{c}\n", "y_1(t) \n", "\\\\ y_2(t) \n", "\\\\ \\vdots \n", "\\\\ y_m(t)\n", "\\end{array}\n", "\\right),\n", "\\qquad\n", "\\mathbf{f}(t,\\mathbf{y}) =\n", "\\left(\n", "\\begin{array}{c}\n", "f_1(t,y_1,y_2,\\dotsc,y_m), \n", "\\\\ f_2(t,y_1,y_2,\\dotsc,y_m), \n", "\\\\ \\vdots \n", "\\\\ f_m(t,y_1,y_2,\\dotsc,y_m),\n", "\\end{array}\n", "\\right),\n", "\\qquad\n", "\\mathbf{y}_0 =\n", "\\left(\n", "\\begin{array}{c}\n", "y_{1,0} \n", "\\\\ y_{2,0} \n", "\\\\ \\vdots \n", "\\\\ y_{m,0}\n", "\\end{array}\n", "\\right).\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Example 3: Lotka-Volterra equation\n", "\n", "The [Lotka-Volterra equation](https://en.wikipedia.org/wiki/Lotka-Volterra_equations) is\n", "a system of two ODEs describing the interaction between preys and\n", "predators over time. The system is given by\n", "\n", "$$\n", "\\begin{align*}\n", "y'(t) &= \\alpha y(t) - \\beta y(t) z(t) \\\\ \n", "z'(t) &= \\delta y(t)z(t) - \\gamma z(t)\n", "\\end{align*}\n", "$$\n", "\n", "where $x$ denotes time, $y(t)$ describes the population of preys and\n", "$z(t)$ the population of predators. The parameters $\\alpha, \\beta,\n", "\\delta$ and $\\gamma$ depends on the populations to be modeled." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Example 4: Spreading of diseases\n", "\n", "Motivated by the ongoing corona virus pandemic, we consider\n", "a (simple!) model for the spreading of an infectious disease,\n", "which goes under the name [SIR model](https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology#The_SIR_model)." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "The SIR models divides the population into three\n", "population classes, namely\n", "* S(t): number individuals **susceptible** for infection,\n", "* I(t): number **infected** individuals, capable of transmitting the disease,\n", "* R(t): number **removed** individuals who cannot be infected due death or to immunity \n", " after recovery\n", "\n", "The model is of the spreading of a disease is based\n", "on moving individual from $S$ to $I$ and then to $R$.\n", "A short derivation can be found in [[LingeLangtangen2016]](#LingeLangtangen2016) (Ch. 4.2)." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "The final ODE system is given by\n", "\n", "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "S' = - \\beta S I\n", "\\label{ode:eq:sir_model_s} \\tag{3}\n", "\\end{equation}\n", "$$\n", "\n", "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "I' = \\beta S I - \\gamma I\n", "\\label{ode:eq:sir_model_i} \\tag{4}\n", "\\end{equation}\n", "$$\n", "\n", "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "R' = \\gamma I,\n", "\\label{ode:eq:sir_model_r} \\tag{5}\n", "\\end{equation}\n", "$$\n", "\n", "where $\\beta$ denotes the infection rate, and $\\gamma$ the removal rate." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Higher order ODEs\n", "An initial value ODE of order $m$ is given by\n", "\n", "$$\n", "u^{(m)} = f(t,u,u',\\dotsc,u^{(m-1)}), \\qquad u(t_0)=u_0, \\quad\n", "u'(t_0)=u'_0,\\quad \\dotsc, \\quad u^{(m-1)}(t_0) = u^{(m-1)}_0.\n", "$$\n", "\n", "Here $u^{(1)} =u'$ and $u^{(m+1)}=\\frac{du^{(m)}}{dx}$, for $m>0$." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Example 5:\n", "\n", "
\n", "\n", "[Van der Pol's equation](https://en.wikipedia.org/wiki/Van_der_Pol_oscillator)\n", "is a second order differential equation, given by:\n", "\n", "$$\n", "u^{(2)} = \\mu (1-u^2)u' - u, \\qquad u(0)=u_0, \\quad u'(0)=u'_0.\n", "$$\n", "\n", "where $\\mu>0$ is some constant. As initial values $u_0=2$ and\n", "$u'_0=0$ are common choices.\n", "\n", "Later in the note we will see how such equations can be rewritten as a\n", "system of first order ODEs. Systems of higher order ODEs can be\n", "treated similarly." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Euler's method\n", "Now we turn to our first numerical method,\n", "namely\n", "[Euler's method](https://wiki.math.ntnu.no/tma4100/tema/differentialequations?&#numeriske_losninger),\n", "known from Mathematics 1.\n", "We quickly review two alternative derivations,\n", "namely one based on *numerical differentiation*\n", "and one on *numerical integration*." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Derivation of Euler's method\n", "Euler's method is the simplest example of a so-called\n", "**one step method (OSM)**.\n", "Given the IVP\n", "\n", "$$\n", "y'(t) = f(t,y(t)), \\qquad y(t_0)=y_0,\n", "$$\n", "\n", "and some final time $T$,\n", "we want to compute to an approximation of $y(t)$\n", "on $[t_0, T]$." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "We start from $t_0$ and choose some (usually small) time step size\n", "$\\tau_0$ and set the new time $t_1 = t_0 + \\tau_0$. The goal is to\n", "compute a value $y_1$ serving as approximation of $y(t_1)$.\n", "\n", "To do so, we Taylor expand the exact (but unknown) solution\n", "$y(t_0+\\tau)$ around $x_0$:\n", "\n", "$$\n", "y(t_0+\\tau) = y(t_0) + \\tau y'(t_0) + \\frac{1}{2}\\tau^2 y''(t_0) + \\dotsm.\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Assume the step size $\\tau$ to be small, such that the solution is\n", "dominated by the first two terms.\n", "In that case, these can be used as\n", "the numerical approximation in the next step:\n", "\n", "$$\n", "y(t_0+\\tau) \\approx y(t_0) + \\tau y'(t_0) = y_0 + \\tau f(t_0, y_0)\n", "$$\n", "\n", "giving\n", "\n", "$$\n", "y_1 = y_0 + \\tau_0 f(t_0,y_0).\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Now we can repeat this procedure and choose the next\n", "(possibly different) time\n", "step $\\tau_1$ and compute a numerical approximation $y_2$\n", "for $y(t)$ at $t_2 = t_1 + \\tau_1$ by setting\n", "\n", "$$\n", "y_2 = y_1 + \\tau_1 f(t_1,y_1).\n", "$$\n", "\n", "The idea is to repeat this procedure until we reached the\n", "final time $T$ resulting in the following" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "**Recipe: Euler's method.**\n", "\n", "Given a function $f(t,y)$ and an initial value $(t_0,y_0)$.\n", "* Set $t = t_0$.\n", "\n", "* $\\texttt{while } t < T$:\n", "\n", " * Choose $\\tau_k$\n", "\n", " * $\\displaystyle y_{k+1} := y_{k} + \\tau f(t_k, y_k)$ \n", "\n", " * $t_{k+1}:=t_k+\\tau_k$\n", "\n", " * $t := t_{k+1}$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "So we can think of the Euler method as a method\n", "which approximates the continuous but unknown solution\n", "$y(t): [t_0, T] \\to \\mathbb{R}$\n", "by a discrete function\n", "$y_{\\Delta}:\\{t_0, t_1, \\ldots, t_{N_t}\\}$\n", "such that $y_{\\Delta}(t_k) := y_k \\approx y(t_k)$.\n", "\n", "How to choose $\\tau_i$? The simplest possibility\n", "is to set a maximum number of steps $N_{\\mathrm{max}} = N_t$ and then\n", "to chose a *constant time step* $\\tau = (T-t_0)/N_{\\mathrm{max}}$\n", "resulting in $N_{\\mathrm{max}}+1$ equidistributed points.\n", "Later we will also learn, how to choose the\n", "*time step adaptively*, depending on the\n", "solution's behavior." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "**Numerical solution between the nodes.**\n", "\n", "At first we have only an approximation of $y(t)$\n", "at the $N_t +1 $ nodes $y_{\\Delta}:\\{t_0, t_1, \\ldots, t_{N_t}\\}$.\n", "If we want to evaluate the numerical solution between the\n", "nodes, a natural idea is to extend the discrete solution\n", "linearly between each pair of time nodes $t_{k}, t_{k+1}$.\n", "This is compatible with the way the numerical solution can\n", "be plotted, namely by connected each pair\n", "$(t_k, y_k)$ and $(t_{k+1}, y_{k+1})$ with straight lines." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Also, in order to compute an approximation\n", "at the next point $t_{k+1}$,\n", "Euler's method only needs to know $f$, $\\tau_k$\n", "and the solution $y_k$ at the *current* point $t_k$,\n", "but not at earlier points $t_{k-1}, t_{k-2}, \\ldots$\n", "Thus Euler's method\n", "is an prototype of a so-called **One Step Method (OSM)**.\n", "We will formalize this concept later." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "**Interpretation: Euler's method via forward difference operators.**\n", "\n", "After rearranging terms,\n", "we can also interpret the computation of\n", "an approximation $y_1 \\approx y(t_1)$\n", "as replacing the\n", "derivative $y'(t_0) = f(t_0, y_0)$ with a **forward difference operator**\n", "\n", "$$\n", "f(t_0,y_0) = y'(t_0) \\approx\n", " \\dfrac{y(t_1) - y(t_0)}{\\tau}\n", "$$\n", "\n", "Thus *Euler's method replace the differential quotient\n", "by a difference quotient.*\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "**Alternative derivation via numerical integration.**\n", "First we recall that for a function $f: [a,b] \\to \\mathbb{R}$, we can\n", "approximate its integral $\\int_a^b f(t) {\\,\\mathrm{d}t}$ by\n", "a *very simple* quadrature rule of the form\n", "\n", "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\int_a^b f(t) {\\,\\mathrm{d}t} \\approx (b-a) f(a).\n", "\\label{ode:ex:rectangular_qr} \\tag{6}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "**Question.**\n", "\n", "What is the degree of exactness of the previous quadrature rule?" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Turning to our IVP, we know formally integrate\n", "the ODE $y'(t) = f(t, y(t))$ on the time\n", "interval $I_k = [t_k, t_{k+1}]$\n", "and then applying the quadrature rule ([6](#ode:ex:rectangular_qr))\n", "leading to\n", "\n", "$$\n", "y(t_{k+1}) - y(t_k)\n", "=\n", "\\int_{t_k}^{t_{k+1}} y'(t) {\\,\\mathrm{d}t}\n", "=\n", "\\int_{t_k}^{t_{k+1}} f(t,y(t)) {\\,\\mathrm{d}t}\n", "\\approx\n", "\\underbrace{(t_{k+1}-t_{k})}_{\\tau_k}f(t_k, y(t_k))\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "Sorting terms gives us back Euler's method\n", "\n", "$$\n", "y(t_{k+1}) \\approx\n", "y(t_k) + \\tau_k f(t_k, y(t_k)).\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Implementation of Euler's method\n", "Euler's method can be implemented in only a few lines of code:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "def explicit_euler(y0, t0, T, f, Nmax):\n", " ys = [y0]\n", " ts = [t0]\n", " dt = (T - t0) / Nmax\n", " while ts[-1] < T:\n", " t, y = ts[-1], ys[-1]\n", " ys.append(y + dt * f(t, y))\n", " ts.append(t + dt)\n", " return (np.array(ts), np.array(ys))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Let's test Euler's method with the simple IVP given in \n", "Example [Example 1: Population growth and decay processes](#ode:exa:exponential)." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "t0, T = 0, 1\n", "y0 = 1\n", "lam = 1\n", "Nmax = 4\n", "\n", "# rhs of IVP\n", "f = lambda t, y: lam * y\n", "\n", "# Compute numerical solution using Euler\n", "ts, ys_eul = explicit_euler(y0, t0, T, f, Nmax)\n", "\n", "# Exact solution to compare against\n", "y_ex = lambda t: y0 * np.exp(lam * (t - t0))\n", "ys_ex = y_ex(ts)\n", "\n", "# Plot it\n", "plt.plot(ts, ys_ex)\n", "plt.plot(ts, ys_eul, \"ro-\")\n", "plt.legend([\"$y_{ex}$\", \"y\"])" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Plot the solution for various $N_t$,\n", "say $N_t = 4, 8, 16, 32$ against the exact solution given in \n", "Example [Example 1: Population growth and decay processes](#ode:exa:exponential)." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "\n", "## Exercise 1: Error study for the Euler's method\n", "
\n", "\n", "We observed that the more we decrease the constant step size $\\tau$ (or increase $N_{\\mathrm{max}}$),\n", "the closer the numerical solution gets to the exact solution.\n", "\n", "Now we ask you to quantify this. More precisely,\n", "write some code to compute the error\n", "\n", "$$\n", "\\max_{i \\in \\{0, \\ldots, N_{\\mathrm{max}}\\}} |y(t_i) - y_i|\n", "$$\n", "\n", "for $N_{\\mathrm{max}} = 4, 8, 16, 32, 64, 128$.\n", "How does the error reduces if you double the number of points?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "# Insert your code here.\n", "\n", "\n", "def error_study(y0, t0, T, f, Nmax_list, solver, y_ex):\n", " max_errs = []\n", " for Nmax in Nmax_list:\n", " # Compute list of timestep ts and computed solution ys\n", "\n", " # Compute y_ex on ts\n", "\n", " # Compute error\n", "\n", " print(f\"For Nmax = {Nmax:3}, max ||y(t_i) - y_i||= {max_errs[-1]:.3e}\")\n", " print(\"The computed error reduction rates are\")\n", " max_errs = np.array(max_errs)\n", " print(max_errs[:-1] / max_errs[1:])\n", "\n", "\n", "# Define list for N_max\n", "\n", "# Run error study\n", "error_study(y0, t0, T, f, Nmax_list, explicit_euler, y_ex)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "**Solution.**" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "def error_study(y0, t0, T, f, Nmax_list, solver, y_ex):\n", " max_errs = []\n", " for Nmax in Nmax_list:\n", " ts, ys = solver(y0, t0, T, f, Nmax)\n", " ys_ex = y_ex(ts)\n", " errors = ys - ys_ex\n", " max_errs.append(np.abs(errors).max())\n", " print(f\"For Nmax = {Nmax:3}, max ||y(t_i) - y_i||= {max_errs[-1]:.3e}\")\n", " print(\"The computed error reduction rates are\")\n", " max_errs = np.array(max_errs)\n", " print(max_errs[:-1] / max_errs[1:])\n", "\n", "\n", "Nmax_list = [4, 8, 16, 32, 64, 128]\n", "error_study(y0, t0, T, f, Nmax_list, explicit_euler, y_ex)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "\n", "\n", "\n", "\n", "# Heun's method\n", "Before we start to look at more exciting examples,\n", "we will derive a one step method which is more accurate then Euler's\n", "method.\n", "Note that Euler's method can be interpreted as being based on a quadrature\n", "rule with degree of exactness equal to $0$. \n", "Let's try to use a better quadrature rule!" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Again, we start from the *exact representation*,\n", "but this time we use the trapezoidal rule,\n", "which has degree of exactness equal to $1$,\n", "yielding\n", "\n", "$$\n", "\\begin{align*}\n", "y(t_{k+1}) - y(t_k)\n", "&=\n", "\\int_{t_k}^{t_{k+1}} f(t,y(t)) {\\,\\mathrm{d}t}\n", "\\approx\n", "\\dfrac{\\tau_k}{2}\n", "\\left(\n", "f(t_{k+1}, y(t_{k+1})\n", "+\n", "f(t_{k}, y(t_{k})\n", "\\right)\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "This suggest to consider the scheme\n", "\n", "$$\n", "\\begin{align*}\n", "y_{k+1} - y_k\n", "&=\n", "\\dfrac{\\tau_k}{2}\n", "\\left(\n", "f(t_{k+1}, y_{k+1})\n", "+\n", "f(t_{k}, y_{k})\n", "\\right)\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "But note that starting from $y_k$, we cannot immediately compute $y_{k+1}$\n", "as it appears also in the expression $f(t_{k+1}, y_{k+1})$!\n", "This is an example of an **implicit method**!\n", "\n", "To turn this scheme into an explicit scheme,\n", "the idea is now to approximate $y_{k+1}$ appearing in $f$ with an explicit Euler\n", "step:" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "$$\n", "\\begin{align*}\n", "y_{k+1}\n", "&=\n", "y_k +\n", "\\dfrac{\\tau_k}{2}\n", "\\left(\n", "f\\bigl(t_{k+1}, y_{k}+\\tau_k f(t_k, y_k)\\bigr)\n", "+\n", "f(t_{k}, y_k)\n", "\\right).\n", "\\end{align*}\n", "$$\n", "\n", "Observe that we have now nested evaluations of $f$. This can be best\n", "arranged by computing the nested expression in stages, first\n", "the inner one and then the outer one.\n", "This leads to the following recipe." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "**Recipe: Heun's method.**\n", "\n", "Given a function $f(t,y)$ and an initial value $(t_0,y_0)$.\n", "* Set $t = t_0$.\n", "\n", "* $\\texttt{while } t < T$:\n", "\n", " * Choose $\\tau_k$\n", "\n", " * Compute stage $k_1 := f(t_k, y_k)$\n", "\n", " * Compute stage $k_2 := f(t_k+\\tau_k, y_k+\\tau_k k_1)$\n", "\n", " * $\\displaystyle y_{k+1} := y_{k} + \\tfrac{\\tau_k}{2}(k_1 + k_2)$ \n", "\n", " * $t_{k+1}:=t_k+\\tau_k$\n", "\n", " * $t := t_{k+1}$\n", "\n", "\n", "\n", "The function `heun` is implemented in `ode.py`:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "def heun(y0, t0, T, f, Nmax):\n", " ys = [y0]\n", " ts = [t0]\n", " dt = (T - t0) / Nmax\n", " while ts[-1] < T:\n", " t, y = ts[-1], ys[-1]\n", " k1 = f(t, y)\n", " k2 = f(t + dt, y + dt * k1)\n", " ys.append(y + 0.5 * dt * (k1 + k2))\n", " ts.append(t + dt)\n", " return (np.array(ts), np.array(ys))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "\n", "## Exercise 2: Comparing Heun with Euler\n", "
\n", "\n", "\n", "**a)**\n", "Redo the Example [Example 1: Population growth and decay processes](#ode:exa:exponential) with Heun, and plot\n", "both the exact solution, $y_{eul}$ and $y_{heun}$\n", "for $N_t = 4, 8, 16, 32$." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "# Insert code here." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "**Solution.**" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "t0, T = 0, 1\n", "y0 = 1\n", "lam = 1\n", "Nmax = 4\n", "\n", "# rhs of IVP\n", "f = lambda t, y: lam * y\n", "\n", "# Compute numerical solution using Euler and Heun\n", "ts, ys_eul = explicit_euler(y0, t0, T, f, Nmax)\n", "ts, ys_heun = heun(y0, t0, T, f, Nmax)\n", "\n", "# Exact solution to compare against\n", "y_ex = lambda t: y0 * np.exp(lam * (t - t0))\n", "ys_ex = y_ex(ts)\n", "\n", "# Plot it\n", "plt.plot(ts, ys_ex)\n", "plt.plot(ts, ys_eul, \"ro-\")\n", "plt.plot(ts, ys_heun, \"b+-\")\n", "plt.legend([\"$y_{ex}$\", \"$y$ Euler\", \"$y$ Heun\"])" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "\n", "**b)**\n", "Redo [Exercise 1: Error study for the Euler's method](#ode:exe:euler_error_study) with Heun." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "# Insert code here." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "\n", "**Solution.**" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "Nmax_list = [4, 8, 16, 32, 64, 128]\n", "error_study(y0, t0, T, f, Nmax_list, heun, y_ex)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "\n", "\n", "\n", "\n", "# Applying Heun's and Euler's method" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Example 6: The Lotka-Volterra equation revisited\n", "\n", "Solve the Lotka-Volterra equation\n", "\n", "$$\n", "\\begin{align*}\n", "y'(t) &= \\alpha y(t) - \\beta y(t) z(t) \\\\ \n", "z'(t) &= \\delta y(t)z(t) - \\gamma z(t)\n", "\\end{align*}\n", "$$\n", "\n", "In this example, use the parameters and initial values\n", "\n", "$$\n", "\\alpha=2,\\quad \\beta=1, \\quad \\delta=0.5,\\quad \\gamma=1, \\qquad y_{1,0}=2,\n", "\\quad y_{2,0} = 0.5.\n", "$$\n", "\n", "User Euler's method to solve the equation over the interval $[0,20]$,\n", "and use $\\tau=0.02$.\n", "Try also other step sizes, e.g. $\\tau=0.1$ and $\\tau=0.002$. \n", "\n", "**NB!** In this case, the exact solution is not known. What is known is\n", "that the solutions are periodic and positive. Is this the case here?\n", "Check for different values of $\\tau$." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "# Reset plotting parameters\n", "plt.rcParams.update({\"figure.figsize\": (12, 6)})\n", "\n", "# Define rhs\n", "\n", "\n", "def lotka_volterra(t, y):\n", " # Set parameters\n", " alpha, beta, delta, gamma = 2, 1, 0.5, 1\n", " # Define rhs of ODE\n", " dy = np.array(\n", " [alpha * y[0] - beta * y[0] * y[1], delta * y[0] * y[1] - gamma * y[1]]\n", " )\n", " return dy\n", "\n", "\n", "t0, T = 0, 20 # Integration interval\n", "y0 = np.array([2, 0.5]) # Initital values\n", "\n", "# Solve the equation\n", "tau = 0.001\n", "Nmax = int(20 / tau)\n", "print(f\"Nmax = {Nmax:4}\")\n", "ts, ys_eul = explicit_euler(y0, t0, T, lotka_volterra, Nmax)\n", "\n", "# Plot results\n", "plt.plot(ts, ys_eul)\n", "\n", "# Solve the equation\n", "tau = 0.1\n", "Nmax = int(20 / tau)\n", "print(f\"Nmax = {Nmax:4}\")\n", "ts, ys_heun = heun(y0, t0, T, lotka_volterra, Nmax)\n", "\n", "plt.plot(ts, ys_heun)\n", "\n", "plt.xlabel(\"t\")\n", "plt.legend(\n", " [\"$y_0(t)$ - Euler\", \"$y_1(t)$ - Euler\", \"$y_0(t)$ - Heun\", \"$y_1(t)$ - Heun\"],\n", " loc=\"upper right\",\n", ")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "\n", "## Exercise 3: Solving the Lotka-Volterra system using Heun's method\n", "
\n", "\n", "Redo the last example with Heun's method and compare the solutions\n", "generated by Euler's and Heun's method.\n", "\n", "" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Higher order ODEs\n", "How can we numerically solve higher order ODEs\n", "using, e.g., Euler's or Heun's method?\n", "\n", "Given the $m$-th order ODE\n", "\n", "$$\n", "u^{(m)}(t) = f\\big(t, u(t), u'(x), \\dotsc, u^{(m-1)}\\big).\n", "$$\n", "\n", "For a unique solution, we assume that the initial values\n", "\n", "$$\n", "u(t_0), u'(t_0), u''(t_0), \\dotsc, u^{(m-1)}(t_0)\n", "$$\n", "\n", "are known." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Such equations can be written as a system of first order ODEs by the\n", "following trick:\n", "Let\n", "\n", "$$\n", "y_1(x) = u(x), \\quad y_2(x) = u'(x), \\quad\n", "y_3(x) = u^{(2)}(x), \\quad \\dotsc \\quad, y_{m}(x) = u^{(m-1)}(x)\n", "$$\n", "\n", "such that\n", "\n", "$$\n", "\\begin{align*}\n", " y_1' &= y_2, & y_1(a) &= u(a) \\\\ \n", " y_2' &= y_3, & y_2(a) &= u'(a) \\\\ \n", " & \\vdots && \\vdots\\\\ \n", " y_{m-1}' &= y_m, & y_{m-1}(a) &= u^{(m-2)}(a) \\\\ \n", " y_m' &= f(t, y_1, y_2, \\ldots, y_{m-1},y_m), & y_m(a) &= u^{(m-1)}(a)\n", "\\end{align*}\n", "$$\n", "\n", "which is nothing but a system of first order ODEs, and can be solved numerically\n", "exactly as before. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Exercise 4: Numerical solution of Van der Pol's equation\n", "\n", "Recalling Example [Example 5:](#ode:exa:vanderpol), the Van der Pol oscillator\n", "is described by the second order differential equation\n", "\n", "$$\n", "u'' = \\mu (1-u^2)u' - u, \\qquad u(0)=u_0, \\quad u'(0)=u_0'.\n", "$$\n", "\n", "It can be rewritten as a system of first order ODEs:\n", "\n", "$$\n", "\\begin{align*}\n", "y_1' &= y_2, & y_1(0) &= u_0, \\\\ \n", "y_2' &= \\mu(1-y_1^2)y_2 - y_1, & y_2(0) &= u_0'.\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "**a)**\n", "Let $\\mu=2$, $u(0)=2$ and $u'(0)=0$ and solve the equation over the interval\n", "$[0,20]$, using the explicit Euler and $\\tau=0.1$. Play with different\n", "step sizes, and maybe also with different values of $\\mu$.\n", "\n", "**b)**\n", "Repeat the previous numerical experiment with Heun's method.\n", "Try to compare the number of steps you need to perform\n", "with Euler vs Heun to obtain visually the \"same\" solution.\n", "(That is, you measure the difference of the two numerical solutions\n", "in the \"eyeball norm\".)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "# Insert code here." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "**Solution.**" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "# Define the ODE\n", "\n", "\n", "def f(t, y):\n", " mu = 2\n", " dy = np.array([y[1], mu * (1 - y[0] ** 2) * y[1] - y[0]])\n", " return dy\n", "\n", "\n", "# Set initial time, stop time and initial value\n", "t0, T - 0, 20\n", "y0 = np.array([2, 0])\n", "\n", "# Solve the equation using Euler and plot\n", "tau = 0.0001\n", "Nmax = int(20 / tau)\n", "print(f\"Nmax = {Nmax:4}\")\n", "ts, ys_eul = explicit_euler(y0, t0, T, f, Nmax)\n", "plt.plot(ts, ys_eul)\n", "\n", "# Solve the equation using Heun\n", "tau = 0.02\n", "Nmax = int(20 / tau)\n", "print(f\"Nmax = {Nmax:4}\")\n", "ts, ys_heun = heun(y0, t0, T, f, Nmax)\n", "plt.plot(ts, ys_heun)\n", "\n", "plt.xlabel(\"x\")\n", "plt.title(\"Van der Pols ligning\")\n", "plt.legend([\"y1 - Euler\", \"y2 - Euler\", \"y1 - Heun\", \"y2 - Heun\"], loc=\"upper right\");" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Bibliography\n", "1.
**S. Linge and H. P. Langtangen**. \n", " *Programming for Computations - Python*,\n", " 1 edition,\n", " Springer,\n", " 2016." ] } ], "metadata": { "celltoolbar": "Slideshow", "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.9.2" } }, "nbformat": 4, "nbformat_minor": 4 }