Eerste versie nieuwe python opgave KMa 2020

Mathematische_slinger.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Huiswerk 4: Gedempte mathematische slinger\n",
    "#### Klasieke Mechanica a, 28 februari 2020\n",
    "\n",
    "Werk door het onderstaande notebook en beantwoord vraag **(a)** tot en met **(g)**. Voeg waar nodig relevante code en plots toe aan je antwoorden en lever deze **op papier** in **voor het begin van** het hoorcollege van **6 maart 2020**.\n",
    "\n",
    "We bekijken het gedrag van een gedempte mathematische slinger, wat gegeven wordt door de volgende non-lineaire bewegingsvergelijking:\n",
    "\n",
    "$$m\\frac{ds^2}{dt^2} = -c\\frac{ds}{dt} - k \\sin(s)$$\n",
    "\n",
    "Deze bewegingsvergelijking is niet alleen non-lineair, maar het is ook een _tweede-orde_ differentiaalvergelijking, omdat er een tweede orde afgeleide in voorkomt."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**a)** Geef de gelineariseerde differentiaalvergelijking, bereken de gedempte oscillatiefrequentie van het gelineariseerde systeem en bereken de bijbehorende Q factor."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerieke methoden voor gewone differentiaalvergelijkingen\n",
    "We gaan het verschil tussen het gelineairiseerde systeem en de originele differentiaalvergelijking bestuderen. Omdat het analytisch oplossen van non-lineaire differentiaalvergelijkingen bijna altijd onmogelijk is, zullen we dit numeriek gaan doen. Een algemene eerste orde differentiaalvergelijking kan als volgt beschreven worden:\n",
    "$$\\frac{dx}{dt} = f(x,t)$$\n",
    "Voor een algemene functie $f(x,t)$. Gegeven een beginvoorwaarde $x(0) = x_0$, kunnen we numeriek $x(t)$ gaan benaderen met de volgende formule:\n",
    "$$x(t+\\delta t) \\approx x(t) + \\delta t \\cdot f(x(t), t)$$\n",
    "Hiermee kunnen we voor een gekozen _stapgrootte_ $\\delta t$ een reeks posities $x(t_n)$ uitrekenen door iteratie van de volgende recurrente betrekking:\n",
    "$$x(t_{n+1}) \\approx x(t_n) + \\delta t \\cdot f(x(t_n), t_n)$$\n",
    "\n",
    "Het toepassen van deze relatie om opeenvolgende waardes uit te rekenen heet Voorwaarts Euler.\n",
    "Laten we dit ter illustratie doen voor een functie $f(x) = -0.5 \\cdot x$ en beginvoorwaarde $x_0=1$. Dit zou een dalende exponent moeten geven:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def f(x):\n",
    "    return -0.5*x\n",
    "\n",
    "def FE_step(func, x, dt):\n",
    "    \"\"\"Perform a single forward Euler step\"\"\"\n",
    "    return x + dt * func(x)\n",
    "\n",
    "def RK4_step(func, x, dt):\n",
    "    \"\"\"perform a single Runge Kutta 4 step\n",
    "    for time independent func\"\"\"\n",
    "    k1 = dt * func(x)\n",
    "    k2 = dt * func(x + k1/2)\n",
    "    k3 = dt * func(x + k2/2)\n",
    "    k4 = dt * func(x + k3)\n",
    "    return x + (k1+2*k2+2*k3+k4)/6\n",
    "\n",
    "x = np.empty(100)\n",
    "x[0] = 1\n",
    "dt = 0.1\n",
    "\n",
    "\n",
    "for n in range(100 - 1):\n",
    "    x[n+1] = FE_step(f, x[n], dt)\n",
    "    # equivalent to:\n",
    "    # x[n+1] = x[n] + dt * f(x[n])\n",
    "\n",
    "t = np.arange(0, 100*dt, dt)\n",
    "    \n",
    "plt.plot(t, x)\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('x');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hogere orde differentiaal vergelijkingen\n",
    "Voor tweede orde differentiaal vergelijkingen (dus met een tweede afgeleide erin) werkt dit niet direct, omdat we de tweede afgeleide niet zo kunnen omschrijven. We kunnen echter wel een $n$-de orde differentiaal vergelijking omschrijven in een stelsel van $n$ verschillende eerste orde differentiaalvergelijkingen door iedere afgeleide als een aparte variabele op te vatten, en dan dezelfde numerieke methode toe te passen. Hiervoor introduceren we $v = \\dot{s}$. Omschrijven geeft:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align}\n",
    "\\frac{ds}{dt} &= \\dot{s} = v\\\\\n",
    "S &= \\binom{s}{\\dot{s}} = \\binom{s}{v}\\\\\n",
    "\\dot{S} &= \\binom{\\dot{s}}{\\ddot{s}} = \\binom{\\dot{s}}{\\dot{v}} = \\binom{v}{-\\frac{c}{m}v - \\frac{k}{m} \\sin(s)}\\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "De laatste regel is nu in een vorm geschikt om Voorwaarts Euler op toe te passen om een numerieke oplossing te verkrijgen. Nu moeten we dit nog omzetten in een Python functie. Hieronder is de functie voor de mathematische slinger gegeven.\n",
    "\n",
    "**b)** Vul zelf de functie voor het gelineariseerde systeem in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def math_slinger(S, c=0.05, k=3, m=1):\n",
    "    s,v = S\n",
    "    dv = - c/m * v - k/m * np.sin(s)\n",
    "    ds = v\n",
    "    return np.array((ds, dv))\n",
    "\n",
    "#opgave: Schrijf de functie voor het gelineariseerde systeem hieronder\n",
    "def math_slinger_lin(S, c=0.05, k=3, m=1):\n",
    "    s,v = S\n",
    "    dv = - c/m * v - k/m * s\n",
    "    ds = v\n",
    "    return np.array((ds, dv))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nu kunnen we voor beide systemen integreren en een plot maken om te visualiseren wat er gebeurd:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dt = 0.02\n",
    "t = np.arange(0, 100, dt)\n",
    "S = np.empty((len(t), 2))\n",
    "S[0][0] = 2\n",
    "S[0][1] = -2\n",
    "S_lin = np.copy(S)\n",
    "for n in range(len(t)-1):\n",
    "    S_lin[n+1] = FE_step(math_slinger_lin, S_lin[n], dt)\n",
    "    S[n+1] = FE_step(math_slinger, S[n], dt)\n",
    "    \n",
    "\n",
    "plt.figure(figsize=[15,6])\n",
    "plt.plot(t, S[:,0], label='mathematical oscillator')\n",
    "plt.plot(t, S_lin[:,0], label='harmonic oscillator')\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('s')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**c)** Wat is het verschil tussen het gelineariseerde systeem en de 'echte' mathematische slinger? Waarom kan dit verschil niet fysisch correct zijn?\n",
    "\n",
    "Wat er gebeurt is dat de numerieke methode instabiel is: de error per stap telt op en veroorzaakt incorrect/onfysisch gedrag. Er zijn twee manieren om dit op te lossen: ofwel we verkleinen de stapgrootte tot het gedrag convergeert, of we gebruiken een geavanceerdere numerieke methode om per stap een betere schatting te maken. We gaan allebei proberen.\n",
    "\n",
    "**d)** Bepaal de maximale waarde voor `dt` waarvoor de eerder gebruikte Voorwaarts Euler methode stabiel is, dat wil zeggen: de maximale waarde voor `dt` waarvoor de relatieve fout niet groeit als functie van de tijd. Voeg een plot toe om je antwoord te illustreren.\n",
    "\n",
    "**e)** Doe hetzelfde als in **d)**, maar voor de Runge-Kutta 4 methode, door `FE_step()` te vervangen door de functie `RK4_step()` die eerder gedefinieerd is. Welke van de twee methodes heeft je voorkeur en waarom?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=[15,6])\n",
    "for dt in [0.02, 0.01, 0.005, 0.001, 0.0005, 0.0001]:\n",
    "    t = np.arange(0, 100, dt)\n",
    "    S = np.empty((len(t), 2))\n",
    "    S[0][0] = 2\n",
    "    S[0][1] = -2\n",
    "    for n in range(len(t)-1):\n",
    "        S[n+1] = FE_step(math_slinger_lin, S[n], dt)\n",
    "        #S[n+1] = S[n] + dt * math_slinger_lin(S[n])\n",
    "    plt.plot(t,S[:,0], label=f'FE {dt:.4f}')\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('s')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=[15,6])\n",
    "for dt in [1, 0.5, 0.1, 0.02, 0.01]:\n",
    "    t = np.arange(0, 100, dt)\n",
    "    S = np.empty((len(t), 2))\n",
    "    S[0][0] = 2\n",
    "    S[0][1] = -2\n",
    "    for n in range(len(t)-1):\n",
    "        S[n+1] = RK4_step(math_slinger_lin, S[n], dt)\n",
    "    plt.plot(t, S[:,0], label=f'RK4 {dt:.4f}')\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('s')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nu we de stabiliteit van onze numerieke methode hebben bepaald aan de hand van het gelineariseerde systeem, is het tijd om het verschil tussen de harmonische oscillator en de mathematische slinger te gaan bestuderen. Kies hieronder je favoriete numerieke methode en een geschikte stepsize `dt` en vergelijk het gedrag van beide systemen voor verschillende beginvoorwaarden."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dt = 0.01 #Vul in\n",
    "\n",
    "t = np.arange(0, 100, dt)\n",
    "S = np.empty((len(t), 2))\n",
    "S[0][0] = 2  # Beginpositie: varieer\n",
    "S[0][1] = -1  # Beginsnelheid: varieer\n",
    "S_lin = np.copy(S)\n",
    "for n in range(len(t)-1):\n",
    "    S_lin[n+1] = RK4_step(math_slinger_lin, S_lin[n], dt)\n",
    "    S[n+1] = RK4_step(math_slinger, S[n], dt)\n",
    "    \n",
    "\n",
    "plt.figure(figsize=[15,6])\n",
    "plt.plot(t, S[:,0], label='mathematical pendulum')\n",
    "plt.plot(t, S_lin[:,0], label='harmonic oscillator')\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('s')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**f)** Beschrijf de geobserveerde verschillen tussen de mathematische slinger en de harmonische oscillator. Voor watvoor beginvoorwaarden zijn de verschillen het grootst? \n",
    "\n",
    "**g)** Bestudeer het verschil in periode en amplitude als functie van initiele uitwijking en snelheid. \n",
    "\n",
    "Hiervoor zijn twee functies aan te raden: [`functools.partial`](https://docs.python.org/3.7/library/functools.html#functools.partial) en [`scipy.signal.argrelmax()`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.argrelmax.html). De eerste laat je de differentiaal functies doorgeven met aangepaste parameters, de tweede kan je gebruiken om de posities en amplitudes van de maxima en minima te vinden, zoals beide hieronder gebruikt.\n",
    "\n",
    "Maak 1 of meer plots die de verschillen tussen de mathematische slinger en de harmonische oscillator als functie van de parameters zo goed mogelijk samenvatten en voeg deze bij je antwoorden."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.signal import argrelmax\n",
    "from functools import partial\n",
    "\n",
    "S = np.empty((len(t), 2))\n",
    "S[0][0] = 1.1\n",
    "S[0][1] = 0.1\n",
    "S_lin = np.copy(S)\n",
    "for n in range(len(t)-1):\n",
    "    S_lin[n+1] = RK4_step(partial(math_slinger_lin, c=0.3), S_lin[n], dt)\n",
    "    S[n+1] = RK4_step(partial(math_slinger, c=0.3), S[n], dt)\n",
    "    \n",
    "extremum_indices = argrelmax(np.abs(S[:,0]))[0]\n",
    "    \n",
    "plt.figure(figsize=[15,6])\n",
    "plt.plot(t, S[:,0], label='mathematical pendulum')\n",
    "plt.plot(t, S_lin[:,0], label='harmonic oscillator')\n",
    "plt.scatter(t[extremum_indices], S[extremum_indices, 0])\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('s')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "from numba import njit\n",
    "nb_math_slinger_lin = njit(math_slinger_lin)\n",
    "@njit()\n",
    "def calc_with_FE(s0, v0, t, dt):\n",
    "    S = np.empty((len(t), 2), dtype=np.float64)\n",
    "    S[0][0] = s0\n",
    "    S[0][1] = v0\n",
    "    for n in range(len(t)-1):\n",
    "        S[n+1] = S[n] + dt * nb_math_slinger_lin(S[n])\n",
    "    return S"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=[15,6])\n",
    "for dt in [0.02, 0.01, 0.005, 0.001, 0.0005, 0.0001]:\n",
    "    t = np.arange(0, 100, dt)\n",
    "    S = calc_with_FE(2, -2, t, dt)\n",
    "    plt.plot(t,S[:,0], label=f'FE {dt:.4f}')\n",
    "plt.legend()"
   ]
  }
 ],
 "metadata": {
  "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}

requirements.txt

numpy
matplotlib
scipy
numba