{
 "metadata": {
  "name": "",
  "signature": "sha256:72fd3024fb1fdbced56d8d02e1572f31cb7d3dc8c980ed59c35f68d62a3b3d71"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Featured Recipe #5: Simulating a Partial Differential Equation: reaction-diffusion systems and Turing patterns"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "> This is a featured recipe from the [**IPython Cookbook**](http://ipython-books.github.io/), the definitive guide to **high-performance scientific computing** and **data science** in Python."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Partial Differential Equations** (PDEs) describe the evolution of dynamical systems involving both time and space. Examples in physics include sound, heat, electromagnetism, fluid flow, elasticity, among others. Examples in biology include tumor growth, population dynamics, and epidemic propagations.\n",
      "\n",
      "PDEs are hard to solve analytically. Therefore, PDEs are often studied via numerical simulations.\n",
      "\n",
      "In this featured recipe, we will illustrate how to simulate a reaction-diffusion system described by a PDE called the **Fitzhugh\u2013Nagumo equation**. A reaction-diffusion system models the evolution of one or several variables subject to two processes: reaction (transformation of the variables into each other) and diffusion (spreading across a spatial region). Some chemical reactions may be described by this type of model, but there are other applications in physics, biology, ecology, and other disciplines.\n",
      "\n",
      "Here, we simulate a system that has been proposed by Alan Turing as a model of animal coat pattern formation. Two chemical substances influencing skin pigmentation interact according to a reaction-diffusion model. This system is responsible for the formation of patterns that are reminiscent of the pelage of zebras, jaguars, and giraffes.\n",
      "\n",
      "We will simulate this system with the finite difference method. This method consists of discretizing time and space, and replacing the derivatives by their discrete equivalents."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## How to do it..."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "1. Let's import the packages."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "%matplotlib inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "2. We will simulate the following system of partial differential equations on the domain $E=[-1,1]^2$:"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "style": "latex"
     },
     "source": [
      "\\begin{align*}\n",
      "\\frac{\\partial u}{\\partial t} &= a \\Delta u + u - u^3 - v + k\\\\\n",
      "\\tau\\frac{\\partial v}{\\partial t} &= b \\Delta v + u - v\\\\\n",
      "\\end{align*}"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The variable $u$ represents the concentration of a substance favoring skin pigmentation, whereas $v$ represents another substance that reacts with the first and impedes pigmentation.\n",
      "\n",
      "At initialization time, we assume that $u$ and $v$ contain independent random numbers on every grid point. Besides, we take **Neumann boundary conditions**: we require the spatial derivatives of the variables with respect to the normal vectors to be null on the boundaries of the domain $E$.\n",
      "\n",
      "Let's define the four parameters of the model."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "a = 2.8e-4\n",
      "b = 5e-3\n",
      "tau = .1\n",
      "k = -.005"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "3. We discretize time and space. The following condition ensures that the discretization scheme we use here is stable:\n",
      "\n",
      "$$dt \\leq \\frac{dx^2}{2}$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "size = 100  # size of the 2D grid\n",
      "dx = 2./size  # space step"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "T = 10.0  # total time\n",
      "dt = .9 * dx**2/2  # time step\n",
      "n = int(T/dt)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "4. We initialize the variables $u$ and $v$. The matrices $U$ and $V$ contain the values of these variables on the vertices of the 2D grid. These variables are initialized with a uniform noise between $0$ and $1$."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "U = np.random.rand(size, size)\n",
      "V = np.random.rand(size, size)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "5. Now, we define a function that computes the discrete Laplace operator of a 2D variable on the grid, using a five-point stencil finite difference method. This operator is defined by:\n",
      "\n",
      "$$\\Delta u(x,y) \\simeq \\frac{u(x+h,y)+u(x-h,y)+u(x,y+h)+u(x,y-h)-4u(x,y)}{dx^2}$$\n",
      "\n",
      "We can compute the values of this operator on the grid using vectorized matrix operations. Because of side effects on the edges of the matrix, we need to remove the borders of the grid in the computation."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def laplacian(Z):\n",
      "    Ztop = Z[0:-2,1:-1]\n",
      "    Zleft = Z[1:-1,0:-2]\n",
      "    Zbottom = Z[2:,1:-1]\n",
      "    Zright = Z[1:-1,2:]\n",
      "    Zcenter = Z[1:-1,1:-1]\n",
      "    return (Ztop + Zleft + Zbottom + Zright - 4 * Zcenter) / dx**2"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "6. Now, we simulate the system of equations using the finite difference method. At each time step, we compute the right-hand sides of the two equations on the grid using discrete spatial derivatives (Laplacians). Then, we update the variables using a discrete time derivative."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# We simulate the PDE with the finite difference method.\n",
      "for i in range(n):\n",
      "    # We compute the Laplacian of u and v.\n",
      "    deltaU = laplacian(U)\n",
      "    deltaV = laplacian(V)\n",
      "    # We take the values of u and v inside the grid.\n",
      "    Uc = U[1:-1,1:-1]\n",
      "    Vc = V[1:-1,1:-1]\n",
      "    # We update the variables.\n",
      "    U[1:-1,1:-1], V[1:-1,1:-1] = \\\n",
      "        Uc + dt * (a * deltaU + Uc - Uc**3 - Vc + k), \\\n",
      "        Vc + dt * (b * deltaV + Uc - Vc) / tau\n",
      "    # Neumann conditions: derivatives at the edges\n",
      "    # are null.\n",
      "    for Z in (U, V):\n",
      "        Z[0,:] = Z[1,:]\n",
      "        Z[-1,:] = Z[-2,:]\n",
      "        Z[:,0] = Z[:,1]\n",
      "        Z[:,-1] = Z[:,-2]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "7. Finally, we display the variable $u$ after a time $T$ of simulation."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.imshow(U, cmap=plt.cm.copper, extent=[-1,1,-1,1]);\n",
      "plt.xticks([]); plt.yticks([]);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "![Turing patterns](images/turing.jpg)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Whereas the variables when completely random at initialization time, we observe the formation of patterns after a sufficiently long simulation time."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## How it works..."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's explain how the finite difference method allowed us to implement the update step. We start from the following system of equations:"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "style": "latex"
     },
     "source": [
      "\\begin{align*}\n",
      "\\frac{\\partial u}{\\partial t}(t;x,y) &= a \\Delta u(t;x,y) + u(t;x,y) - u(t;x,y)^3 - v(t;x,y) + k\\\\\n",
      "\\tau\\frac{\\partial v}{\\partial t}(t;x,y) &= b \\Delta v(t;x,y) + u(t;x,y) - v(t;x,y)\\\\\n",
      "\\end{align*}"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We first use the following scheme for the discrete Laplace operator:"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$$\\Delta u(x,y) \\simeq \\frac{u(x+h,y)+u(x-h,y)+u(x,y+h)+u(x,y-h)-4u(x,y)}{dx^2}$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We also use this scheme for the time derivative of $u$ and $v$:"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$$\\frac{\\partial u}{\\partial t}(t;x,y) \\simeq \\frac{u(t+dt;x,y)-u(t;x,y)}{dt}$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We end up with the following iterative update step:"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "style": "latex"
     },
     "source": [
      "\\begin{align*}\n",
      "u(t+dt;x,y) &= u(t;x,y) + dt \\left( a \\Delta u(t;x,y) + u(t;x,y) - u(t;x,y)^3 - v(t;x,y) + k \\right)\\\\\n",
      "v(t+dt;x,y) &= v(t;x,y) + \\frac{dt}{\\tau} \\left( b \\Delta v(t;x,y) + u(t;x,y) - v(t;x,y) \\right)\\\\\n",
      "\\end{align*}"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Here, our Neumann boundary conditions state that the spatial derivatives with respect to the normal vectors are null on the boundaries of the domain E:"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$$\\forall w \\in \\{u, v\\}, \\, \\forall t \\geq 0, \\, \\forall x, y \\in \\partial E, \\quad \\frac{\\partial w}{\\partial x}(t; -1, y) = \\frac{\\partial w}{\\partial x}(t; 1, y) = \\frac{\\partial w}{\\partial y}(t; x, -1) = \\frac{\\partial w}{\\partial y}(t; x, 1) = 0$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We implement these boundary conditions by duplicating values in matrices $U$ and $V$ on the edges (see code)."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In order to ensure that our numerical scheme converges to a numerical solution that is close to the actual (unknown) mathematical solution, the stability of the scheme needs to be ascertained. One can show that, here, a sufficient condition for the stability is:"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$$dt \\leq \\frac{dx^2}{2}$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## There's more..."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Here are further references on partial differential equations, reaction-diffusion systems, and numerical simulations of those systems."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "* [Partial Differential Equations](http://en.wikipedia.org/wiki/Partial_differential_equation),\n",
      "* [Reaction-diffusion systems](http://en.wikipedia.org/wiki/Reaction%E2%80%93diffusion_system),\n",
      "* [Fitzhugh-Nagumo system](http://en.wikipedia.org/wiki/FitzHugh%E2%80%93Nagumo_equation),\n",
      "* [Neumann boundary conditions](http://en.wikipedia.org/wiki/Neumann_boundary_condition),\n",
      "* [Von Neumann stability analysis](http://en.wikipedia.org/wiki/Von_Neumann_stability_analysis)."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "> This was a featured recipe from the [IPython Cookbook](http://ipython-books.github.io/), by [Cyrille Rossant](http://cyrille.rossant.net), Packt Publishing, 2014."
     ]
    }
   ],
   "metadata": {}
  }
 ]
}