reaction-diffusion_demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<table>\n",
    " <tr align=left><td><img align=left src=\"https://i.creativecommons.org/l/by/4.0/88x31.png\">\n",
    " <td>Text provided under a Creative Commons Attribution license, CC-BY. All code is made available under the FSF-approved MIT license. (c) Kyle T. Mandli</td>\n",
    "</table>\n",
    "Based on an example from https://github.com/ketch/finite-difference-course"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "from __future__ import print_function\n",
    "\n",
    "import numpy\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.sparse as sparse\n",
    "import scipy.sparse.linalg as linalg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mixed Equations In-Class Project"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the reaction-diffusion PDE\n",
    "$$\\begin{aligned}\n",
    "    u_t &= \\sigma D_1 \\nabla^2 u + f(u, v) \\\\\n",
    "    v_t &= \\sigma D_2 \\nabla^2 v + g(u, v)\n",
    "\\end{aligned}$$\n",
    "in two-dimensions (i.e. $\\nabla^2 u = u_{xx} + u_{yy}$) and with the source terms\n",
    "$$\\begin{aligned}\n",
    "    f(u,v) &= \\alpha u (1 - \\tau_1 v^2) + v (1 - \\tau_2 u) \\\\\n",
    "    g(u,v) &= \\beta v + \\alpha \\tau_1 u v^2 + u (\\gamma + \\tau_2 v).\n",
    "\\end{aligned}$$\n",
    "These equations with the appropriate parameters $\\sigma, D_1, D_2, \\alpha, \\beta, \\tau_1, \\tau_2, \\gamma$ can be used to study emergent patterns from seemingly random initial data which we will investigate numerically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f_reaction(U, V, sigma, tau_1, tau_2, alpha, beta, gamma):\n",
    "    return alpha * U * (1.0 - tau_1 * V**2) + V * (1.0 - tau_2 * U)\n",
    "\n",
    "def g_reaction(U, V, sigma, tau_1, tau_2, alpha, beta, gamma):\n",
    "    return beta * V * (1.0 + alpha * tau_1 / beta * U * V) + U * (gamma + tau_2 * V)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spatial Derivative Discretization\n",
    "\n",
    "Let's consider the above PDEs on the a square domain $\\Omega = [-1, 1] \\times [-1, 1]$ with periodic boundary conditions.  First write a function that uses a five-point stencil to represent the Laplacian operator in 2d and returns the appropriate sparse matrix reprsentation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def laplacian_discretization(m):\n",
    "    \"\"\"Constructs a sparse matrix that discretizes the 2d Laplacian\n",
    "    \n",
    "    Uses a five-point stencil and periodic boundary conditions.\n",
    "    \"\"\"\n",
    "    delta_x = 2.0 / (m + 1)\n",
    "    \n",
    "    # Primary discretization\n",
    "    e = numpy.ones(m)\n",
    "    T = sparse.spdiags([e, -4.0 * e, e], [-1, 0, 1], m, m)\n",
    "    S = sparse.spdiags([e, e], [-1, 1], m, m)\n",
    "    I = sparse.eye(m)\n",
    "    A = sparse.kron(I, T) + sparse.kron(S, I)\n",
    "    \n",
    "    # Construct periodic BCs\n",
    "    e = numpy.ones(m**2)\n",
    "    A_periodic = sparse.spdiags([e, e],[m - m**2, m**2 - m], m**2, m**2).tolil()\n",
    "    # Left & right BCs:\n",
    "    for i in range(m):\n",
    "        A_periodic[i * m, (i + 1) * m - 1] = 1.0\n",
    "        A_periodic[(i + 1) * m - 1, i * m] = 1.0\n",
    "    \n",
    "    # Combine two matrices\n",
    "    A = A + A_periodic\n",
    "    A /= delta_x**2\n",
    "    A = A.todia()\n",
    "    \n",
    "    return A\n",
    "\n",
    "A = laplacian_discretization(4)\n",
    "plt.spy(A)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Time Stepping\n",
    "\n",
    "First let's see if we can make a simple explicit method, in this case forward Euler, work for us.  We know this might not be such as great idea due to the diffusion term but maybe the reaction terms will be helpfull.\n",
    "\n",
    "First write a function that uses forward Euler to take a single time step to solve the equations of interest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def forward_euler_step(U, V, delta_t, A, sigma, f, g, D1=0.5, D2=1.0):\n",
    "    \"\"\"Take a single forward Euler step on the reaction-diffusion equation\"\"\"\n",
    "    \n",
    "    U_new = U + delta_t * (sigma * D1 * A * U + f(U, V))\n",
    "    V_new = V + delta_t * (sigma * D2 * A * V + g(U, V))\n",
    "    \n",
    "    return U_new, V_new"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now try to solve the PDE given the parameters\n",
    "$$\n",
    "    \\sigma = 0.0021, ~ \\tau_1 = 3.5, ~ \\tau_2 = 0.0, ~ \\alpha = 0.899, ~ \\beta=-0.91, ~\\gamma=-\\alpha\n",
    "$$\n",
    "with the default values of $D_1 = 0.5$ and $D_2 = 1.0$.  We will also take a random initial condition.\n",
    "\n",
    "Note what step-size we might need here.  For the two-dimensional heat equation we can show that forward Euler is going to require a step size of\n",
    "$$\n",
    "    \\Delta t \\leq \\frac{\\Delta x^2}{4 \\kappa}\n",
    "$$\n",
    "where now $\\kappa$ is the coefficient out in front of the Laplacian.  Here we will take the maximum of the coefficient in front of the Laplacians to remain stable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def forward_euler_coupled_solver(sigma, tau_1, tau_2, alpha, beta, gamma, D_1, D_2):\n",
    "    # Alias reaction functions with the above parameters\n",
    "    f = lambda U, V: f_reaction(U, V, sigma, tau_1, tau_2, alpha, beta, gamma)\n",
    "    g = lambda U, V: g_reaction(U, V, sigma, tau_1, tau_2, alpha, beta, gamma)\n",
    "\n",
    "    # Set up grid\n",
    "    m = 150\n",
    "    delta_x = 2.0 / m\n",
    "    x = numpy.linspace(-1.0, 1.0, m)\n",
    "    y = numpy.linspace(-1.0, 1.0, m)\n",
    "    Y, X = numpy.meshgrid(y, x)\n",
    "\n",
    "    # Initial data\n",
    "    U = numpy.random.randn(m, m) / 2.0\n",
    "    V = numpy.random.randn(m, m) / 2.0\n",
    "    fig = plt.figure()\n",
    "    axes = fig.add_subplot(1, 1, 1, aspect='equal')\n",
    "    plot = axes.pcolor(x, y, U, cmap=plt.get_cmap(\"viridis\"))\n",
    "    fig.colorbar(plot)\n",
    "\n",
    "    # Setup spatial discretization\n",
    "    U = U.reshape(-1)\n",
    "    V = V.reshape(-1)\n",
    "    A = laplacian_discretization(m)\n",
    "\n",
    "    # Time\n",
    "    t = 0.0\n",
    "    t_final = 300.0\n",
    "    delta_t = delta_x**2 / (5.0 * sigma)\n",
    "    num_steps = int(numpy.round(t_final / delta_t))\n",
    "\n",
    "    # Evolve in time\n",
    "    next_output_time = 0.0\n",
    "    for j in range(num_steps):\n",
    "        U, V = forward_euler_step(U, V, delta_t, A, sigma, f, g)\n",
    "        t += delta_t\n",
    "\n",
    "        if t >= next_output_time:\n",
    "            next_output_time += 50.0\n",
    "            U_output = U.reshape((m, m))\n",
    "\n",
    "            fig = plt.figure()\n",
    "            axes = fig.add_subplot(1, 1, 1, aspect='equal')\n",
    "            plot = axes.pcolor(x, y, U_output, cmap=plt.get_cmap(\"viridis\"))\n",
    "            fig.colorbar(plot)\n",
    "            axes.set_title(\"t = %s\" % t)\n",
    "\n",
    "    plt.show()\n",
    "    \n",
    "forward_euler_coupled_solver(sigma=0.0021, tau_1=3.5, tau_2=0, alpha=0.899, beta=-0.91, gamma=-0.899, D_1=0.5, D_2=1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implicit-Explicit Splitting\n",
    "\n",
    "The previous approach was clearly very slow so let's try applying one of our splitting techniques to the problem instead.  IMEX methods are actually pretty ideal for this case so let's try using backwards Euler for the stiff diffusion term and the forward Euler time step for the explicit reaction terms.\n",
    "\n",
    "Implicit:\n",
    "$$\\begin{aligned}\n",
    "    u_t &= \\sigma D_1 \\nabla^2 u \\\\\n",
    "    v_t &= \\sigma D_2 \\nabla^2 v\n",
    "\\end{aligned}$$\n",
    "\n",
    "Explicit:\n",
    "$$\\begin{aligned}\n",
    "    u_t &= f(u, v) \\\\\n",
    "    v_t &= g(u, v)\n",
    "\\end{aligned}$$\n",
    "\n",
    "Numerical method:\n",
    "$$\\begin{aligned}\n",
    "    U^\\ast &= U^n + \\Delta t \\sigma D_1 \\nabla^2 U^\\ast \\\\\n",
    "    V^\\ast &= V^n + \\Delta t \\sigma D_2 \\nabla^2 V^\\ast \\\\\n",
    "    U^{n+1} &= U^\\ast + \\Delta t f(U^\\ast, V^\\ast) \\\\\n",
    "    V^{n+1} &= V^\\ast + \\Delta t g(U^\\ast, V^\\ast) \\\\\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def backward_euler_diffusion_step(U, V, A, delta_t, sigma, D_1, D_2):\n",
    "    U = linalg.spsolve((sparse.eye(A.shape[0]) - delta_t * sigma * D_1 * A), U)\n",
    "    V = linalg.spsolve((sparse.eye(A.shape[0]) - delta_t * sigma * D_2 * A), V)\n",
    "    return U, V\n",
    "\n",
    "def forward_euler_reaction_step(U, V, delta_t, f, g):\n",
    "    U_new = U + delta_t * f(U, V)\n",
    "    V_new = V + delta_t * g(U, V)\n",
    "    return U_new, V_new"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def imex_solver(sigma, tau_1, tau_2, alpha, beta, gamma, D_1, D_2):\n",
    "    # Alias reaction functions with the above parameters\n",
    "    f = lambda U, V: f_reaction(U, V, sigma, tau_1, tau_2, alpha, beta, gamma)\n",
    "    g = lambda U, V: g_reaction(U, V, sigma, tau_1, tau_2, alpha, beta, gamma)\n",
    "\n",
    "    # Set up grid\n",
    "    m = 150\n",
    "    delta_x = 2.0 / m\n",
    "    x = numpy.linspace(-1.0, 1.0, m)\n",
    "    y = numpy.linspace(-1.0, 1.0, m)\n",
    "    Y, X = numpy.meshgrid(y, x)\n",
    "\n",
    "    # Initial data\n",
    "    U = numpy.random.randn(m, m) / 2.0\n",
    "    V = numpy.random.randn(m, m) / 2.0\n",
    "    fig = plt.figure()\n",
    "    axes = fig.add_subplot(1, 1, 1, aspect='equal')\n",
    "    plot = axes.pcolor(x, y, U, cmap=plt.get_cmap(\"viridis\"))\n",
    "    fig.colorbar(plot)\n",
    "\n",
    "    # Setup spatial discretization\n",
    "    U = U.reshape(-1)\n",
    "    V = V.reshape(-1)\n",
    "    A = laplacian_discretization(m)\n",
    "\n",
    "    # Time\n",
    "    t = 0.0\n",
    "    t_final = 30.0\n",
    "    delta_t = delta_x / (10.0 * sigma)\n",
    "    num_steps = int(numpy.round(t_final / delta_t))\n",
    "\n",
    "    # Evolve in time\n",
    "    next_output_time = 0.0\n",
    "    for j in range(num_steps):\n",
    "        U, V = backward_euler_diffusion_step(U, V, A, delta_t, sigma, D_1, D_2)\n",
    "        U, V = forward_euler_step(U, V, delta_t, A, sigma, f, g)\n",
    "        t += delta_t\n",
    "\n",
    "        if t >= next_output_time:\n",
    "            next_output_time += 5.0\n",
    "            U_output = U.reshape((m, m))\n",
    "\n",
    "            fig = plt.figure()\n",
    "            axes = fig.add_subplot(1, 1, 1, aspect='equal')\n",
    "            plot = axes.pcolor(x, y, U_output, cmap=plt.get_cmap(\"viridis\"))\n",
    "            fig.colorbar(plot)\n",
    "            axes.set_title(\"t = %s\" % t)\n",
    "\n",
    "    plt.show()\n",
    "    \n",
    "# Parameters\n",
    "imex_solver(sigma=0.0021, tau_1=3.5, tau_2=0, alpha=0.899, beta=-0.91, gamma=-0.899, D_1=0.5, D_2=1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try playing with the input parameters and see what kind of behavior you see."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma=0.0045; tau1=2.02; tau2=0.; alpha=2.0; beta=-0.91; gamma=-alpha;\n",
    "sigma=0.0005; tau1=2.02; tau2=0.; alpha=2.0; beta=-0.91; gamma=-alpha;\n",
    "sigma=0.0021; tau1=3.5; tau2=0; alpha=0.899; beta=-0.91; gamma=-alpha;\n",
    "sigma=0.0045; tau1=0.02; tau2=0.2; alpha=1.9; beta=-0.85; gamma=-alpha;\n",
    "sigma=0.0001; tau1=0.02; tau2=0.2; alpha=0.899; beta=-0.91; gamma=-alpha;\n",
    "sigma=0.0045; tau1=0.02; tau2=0.2; alpha=1.9; beta=-0.91; gamma=-alpha;"
   ]
  }
 ],
 "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.4"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}