{
"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
}