{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# TD7: Searching for steady states\n",
"\n",
"The goal of this TD is to get familiar with some techniques to find the steady states of a dynamical system. We will be interested in the following equation:\n",
"$$\n",
" \\dot{x} = -x^3 + 5x + \\mu\n",
"$$\n",
"This differential equation depends on the real parameter $\\mu$. In this case, the steady states are 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",
"In the following, we will need the usual Python packages:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib inline"
]
},
{
"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 solve the multiple Cauchy problems associated with a given set of initial conditions. Therefore, in this part we will solve 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 \\in \\mathbb{R}$.\n",
"\n",
"We will use a Runge-Kutta method to solve this system. 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 (x (t))\\\\\n",
" X_2 &= \\Delta t~f_\\mu (x(t) + \\frac{X_1}{2})\\\\\n",
" X_3 &= \\Delta t~f_\\mu (x(t) + \\frac{X_2}{2})\\\\\n",
" X_4 &= \\Delta t~f_\\mu (x(t) + X_3)\\\\\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]$ where $T$ is big enough so that the numerical solution does not evolve significantly anymore. In the following, $T$ will be defined so that $|f_\\mu(x(T))|$ is smaller than a given threshold $\\sigma$. Obviously, the numerical solution has to converge so that $T$ can be defined.\n",
"\n",
"**1)** Build a $RK4$ method 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": [
"We now want to build a **bifurcation diagram**, which represents the steady states $x_\\infty$ as a function of the parameter $\\mu$ (there can be several steady states for a given value of $\\mu$ as mentionned above).\n",
"\n",
"**2)** Plot the bifurcation diagram of this dynamical system. For example, we can use values of $\\mu$ that range from $-10$ to $10$, and for each of these values, use values of $x_0$ that range from $-5$ to $5$. We will need to define a reasonable threshold $\\sigma$ as mentionned above."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**3)** Can we find every steady state 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 unstable steady states. **Newton's method** is a technique that can find zeros of a function independently of their stability within 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 impact the limit value of the previous sequence. We will need to try several initial conditions to find every zero of $f$.\n",
"\n",
"**4)** 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": [
"**5)** 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 a suitable $x_0$ and for a given value of $\\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": [
"**6)** Plot a new bifurcation diagram of our dynamical system using Newton's method. Is it 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 find these roots starting from a given set of initial values within 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 relies on the following recurrent sequence:\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",
"**7)** Calculate the expressions of $f$ and $\\text{D}f$.\n",
"\n",
"**8)** Perform Newton's method on $f$ for a dense enough 2D grid of initial conditions and for a fixed 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": [
"import scipy.linalg as la\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**9)** 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 argument of the final complex term of the Newton sequence. Does Newton's method always converge?"
]
},
{
"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.6.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}