{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TD7: Searching for steady states\n",
    "\n",
    "The goal of this TD is to get familiar with a few techniques that allow you to find the steady states of a dynamical system from a numerical point of view.\n",
    "\n",
    "We will be interested in the following differential equation:\n",
    "$$\n",
    "    \\dot{x} = -x^3 + 5x + \\mu,\n",
    "$$\n",
    "that depends on a real parameter $\\mu$. Here, the steady states are simply the roots of the polynomial $f_\\mu(x) = -x^3 + 5x + \\mu$. Depending on the value of $\\mu$, there can be one, two or three real roots.\n",
    "\n",
    "We will need the usual Python packages:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "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": [
    "## Part 1: Runge-Kutta approach: an uncomplete bifurcation diagram\n",
    "\n",
    "A first approach to find the steady states of a dynamical system is to draw some various initial conditions and integrate the corresponding Cauchy problems for a sufficient long time. This time has to be long enough to reach a steady state, at least approximately.\n",
    "\n",
    "*Note that in general, a solution of a dynamical system is not guaranteed to converge, but for our particular system any solution will converge to a steady state (can you see why?).*\n",
    "\n",
    "We thus want to solve multiple Cauchy problems such as:\n",
    "$$\n",
    "\\begin{cases}\n",
    "    \\dot{x} &= f_\\mu (x)\\\\\n",
    "    x (0) &= x_0\n",
    "\\end{cases}\n",
    "$$\n",
    "with $x_0$ chosen from a set of relevant values $\\{x_0^{(1)}, ..., x_0^{(m)}\\}$.\n",
    "\n",
    "We will use a Runge-Kutta method to integrate these systems. We recall the RK4 method (Runge-Kutta method of fourth order) that computes the value of the numerical solution $x (t)$ as follows:\n",
    "$$\n",
    "\\begin{align}\n",
    "    X_1 &= \\Delta t~f_\\mu\\left(x (t)\\right)\\\\\n",
    "    X_2 &= \\Delta t~f_\\mu\\left(x(t) + \\frac{X_1}{2}\\right)\\\\\n",
    "    X_3 &= \\Delta t~f_\\mu\\left(x(t) + \\frac{X_2}{2}\\right)\\\\\n",
    "    X_4 &= \\Delta t~f_\\mu\\left(x(t) + X_3\\right)\\\\\n",
    "    x (t + \\Delta t) &= x (t) + \\frac1{6}(X_1 + 2X_2 + 2X_3 + X_4)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "For a given initial condition $x_0$, we will solve the Cauchy problem on $[0, T]$, with $T$ big enough to claim that the numerical solution has approximately converged. In practice, we can use the following criterion to define $T$: let $T$ be the smallest positive number such as $|f_\\mu(x(T))| \\leq \\sigma$, with $\\sigma$ an arbitrary threshold .\n",
    "\n",
    "**1)** Define a $rk4$ function which returns the value of $x(t+\\Delta t)$ in function of $x(t)$, $f$, and $\\Delta t$ using the previous Runge-Kutta method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**2)** For given values of $\\mu$, $x_0$ and $\\sigma$, integrate the system on $[0,T]$. We will plot the corresponding solution. Give an approximate value of the steady state $x_\\infty$ reached by this solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we know how to integrate our system and to approximate its steady states, we want to describe the dependence of the steady states to the $\\mu$ parameter. We want to build a **bifurcation diagram**, which represents the steady states $x_\\infty$ as a function of $\\mu$ (there can be several steady states for a given value of $\\mu$ as mentioned above).\n",
    "\n",
    "**3)** Plot the bifurcation diagram of this dynamical system. We can typically probe values of $\\mu$ that range from $-10$ to $10$, and for each of these values, probe values of $x_0$ that range from $-5$ to $5$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**4)** Can we find all the steady states of our system for a given value of $\\mu$ with this method?\n",
    "\n",
    "## Part 2: Newton's method: a complete bifurcation diagram\n",
    "\n",
    "We need another technique if we want to find the unstable steady states. **Newton's method** is a general technique that find the zeros of a function independently of how stable these zeros are in the context of a dynamical system.\n",
    "\n",
    "For a function $f$ defined over the real numbers, starting from a real $x_0$, Newton's method converges towards a zero of $f$ through the following recurrent sequence:\n",
    "$$\n",
    "    x_{n+1} = x_n - \\frac{f(x_n)}{f^\\prime(x_n)}\n",
    "$$\n",
    "Obviously the choice of $x_0$ will be determinant to the limit value of the previous sequence. We will need to try several initial conditions to find every zero of $f$.\n",
    "\n",
    "**5)** Build a *newton* routine that computes the value of $x_{n+1}$ as a function of $f$, $f^\\prime$ and $x_n$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**6)** Observe Newton's method quadratic convergence to a zero of $f_\\mu$ polynomial printing $\\left(f_\\mu\\left(x_n\\right)\\right)_n$ first terms for given values of $x_0$ and $\\mu$. Do we always have a quadratic convergence of this method?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**7)** Plot a new bifurcation diagram of our dynamical system using Newton's method. Does it seem complete now?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bonus part: The Newton fractal\n",
    "\n",
    "The [Newton fractal](https://en.wikipedia.org/wiki/Newton_fractal) is a fractal figure which stems from the application of Newton's method to the following complex polynomial:\n",
    "$$\n",
    "    f(z) = z^3 - 1\n",
    "$$\n",
    "The roots of the previous polynomial are well known, these are the third roots of unity: $\\{1, e^{\\frac{2i\\pi}{3}}, e^{-\\frac{2i\\pi}{3}}\\}$. Still, it is interesting to approximate these roots with a Newton's method applied to a regulard grid of initial values in the complex plane as spectacular properties of convergence come out from Newton's algorithm.\n",
    "\n",
    "In the following, we will identify $\\mathbb{C}$ to $\\mathbb{R}^2$ for pedagogical purposes. Therefore Newton's method is to be applied to a function $f:\\mathbb{R}^2\\rightarrow\\mathbb{R}^2$ defined as follows:\n",
    "$$\n",
    "    f (x, y) = \\left( \\begin{array}{c} f_1 (x,y) \\\\ f_2 (x,y)\\end{array}\\right) = \\left( \\begin{array}{c} \\Re\\left(\\left(x+iy\\right)^3 - 1\\right) \\\\ \\Im\\left(\\left(x+iy\\right)^3 - 1\\right)\\end{array}\\right)\n",
    "$$\n",
    "where $\\Re$ and $\\Im$ respectively denotes the real and imaginary part of a given complex number.\n",
    "\n",
    "For a vector-valued function such as above $f$ function, Newton's method can be written as follows:\n",
    "$$\n",
    "    X_{n+1} = X_n - \\left[\\text{D}f(X_n)\\right]^{-1}f(X_n)\n",
    "$$\n",
    "where $X_n$ is a vector (of dimension 2 for above $f$ function) and $\\text{D}f(X_n)$ is the Jacobian matrix of $f$ at $X_n$. We recall the general expression of the Jacobian matrix for a function $f:\\mathbb{R}^m\\rightarrow\\mathbb{R}^m$:\n",
    "$$\n",
    "    \\text{D}f(X) = \\left(\\frac{\\partial f_i}{\\partial x_j}(X)\\right)_{i, j\\in\\{1,...,m\\}}\n",
    "$$\n",
    "\n",
    "**8)** What are the analytical expressions of $f$ and $\\text{D}f$?\n",
    "\n",
    "**9)** Perform Newton's method on $f$ for a dense enough 2D grid of initial conditions and for a uniform number of Newton iterations. We will use [scipy.linalg.inv](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.inv.html) function to invert the Jacobian matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**10)** Plot using [pcolor](https://matplotlib.org/api/_as_gen/matplotlib.pyplot.pcolor.html) function for each point of the previous grid of initial conditions the complex argument of the final term of the Newton sequence (we can use [numpy.angle](https://docs.scipy.org/doc/numpy/reference/generated/numpy.angle.html) function to do that). You should observe a fractal figure called the *Newton fractal*. Does Newton's method always converge in this case?"
   ]
  },
  {
   "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
}
