{ "cells": [ { "cell_type": "markdown", "metadata": { "toc": "true" }, "source": [ "# Table of Contents\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Runge-Kutta methods for ODE integration in Python\n", "\n", "- I want to implement and illustrate the [Runge-Kutta method](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods) (actually, different variants), in the [Python programming language](https://www.python.org/).\n", "\n", "- The Runge-Kutta methods are a family of numerical iterative algorithms to approximate solutions of [Ordinary Differential Equations](https://en.wikipedia.org/wiki/Ordinary_differential_equation). I will simply implement them, for the mathematical descriptions, I let the interested reader refer to the Wikipedia page, or [any](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods#References) [good](https://www.directtextbook.com/isbn/9780521007948) [book](https://www.decitre.fr/livres/analyse-numerique-et-equations-differentielles-9782868838919.html) or [course](https://courses.maths.ox.ac.uk/node/4294) on numerical integration of ODE.\n", "- I will start with the order 1 method, then the order 2 and the most famous order 4.\n", "- They will be compared on different ODE." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Preliminary" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2017-11-23T19:18:23+01:00\n", "\n", "CPython 3.6.3\n", "IPython 6.2.1\n", "\n", "compiler : GCC 7.2.0\n", "system : Linux\n", "release : 4.13.0-16-generic\n", "machine : x86_64\n", "processor : x86_64\n", "CPU cores : 4\n", "interpreter: 64bit\n" ] } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "%load_ext watermark\n", "%watermark" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "from scipy.integrate import odeint # for comparison" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "I will use as a first example the one included in [the scipy documentation for this `odeint` function](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html).\n", "\n", "$$\\theta''(t) + b \\theta'(t) + c \\sin(\\theta(t)) = 0.$$\n", "\n", "If $\\omega(t) = \\theta'(t)$, this gives\n", "$$ \\begin{cases}\n", "\\theta'(t) = \\omega(t) \\\\\n", "\\omega'(t) = -b \\omega(t) - c \\sin(\\theta(t))\n", "\\end{cases} $$\n", "\n", "Vectorially, if $y(t) = [\\theta(t), \\omega(t)]$, then the equation is $y' = f(t, y)$ where $f(t, y) = [y_2(t), -b y_2(t) - c \\sin(y_1(t))]$." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def pend(y, t, b, c):\n", " return np.array([y[1], -b*y[1] - c*np.sin(y[0])])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We assume the values of $b$ and $c$ to be known, and the starting point to be also fixed:" ] }, { "cell_type": "code", "execution_count": 44, "metadata": {}, "outputs": [], "source": [ "b = 0.25\n", "c = 5.0\n", "y0 = np.array([np.pi - 0.1, 0.0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `odeint` function will be used to solve this ODE on the interval $t \\in [0, 10]$, with $101$ points." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "t = np.linspace(0, 10, 101)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is used like this, and our implementations will follow this signature." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "sol = odeint(pend, y0, t, args=(b, c))" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[