{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "\n", "\n", "\n", "# Numerical solution of ordinary differential equations: High order Runge-Kutta methods\n", "\n", " \n", "**André Massing**\n", "\n", "Date: **March 17, 2021**\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": "fragment" } }, "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 import some important Python modules." ] }, { "cell_type": "code", "execution_count": 2, "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": [ "# Runge-Kutta Methods\n", "In the previous lectures we introduced\n", "*Euler's method* and *Heun's method* as\n", "particular instances of the *One Step Methods*,\n", "and we presented the general error theory\n", "for one step method.\n", "\n", "In this Lecture, we introduce a large family\n", "of the one step methods which go under the name\n", "**Runge-Kutta methods (RKM)**. We will see that Euler's method\n", "and Heun's method are instance of RKMs." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Derivation of Runge-Kutta Methods" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "For a given time interval $I_i = [t_i, t_{i+1}]$ we\n", "want to compute $y_{i+1}$ assuming that $y_i$ is given.\n", "Start from the exact expression\n", "\n", "$$\n", "y(t_{i+1}) - y(t_i) = \\int_{t_i}^{t_{i+1}} f(t, y(t)){\\,\\mathrm{d}t},\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "The idea is now to approximate the integral by some quadrature\n", "rule $\\mathrm{Q}[\\cdot](\\{\\xi_j\\}_{j=1}^s,\\{b_j\\}_{j=1}^s)$ defined on $I_i$.\n", "Then we get" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "y(t_{i+1}) - y(t_i) = \\int_{t_i}^{t_{i+1}} f(t, y(t)){\\,\\mathrm{d}t}\n", "\\label{_auto1} \\tag{1}\n", "\\end{equation}\n", "$$\n", "\n", "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "\\approx \\tau \\sum_{j=0}^s b_j f(\\xi_j, y(\\xi_j))\n", "\\label{_auto2} \\tag{2}\n", "\\end{equation}\n", "$$\n", "\n", "Now we can define $\\{c_j\\}_{j=1}^s$ such that $\\xi_j = t_{i} + c_j \\tau$\n", "for $j=1,\\ldots,s$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "\n", "## Exercise 1: A first condition on $b_j$\n", "\n", "\n", "\n", "\n", "\n", "**Question:** What value do you expect for $\\sum_{j=1}^s b_{j}$?\n", "\n", "**Choice A:**\n", " $\\sum_{j=1}^s b_{j} = \\tau$\n", "\n", "**Choice B:**\n", " $\\sum_{j=1}^s b_{j} = 0$\n", "\n", "**Choice C:**\n", " $\\sum_{j=1}^s b_{j} = 1$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "In contrast to pure numerical integration, we don't know the values\n", "of $y(\\xi_j)$. Again, we could use the same idea to approximate\n", "\n", "$$\n", "y(\\xi_j) - y(t_i) = \\int_{t_i}^{t_i+c_j \\tau} f(t, y(t)){\\,\\mathrm{d}t}\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "But then again we get a closure problem if we choose new quadrature points.\n", "The idea is now to *not introduce even more new quadrature points* but to\n", "use same $y(\\xi_j)$ to avoid the closure problem.\n", "Note that this leads to an approximation of the integrals $\\int_{t_i}^{t_i+c_j \\tau}$\n", "with possible nodes *outside* of $[t_i, t_i + c_j \\tau ]$." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "This leads us to\n", "\n", "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "y(\\xi_j) - y(t_i) = \\int_{t_i}^{t_i+c_j \\tau} f(t, y(t)){\\,\\mathrm{d}t}\n", "\\label{_auto3} \\tag{3}\n", "\\end{equation}\n", "$$\n", "\n", "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "\\approx c_j \\tau \\sum_{l=1}^{s}\n", "\\tilde{a}_{jl}\n", "f(\\xi_l, y(\\xi_l))\n", "\\label{_auto4} \\tag{4}\n", "\\end{equation}\n", "$$\n", "\n", "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "=\n", "\\tau \\sum_{l=1}^{s}\n", "{a}_{jl}\n", "f(\\xi_l, y(\\xi_l))\n", "\\label{_auto5} \\tag{5}\n", "\\end{equation}\n", "$$\n", "\n", "where we set $ c_j \\tilde{a}_{jl} = a_{jl}$." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## Exercise 2: A first condition on $a_{jl}$\n", "\n", "\n", "\n", "\n", "\n", "**Question:** What value do you expect for $\\sum_{l=1}^s a_{jl}$?\n", "\n", "**Choice A:**\n", " $\\sum_{l=1}^s a_{jl} = \\tfrac{1}{c_j}$\n", "\n", "**Choice B:**\n", " $\\sum_{l=1}^s a_{jl} = c_j $\n", "\n", "**Choice C:**\n", " $\\sum_{l=1}^s a_{jl} = 1 $\n", "\n", "**Choice D:**\n", " $\\sum_{l=1}^s a_{jl} = \\tau $" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "\n", "\n", "## Definition 1: Runge-Kutta methods\n", "\n", "
\n", "\n", "Given $b_j$, $c_j$, and $a_{jl}$ for $j,l = 1,\\ldots s$, the Runge-Kutta method is\n", "defined by the recipe\n", "\n", "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "Y_{j}\n", "= y_i + \\tau \\sum_{l=1}^{s} {a}_{jl}\n", "f(t_i + c_l \\tau, Y_l) \\quad \\text{for } j = 1,\\ldots s,\n", "\\label{eq:rk-stages} \\tag{6}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "\\label{eq:rk-final} \\tag{7}\n", "y_{i+1} = y_i + \\tau \\sum_{j=1}^s b_j f(t_i + c_j \\tau, Y_j)\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Runge-Kutta schemes are often specified in the form of a **Butcher table**:\n", "\n", "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\label{ode:eq:butchertable} \\tag{8}\n", "\\begin{array}{c|ccc}\n", "c_1 & a_{11} & \\cdots & a_{1s}\n", "\\\\ \n", "\\vdots & \\vdots & & \\vdots\n", "\\\\ \n", "c_s & a_{s1} & \\cdots & a_{ss}\n", "\\\\ \n", "\\hline\n", "& b_1 & \\cdots & b_s\n", "\\end{array}\n", "\\end{equation}\n", "$$\n", "\n", "If $a_{ij} = 0$ for $j \\geqslant i$ the Runge-Kutta method is called **explicit**.\n", "(Why?)\n", "\n", "Note that in the final step, all the function evaluation we need\n", "to perform have already been performed when computing the **stages** $Y_j$." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Therefore one often rewrite the scheme by introducing **stage derivatives**\n", "\n", "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "k_l\n", "= f(t_i + c_l \\tau, Y_l)\n", "\\label{_auto6} \\tag{9}\n", "\\end{equation}\n", "$$\n", "\n", "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", " = f(t_i + c_l \\tau, y_i + \\tau \\sum_{j=1}^{s} {a}_{lj}\n", "k_j) \\quad\n", "j = 1,\\ldots s,\n", "\\label{_auto7} \\tag{10}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "so the resulting scheme will be (swapping index $l$ and $j$)\n", "\n", "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "k_{j} =\n", "f(t_i + c_j \\tau, y_i + \\tau \\sum_{l=1}^{s} {a}_{jl} k_l)\n", "\\quad\n", "j = 1,\\ldots s,\n", "\\label{eq:rk-stage-derivatives} \\tag{11}\n", "\\end{equation}\n", "$$\n", "\n", "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "\\label{eq:rk-final-stage-derivatives} \\tag{12}\n", "y_{i+1} = y_{i} + \\tau \\sum_{j=1}^s b_j k_j\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "\n", "## Exercise 3: Butcher table for the explicit Euler\n", "\n", "Write down the Butcher table for the explicit Euler." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "**Solution.**\n", "Define\n", "$k_1 = f(t_i, y_i) = f(t_i + 0\\cdot \\tau, y_i + \\tau \\cdot 0 \\cdot k_1)$.\n", "Then the explicit Euler step\n", "$y_{i+1} = y_i + \\tau k_1 =\n", "y_i + \\tau \\cdot 1 \\cdot k_1$, \n", "and thus the Butcher table is given by\n", "\n", "$$\n", "\\begin{array}{c|c}\n", "0 & 0\n", "\\\\ \n", "\\hline\n", "& 1\n", "\\end{array}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## Exercise 4: The improved explicit Euler method\n", "\n", "We formally derive the **explicit midpoint rule** or\n", "**improved explicit Euler method**.\n", "Applying the midpoint rule to our integral representatio yields" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "y(t_{i+1}) - y(t_i)\n", "= \\int_{t_i}^{t_{i+1}} f(t, y(t)){\\,\\mathrm{d}t}\n", "\\label{_auto8} \\tag{13}\n", "\\end{equation}\n", "$$\n", "\n", "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "\\approx \\tau f(t_i + \\tfrac{1}{2}\\tau, y(t_i + \\tfrac{1}{2}\\tau))\n", "\\label{_auto9} \\tag{14}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "Since we cannot determine the value $y(t_i + \\tfrac{1}{2}\\tau)$ from this system,\n", "we approximate\n", "it using a half Euler step\n", "\n", "$$\n", "y(t_i + \\tfrac{1}{2}\\tau) \\approx\n", "y(t_i) + \\tfrac{1}{2}\\tau f(t_i, y(t_i))\n", "$$\n", "\n", "leading to the scheme\n", "\n", "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "y_{i+1/2} := y_i + \\tfrac{1}{2}\\tau f(t_i, y_i)\n", "\\label{ode:eq:improved_euler_step_1} \\tag{15}\n", "\\end{equation}\n", "$$\n", "\n", "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "y_{i+1} := y_i + \\tau f(t_i + \\tfrac{1}{2}\\tau, y_{i+1/2})\n", "\\label{ode:eq:improved_euler_step_2} \\tag{16}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "**a)**\n", "Is this a one-step function? Can you define the increment function $\\Phi$?" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "\n", "**Solution.**\n", "Yes it is, and it's increment function is given by\n", "\n", "$$\n", "\\Phi(t_i, y_i, y_{i+1}, \\tau) = f(t_i + \\tfrac{1}{2}\\tau, y_i + \\tfrac{1}{2}\\tau f(t_i, y_i))\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "\n", "**b)**\n", "Can you rewrite this as a Runge-Kutta method?\n", "If so, determine the Butcher table of it." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "**Solution.**\n", "Define $k_1$ and $k_2$ as follows," ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "y_{i+1/2} := y_i + \\tfrac{1}{2}\\tau \\underbrace{f(t_i, y_i)}_{=: k_1}\n", "\\label{ode:eq:improved_euler_k1} \\tag{17}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "y_{i+1}\n", ":= y_i + \\tau f(t_i + \\tfrac{1}{2}\\tau, y_{i+1/2})\n", "= y_i + \\tau \\underbrace{f(t_i + \\tfrac{1}{2}\\tau, y_{i} + \\tau \\tfrac{1}{2} k_1)}_{:= k_2}.\n", "\\label{_auto10} \\tag{18}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "\n", "
\n", "\n", "\n", "Then\n", "\\begin{equation}\n", "y_{i+1} = y_i + \\tau k_2\n", "\\label{ode:eq:improved_euler_k2} \\tag{19}\n", "\\end{equation}" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Thus the Butcher table is given by\n", "\n", "$$\n", "\\begin{array}{c|c c}\n", "0 & 0 & 0\n", "\\\\ \n", "\\tfrac{1}{2} & \\tfrac{1}{2} & 0\n", "\\\\ \n", "\\hline\n", "& 0 & 1\n", "\\end{array}\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "\n", "\n", "\n", "\n", "## Implementation of explicit Runge-Kutta methods" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Below you will find the implementation a general\n", "solver class\n", "`Explicit_Runge_Kutta`\n", "which at its initialization takes\n", "in a Butcher table and has `__call__` function\n", "\n", "```Python\n", " def __call__(self, y0, f, t0, T, n):\n", "```" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "and can be used like this\n", "\n", "```Python\n", " # Define Butcher table\n", " a = np.array([[0, 0, 0],\n", " [1.0/3.0, 0, 0],\n", " [0, 2.0/3.0, 0]])\n", " \n", " b = np.array([1.0/4.0, 0, 3.0/4.0])\n", " \n", " c = np.array([0, 1.0/3.0, 2.0/3.0])\n", "```" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "```Python\n", " # Create solver\n", " rk3 = Explicit_Runge_Kutta(a, b, c)\n", " \n", " # Solve problem (applies __call__ function)\n", " ts, ys = rk3(y0, t0, T, f, Nmax)\n", "```" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "The complete implementation is given here:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "class Explicit_Runge_Kutta:\n", " def __init__(self, a, b, c):\n", " self.a = a\n", " self.b = b\n", " self.c = c\n", "\n", " def __call__(self, y0, t0, T, f, Nmax):\n", " # Extract Butcher table\n", " a, b, c = self.a, self.b, self.c\n", "\n", " # Stages\n", " s = len(b)\n", " ks = [np.zeros_like(y0, dtype=np.double) for s in range(s)]\n", "\n", " # Start time-stepping\n", " ys = [y0]\n", " ts = [t0]\n", " dt = (T - t0) / Nmax\n", "\n", " while ts[-1] < T:\n", " t, y = ts[-1], ys[-1]\n", "\n", " # Compute stages derivatives k_j\n", " for j in range(s):\n", " t_j = t + c[j] * dt\n", " dY_j = np.zeros_like(y, dtype=np.double)\n", " for l in range(j):\n", " dY_j += dt * a[j, l] * ks[l]\n", "\n", " ks[j] = f(t_j, y + dY_j)\n", "\n", " # Compute next time-step\n", " dy = np.zeros_like(y, dtype=np.double)\n", " for j in range(s):\n", " dy += dt * b[j] * ks[j]\n", "\n", " ys.append(y + dy)\n", " ts.append(t + dt)\n", "\n", " return (np.array(ts), np.array(ys))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Example 1: Implementation and testing of the improved Euler method\n", "\n", "\n", "We implement the **improved explicit Euler** from above and\n", "plot the analytical and the numerical solution. \n", "Finally, we determine the convergence order." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "def compute_eoc(y0, t0, T, f, Nmax_list, solver, y_ex):\n", " errs = []\n", " for Nmax in Nmax_list:\n", " ts, ys = solver(y0, t0, T, f, Nmax)\n", " ys_ex = y_ex(ts)\n", " errs.append(np.abs(ys - ys_ex).max())\n", " print(f\"For Nmax = {Nmax:3}, max ||y(t_i) - y_i||= {errs[-1]:.3e}\")\n", "\n", " errs = np.array(errs)\n", " Nmax_list = np.array(Nmax_list)\n", " dts = (T - t0) / Nmax_list\n", "\n", " eocs = np.log(errs[1:] / errs[:-1]) / np.log(dts[1:] / dts[:-1])\n", " return errs, eocs" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "For Nmax = 4, max ||y(t_i) - y_i||= 2.343e-02\n", "For Nmax = 8, max ||y(t_i) - y_i||= 6.441e-03\n", "For Nmax = 16, max ||y(t_i) - y_i||= 1.688e-03\n", "For Nmax = 32, max ||y(t_i) - y_i||= 4.322e-04\n", "For Nmax = 64, max ||y(t_i) - y_i||= 1.093e-04\n", "For Nmax = 128, max ||y(t_i) - y_i||= 2.749e-05\n", "[2.34261385e-02 6.44058991e-03 1.68830598e-03 4.32154479e-04\n", " 1.09316895e-04 2.74901378e-05]\n", "[1.86285442 1.93161644 1.96595738 1.98303072 1.99153035]\n" ] } ], "source": [ "# Define Butcher table for improved Euler\n", "a = np.array([[0, 0], [0.5, 0]])\n", "b = np.array([0, 1])\n", "c = np.array([0, 0.5])\n", "\n", "# Define rk2\n", "rk2 = Explicit_Runge_Kutta(a, b, c)\n", "\n", "t0, T = 0, 1\n", "y0 = 1\n", "lam = 1\n", "\n", "# rhs of IVP\n", "f = lambda t, y: lam * y\n", "\n", "# Exact solution to compare against\n", "y_ex = lambda t: y0 * np.exp(lam * (t - t0))\n", "\n", "# EOC test\n", "Nmax_list = [4, 8, 16, 32, 64, 128]\n", "\n", "errs, eocs = compute_eoc(y0, t0, T, f, Nmax_list, rk2, y_ex)\n", "print(errs)\n", "print(eocs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercise 5\n", "\n", "If you have time, determine the experimental order of convergence **EOC**\n", "for the following methods:\n", "\n", "**Heun's method with 3 stages**\n", "$$\n", "\\begin{array}{c|c c c}\n", "0 & 0 & 0 & 0\n", "\\\\ \n", "\\tfrac{1}{3} & \\tfrac{1}{3} & 0 & 0\n", "\\\\ \n", "\\tfrac{2}{3} & 0 & \\tfrac{2}{3} & 0 \n", "\\\\\n", "\\hline\n", "& \\tfrac{1}{4} & 0 & \\tfrac{3}{4} \n", "\\end{array}\n", "$$" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Insert code here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**The classical Runge-Kutta method with 4 stages**\n", "$$\n", "\\begin{array}{c|c c c c}\n", "0 & 0 & 0 & 0 & 0\n", "\\\\ \n", "\\tfrac{1}{2} & \\tfrac{1}{2} & 0 & 0 & 0 \n", "\\\\ \n", "\\tfrac{1}{2} & 0 & \\tfrac{1}{2} & 0 & 0\n", "\\\\\n", "1 & 0 & 0 & 1 & 0\n", "\\\\\n", "\\hline\n", "& \\tfrac{1}{6} & \\tfrac{1}{3}& \\tfrac{1}{3} & \\tfrac{1}{6}\n", "\\end{array}\n", "$$" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Insert code here" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "\n", "\n", "\n", "\n", "## Order conditions for Runge-Kutta Methods\n", "The convergence theorem for one-step methods\n", "gave us some necessary conditions to guarantee\n", "that a method is convergent order of $p$:\n", "* \"consistency order $p$\" + \"Increment function satisfies a Lipschitz condition\"\n", " $\\Rightarrow$ \"convergence order $p$.\n", " \n", "To put it in slightly different words:\n", " \n", "* \"local truncation error behaves like $\\mathcal{O}(\\tau^{p+1})$\" + \"Increment function \n", " satisfies a Lipschitz condition\" \n", " $\\Rightarrow$ \"global truncation error behaves like $\\mathcal{O}(\\tau^{p})$\"\n", "\n", "It turns out that for $f$ is at least $C^1$ with respect to all\n", "its arguments then the increment function $\\Phi$\n", "associated with any Runge-Kutta methods satisfies\n", "a Lipschitz condition." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Theorem 1: Order conditions for Runge-Kutta methods\n", "\n", "
\n", "A Runge - Kutta method has consistency order $p$ if and only if all the\n", "conditions up to and including $p$ in the table below are satisfied.\n", "\n", "$$\n", "\\begin{array}{c|c|c} \n", " p & \\text{conditions} \\\\ \\hline \n", " 1 & \\sum b_i = 1 \\\\ \\hline \n", " 2 & \\sum b_i c_i = 1/2 \\\\ \\hline \n", " 3 & \\sum b_i c_i^2 = 1/3\\\\ \n", " & \\sum b_i a_{ij} c_j = 1/6 \n", " \\\\ \\hline \n", " 4 & \\sum b_ic_i^3=1/4 \\\\ \n", " & \\sum b_i c_i a_{ij}c_j=1/8 \\\\ \n", " & \\sum b_i a_{ij}c_j^2=1/12 \\\\ \n", " & \\sum b_i a_{ij} a_{jk} c_k = 1/24 \\\\ \\hline \n", "\\end{array}\n", "$$\n", "\n", "where sums are taken over all the indices from 1 to $s$.\n", "\n", " \n", "**Proof.**\n", "Without proof." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Theorem 2: Convergence theorem for Runge-Kutta methods\n", "\n", "
\n", "Given the IVP ${\\boldsymbol y}' = {\\boldsymbol f}(t, {\\boldsymbol y}), {\\boldsymbol y}(0) = {\\boldsymbol y}_0$.\n", "Assume $f \\in C^1$ and that a given Runge-Kutta method satisfies\n", "the order conditions\n", "from Theorem [Theorem 1: Order conditions for Runge-Kutta methods](#thm:rk-order-conditions)\n", "up to order $p$.\n", "Then the Runge-Kutta method is convergent of order $p$." ] } ], "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 }