{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TD1: Temporal Discretization\n",
    "\n",
    "With temporal discretization methods we can numerically solve Cauchy problems such as:\n",
    "$$\n",
    "    \\begin{cases}\n",
    "        u^\\prime = f(u)\\\\\n",
    "        u(0) = u_0\n",
    "    \\end{cases}\n",
    "$$\n",
    "where $u$ is typically a function $u : \\mathbb{R} \\rightarrow \\mathbb{R}^d$ with $d$ the dimension of the problem (for $d = 2$ we may also find $u : \\mathbb{R} \\rightarrow \\mathbb{C}$).\n",
    "In this first TD we will review some of the basic methods to solve Ordinary Differential Equations (ODEs) before studying one of the most famous chaotic dynamical system: the Lorenz system."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 1: The Exponential system\n",
    "\n",
    "We want to study in this part several linear one-stage methods to solve the following system on $[0,T]$:\n",
    "$$\\begin{cases}\n",
    "    u^\\prime=\\alpha u \\\\\n",
    "    u (0) = u_0\n",
    "\\end{cases}\n",
    "$$\n",
    "with $u : [0,T] \\rightarrow \\mathbb{C}$, $\\alpha \\in \\mathbb{C}$ and $u_0 \\in \\mathbb{C}$.\n",
    "We discretize the interval $[0, T]$ into $N + 1$ equally spaced points $\\{0, \\Delta t, ..., N \\Delta t\\}$ defining the time step $\\Delta t$ as $\\Delta t = \\frac{T}{N}$. The solution of this Cauchy problem is analytical and we will note $u(t) = u_0 e^{\\alpha t}$ its unique solution.\n",
    "\n",
    "We will make a clear distinction between the exact solution $u$ and its numerical approximations $u_{\\text{num},*}$ where the * character refers to the selected numerical scheme.\n",
    "\n",
    "More precisely, $u_{\\text{num},*} = (u_{\\text{num},*}^{(n)})_{n\\in\\{ 0 ,..., N\\}}$ will be a finite sequence of N + 1 (potentially complex) numbers such that $u_{\\text{num},*}^{(n)}$ approximates the value of the exact solution $u (n \\Delta t)$.\n",
    "\n",
    "In the following, Python codes will only need *numpy* and *matplotlib* packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# To draw matplotlib plots within this notebook.\n",
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will also use the following parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = 10.0\n",
    "N = 100\n",
    "dt = T / N"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**1) Forwards Euler:** We recall that Forwards Euler method (FE) is a first order method which computes the value of $u (t + \\Delta t)$ given $u(t)$ using the following Taylor expansion:\n",
    "$$\n",
    "    u (t + \\Delta t) = u(t) + \\Delta t~f (u(t)) + \\mathcal{O}(\\Delta t ^2)\n",
    "$$\n",
    "Express $u_{\\text{num},\\text{FE}}^{(n + 1)}$ as a function of $u_{\\text{num},\\text{FE}}^{(n)}$. This relation is the reccurence relation of Forwards Euler method.\n",
    "\n",
    "Then, what is the expression of $u_{\\text{num},\\text{FE}}^{(n)}$ as a function of $u_{\\text{num},\\text{FE}}^{(0)} = u_0$?\n",
    "\n",
    "**2)** We define $q = \\alpha \\Delta t$ and pick $u_0 = 1$. Program a Forwards Euler based solving method (using the recurrence relation) of the exponential system on $[0,T]$ and plot the numerical solution and the exact solution for several real values of $q$ (both positive and negative values). Can you observe some numerical instabilities?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**3) Backwards Euler / Crank-Nicolson:** Same questions for Backwards Euler method and Crank-Nicolson method. We recall that Backwards Euler method relies on the following expansion:\n",
    "$$\n",
    "    u (t + \\Delta t) = u(t) + \\Delta t~f(u(t + \\Delta t)) + \\mathcal{O}(\\Delta t ^2)\n",
    "$$\n",
    "and that Crank-Nicolson relies on the following expansion:\n",
    "$$\n",
    "    u (t + \\Delta t) = u(t) + \\frac{\\Delta t}{2}\\left(f\\left(u\\left(t\\right)\\right) + f\\left(u\\left(t + \\Delta t\\right)\\right)\\right) + \\mathcal{O}(\\Delta t ^3)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**4)** Compare on this example the three previous methods by plotting the three numerical solutions next to the exact solution. We will use a value of $q$ which is included in the absolute stability region of every method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**5)** Show theoretically that Crank-Nicolson scheme is a second-order method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2: The Harmonic Oscillator\n",
    "\n",
    "We now want to solve the following Cauchy problem on $[0,T]$:\n",
    "$$\\begin{cases}\n",
    "    u^{\\prime\\prime} +\\omega^2 u =0\\\\\n",
    "    u(0) = u_0\\\\\n",
    "    u^\\prime(0) = u^\\prime_0\n",
    "\\end{cases}\n",
    "$$\n",
    "where $u : [0,T] \\rightarrow \\mathbb{R}$, $\\omega\\in\\mathbb{R}$, $u_0 \\in \\mathbb{R}$ and $u^\\prime_0 \\in \\mathbb{R}$.\n",
    "\n",
    "**1)** Introducing $v = \\frac{u^\\prime}{\\omega}$, show that this sytem is equivalent to a complex exponential system such as defined in the previous part. We will express $\\alpha$ in terms of $\\omega$.\n",
    "\n",
    "**2)** Solve numerically this system for $\\omega = 1$ using Forwards Euler method, Backwards Euler method, Crank-Nicolson method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**3)** Plot the previous numerical solutions in the $(u, v)$ plan. If we want to enforce energy conservation in this system, which numerical method should we use? Explain theoretically what we observe following the selected numerical method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 3: A chaotic dynamical system: the Lorenz system\n",
    "\n",
    "We are now interested in the Lorenz system, first studied by the mathematician and meteorologist Edward Lorenz to model the atmospheric convection (we will refer to the following [Wikipedia page](https://en.wikipedia.org/wiki/Lorenz_system)). The Lorenz system is defined as follows:\n",
    "$$\\begin{cases}\n",
    "x^\\prime = \\sigma(y - x)\\\\\n",
    "y^\\prime = x (\\rho - z) - y\\\\\n",
    "z^\\prime = xy - \\beta z\n",
    "\\end{cases}\n",
    "$$\n",
    "where $x, y, z : \\mathbb{R} \\rightarrow \\mathbb{R}$ and $\\sigma, \\rho, \\beta \\in \\mathbb{R}$ are some numerical parameters. Lorenz first used the following values of the parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rho = 28.0\n",
    "sigma = 10.0\n",
    "beta = 8/3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to solve the Lorenz system using Runge-Kutta multistage methods. We will call RK2, the following Runge-Kutta method:\n",
    "$$\n",
    "    U_1 = u(t)\\\\\n",
    "    U_2 = u(t) + \\Delta t~f(U_1)\\\\\n",
    "    u(t + \\Delta t) = u(t) + \\frac{\\Delta t}{2}(f(U_1) + f(U_2)) + \\mathcal{O}(\\Delta t^3)\n",
    "$$\n",
    "\n",
    "**1)** Show theoretically that RK2 is a second-order method.\n",
    "\n",
    "**2)** Solve numerically the Lorenz system on the interval $[0,T]$ using RK2 method with the following initial values: $x (0) = y (0) = z (0) = 1$. We will use at least $T = 100$ and $N = 10000$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**3)** Plot the solution in the 3D-space $(x,y,z)$. We can use *[mpl_toolkits.mplot3d](https://matplotlib.org/mpl_toolkits/mplot3d/tutorial.html)* package for 3D plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We call this figure the *Lorenz attractor*.\n",
    "The Lorenz system exhibits chaotic behaviors for some specific values of its parameters: this means that solutions are highly senstive to initial conditions.\n",
    "\n",
    "**4)** For the same values of the parameters, slightly change the initial conditions and observe this chaotic behavior."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n"
   ]
  }
 ],
 "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
