{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0b88f670",
   "metadata": {},
   "source": [
    "# Solving differential equations\n",
    "\n",
    "- by Börge Göbel"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ed1b45a",
   "metadata": {},
   "source": [
    "## 1. Euler method"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86b7b82b",
   "metadata": {},
   "source": [
    "## 1.1 First order differential equation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e05912fb",
   "metadata": {},
   "source": [
    "We try to solve the following type of differential equation\n",
    "\n",
    "\\\\( \\frac{\\mathrm{d}y}{\\mathrm{d}t} = f(t,y)\\\\)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0127bce0",
   "metadata": {},
   "source": [
    "Since \\\\( \\frac{\\mathrm{d}y}{\\mathrm{d}t} = \\frac{y(t+h)-y(t)}{h}\\\\), we know that \\\\( y(t+h) = \\frac{\\mathrm{d}y}{\\mathrm{d}t}h + y(t)\\\\).\n",
    "\n",
    "Therefore, we can repetitively iterate the propagation: \n",
    "\n",
    "From the value \\\\( y_n \\\\) at step \\\\( n \\\\), corresponding to the time \\\\( t \\\\), we can calculate the value \\\\( y_{n+1} \\\\) at step \\\\( (n+1) \\\\), corresponding to the time \\\\( (t+h) \\\\):\n",
    "\n",
    "\\\\( y_{n+1} = y_n + \\frac{\\mathrm{d}y}{\\mathrm{d}t}h \\\\) which is \n",
    "\n",
    "\\\\( y_{n+1} = y_n + f(t,y_n) h \\\\)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eaf63e10",
   "metadata": {},
   "source": [
    "### Example 1) Radioactive decay"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82db4e00",
   "metadata": {},
   "source": [
    "\\\\( \\dot{y} = -y\\\\) or\n",
    "\n",
    "\\\\( \\frac{\\mathrm{d}y}{\\mathrm{d}t} = f(t,y) = -y\\\\)\n",
    "\n",
    "Analytical solution: \\\\( y(t)=y_0 \\exp(-t)\\\\)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "0ae494cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import integrate\n",
    "import matplotlib.pyplot as plt "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ccd3a0f0",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "5af1c96a",
   "metadata": {},
   "source": [
    "### Define a function \"eulerODE\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4753714b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "b51cca27",
   "metadata": {},
   "source": [
    "### Example 2) Time-amplified decay"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89295289",
   "metadata": {},
   "source": [
    "\\\\( \\dot{y} = -ayt\\\\) or\n",
    "\n",
    "\\\\( \\frac{\\mathrm{d}y}{\\mathrm{d}t} = f(t,y) = -ayt\\\\)\n",
    "\n",
    "Analytical solution: \\\\( y(t)=y_0 \\exp(-t^2a/2)\\\\)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3448bd17",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "83b20246",
   "metadata": {},
   "source": [
    "### 1.2 Higher-order differential equations"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44f1a8b2",
   "metadata": {},
   "source": [
    "Example: Second-order differential equation: \\\\( y''(t) = f\\left(t,y(t),y'(t)\\right)\\\\)\n",
    "\n",
    "Introduce: \\\\( z_0(t) = y(t)\\\\) and \\\\( z_1(t) = y'(t)\\\\)\n",
    "\n",
    "\\\\( \\begin{pmatrix}z_0'(t)\\\\z_1'(t)\\end{pmatrix}=\\begin{pmatrix}z_1(t)\\\\f\\left(t,z_0(t),z_1(t)\\right)\\end{pmatrix}\\\\)\n",
    "\n",
    "Therefore, we can describe the second-order differential equation by a set of two first-order differential equations. We can solve both with our Euler method\n",
    "\n",
    "\\\\( z_0^{(n+1)} = z_0^{(n)} + z_1^{(n)} h \\\\)\n",
    "\n",
    "\\\\( z_1^{(n+1)} = z_1^{(n)} + f\\left(t,z_0^{(n)},z_1^{(n)}\\right) h \\\\)\n",
    "\n",
    "Or, going back to our initial nomenclature:\n",
    "\n",
    "\\\\( y_{n+1} = y_{n} + y'_{n} h \\\\)\n",
    "\n",
    "\\\\( y'_{n+1} = y'_{n} + f\\left(t,y_{n},y'_{n}\\right) h \\\\)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c0de8cd",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "e533b5c2",
   "metadata": {},
   "source": [
    "### Example 3) Free fall"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4378cc6a",
   "metadata": {},
   "source": [
    "\\\\( \\ddot{y} = -g\\\\) or\n",
    "\n",
    "\\\\( \\frac{\\mathrm{d}^2y}{\\mathrm{d}t^2} = f(t,y,\\dot{y}) = -g\\\\)\n",
    "\n",
    "Analytical solution: \\\\( y(t)=-\\frac{g}{2}t^2+v_0t+y_0\\\\)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "968e9392",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "13ec1509",
   "metadata": {},
   "source": [
    "### Example 4) Harmonic oscillator"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8174ec35",
   "metadata": {},
   "source": [
    "\\\\( \\theta''(t) + b\\theta'(t) + c\\sin(\\theta(t)) = 0 \\\\)\n",
    "\n",
    "Here, \\\\( b \\\\) is the damping parameter and \\\\( c \\\\) is determined by the pendulum length \\\\( c = \\frac{g}{l} \\\\)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "017fc71a",
   "metadata": {},
   "source": [
    "### Small-angle approximation\n",
    "\n",
    "For small angles \\\\( \\theta\\ll 1 \\\\) and without damping b = 0, we have \n",
    "\n",
    "\\\\( \\theta''(t) = - \\frac{g}{l}\\theta(t) \\\\) with the solution (for \\\\( \\theta'(0) = 0 \\\\))\n",
    "\n",
    "\\\\( \\theta(t) = \\theta_0\\cos\\left(\\sqrt{\\frac{g}{l}}t\\right) \\\\) and a period of \\\\( T = 2\\pi\\sqrt{\\frac{l}{g}} \\\\)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc04ebb9",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "e00f08a4",
   "metadata": {},
   "source": [
    "### Actual equation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a50af7fa",
   "metadata": {},
   "source": [
    "- Small starting angle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee045fd4",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "3757d0b4",
   "metadata": {},
   "source": [
    "- Large starting angle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c4983b8",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "8d85c939",
   "metadata": {},
   "source": [
    "- With damping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61b3a93d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "5b0ef086",
   "metadata": {},
   "source": [
    "- Driven oscillator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "747be128",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "2c204ea4",
   "metadata": {},
   "source": [
    "# 2. Improved methods"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05c4a7f2",
   "metadata": {},
   "source": [
    "The exist two useful solvers:\n",
    "- Old: scipy.integrate.oldeint\n",
    "- New: scipy.integrate.solve_ivp"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19a8d176",
   "metadata": {},
   "source": [
    "### Example 2) Time-amplified decay"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c16e98a1",
   "metadata": {},
   "source": [
    "\\\\( \\dot{y} = -ayt\\\\) or\n",
    "\n",
    "\\\\( \\frac{\\mathrm{d}y}{\\mathrm{d}t} = f(t,y) = -ayt\\\\)\n",
    "\n",
    "Analytical solution: \\\\( y(t)=y_0 \\exp(-t^2a/2)\\\\)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dbd763f0",
   "metadata": {},
   "source": [
    "- Our old results (Euler method)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb9000b8",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "d6815209",
   "metadata": {},
   "source": [
    "- New results (solve_ivp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "798b4b51",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "b6002995",
   "metadata": {},
   "source": [
    "### Example 3) Free fall"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95a05a86",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "4ddae053",
   "metadata": {},
   "source": [
    "### Example 4) Driven pendulum"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2414cfb",
   "metadata": {},
   "source": [
    "- Our old results (Euler method)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c20e1e3e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "2b4e0271",
   "metadata": {},
   "source": [
    "- New results (solve_ivp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f69cb338",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "ea653680",
   "metadata": {},
   "source": [
    "### Compare more methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "b9c20a0c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# [https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html#scipy.integrate.solve_ivp]\n",
    "# methods:\n",
    "# RK45\n",
    "# RK23\n",
    "# DOP853\n",
    "# Radau\n",
    "# BDF\n",
    "# LSODA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fbb75ec2",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "d8a5d37a",
   "metadata": {},
   "source": [
    "## 3. Theory of the Runge-Kutta methods"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc614b20",
   "metadata": {},
   "source": [
    "There exist several different Runge-Kutta methods"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd89d4bd",
   "metadata": {},
   "source": [
    "Derivation is difficult: https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods#Derivation_of_the_Runge%E2%80%93Kutta_fourth-order_method"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "870730fb",
   "metadata": {},
   "source": [
    "### 3.1 Implementation of RK4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "295c4b3f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# [https://en.wikipedia.org/wiki/List_of_Runge%E2%80%93Kutta_methods#Classic_fourth-order_method]\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae570923",
   "metadata": {},
   "source": [
    "### 3.2 Implementation of RK45"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "1f76f261",
   "metadata": {},
   "outputs": [],
   "source": [
    "# [https://en.wikipedia.org/wiki/List_of_Runge%E2%80%93Kutta_methods#Fehlberg]\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "377a45b8",
   "metadata": {},
   "source": [
    "### 3.3 Comparison with Euler method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a57edfc9",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}