8_A_reaction-diffusion_model_for_pattern_formation.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:e5673bf2886f9cfd2053a9affdb459d30869ed7b879d0fc6e4e9c0c84109e9ec"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline\n",
      "import numpy as np\n",
      "import scipy.optimize\n",
      "import scipy.sparse\n",
      "import matplotlib.pyplot as plt\n",
      "from IPython.display import clear_output\n",
      "import time"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# The model\n",
      "\n",
      "In this notebook we'll solve a system of reaction-diffusion PDEs in two dimensions:\n",
      "\n",
      "\\begin{align}\n",
      "u_t & = \\delta D_1 \\nabla^2 u + f(u,v) \\\\\\\\\n",
      "v_t & = \\delta D_2 \\nabla^2 v + g(u,v)\n",
      "\\end{align}\n",
      "where $\\nabla^2 u = u_{xx} + u_{yy}$ denotes the Laplacian and $f,g$ represent reaction terms.\n",
      "\n",
      "For simplicity, we'll consider the square domain $[-1,1]\\times[-1,1]$ with periodic boundary conditions; i.e., the conditions on $u(x,y,t)$ are\n",
      "\\begin{align}\n",
      "u(-1,y,t) & = u(1,y,t) \\\\\\\\\n",
      "u(x,-1,t) & = u(x,1,t)\n",
      "\\end{align}\n",
      "with corresponding conditions on $v$."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Reaction equations\n",
      "The reaction terms we will use are\n",
      "\\begin{align}\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{align}\n",
      "Once we have things working, we'll investigate how the constants $D_1, D_2, \\alpha, \\beta, \\tau_1, \\tau_2, \\gamma$ change the kinds of solutions obtained."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "delta=0.0045; tau1=0.02; tau2=0.2; alpha=0.899; beta=-0.91; gamma=-alpha;\n",
      "\n",
      "def f(u,v):\n",
      "    return alpha*u*(1-tau1*v**2) + v*(1-tau2*u);\n",
      "\n",
      "def g(u,v):\n",
      "    return beta*v*(1+alpha*tau1/beta*u*v) + u*(gamma+tau2*v);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Diffusion terms\n",
      "Next let's write a discretization of the diffusion terms.  We'll use the 5-point Laplacian approximation:\n",
      "\n",
      "\\begin{align}\n",
      "\\nabla^2 u(x_i,y_j) & \\approx \\frac{U_{i-1,j} - 2U_{ij} + U_{i+1,j}}{\\Delta x}\n",
      "                           +  \\frac{U_{i,j-1} - 2U_{ij} + U_{i,j+1}}{\\Delta y}\n",
      "\\end{align}\n",
      "For simplicity we'll assume $\\Delta x = \\Delta y = h$.  Examine the code below carefully, noting how the periodic boundary conditions are imposed."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def five_pt_laplacian_sparse_periodic(m,a,b):\n",
      "    \"\"\"Construct a sparse matrix that applies the 5-point laplacian discretization\n",
      "       with periodic BCs on all sides.\"\"\"\n",
      "    e=np.ones(m**2)\n",
      "    e2=([1]*(m-1)+[0])*m\n",
      "    e3=([0]+[1]*(m-1))*m\n",
      "    h=(b-a)/(m+1)\n",
      "    A=scipy.sparse.spdiags([-4*e,e2,e3,e,e],[0,-1,1,-m,m],m**2,m**2)\n",
      "    # Top & bottom BCs:\n",
      "    A_periodic = scipy.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.\n",
      "        A_periodic[(i+1)*m-1,i*m] = 1.\n",
      "    A = A + A_periodic\n",
      "    A/=h**2\n",
      "    A = A.todia()\n",
      "    return A"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "A = five_pt_laplacian_sparse_periodic(4,-1.,1.)\n",
      "plt.spy(A)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Forward Euler time stepping\n",
      "Since we are primarily interested in the final, steady-state solution, we can use a first-order accurate method in time.  Let's start by trying Euler's method:\n",
      "$$\\begin{align}\n",
      "U^{n+1} & = U^n + k ( \\delta D_1 \\nabla^2_h U^n  + f(U^n,V^n)) \\\\\\\\\n",
      "V^{n+1} & = V^n + k ( \\delta D_2 \\nabla^2_h V^n  + g(U^n,V^n)).\n",
      "\\end{align}$$\n",
      "In the cell below, write a function that advances the solution one step using this approach."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def one_step(u,v,k,A,delta,D1=0.5,D2=1.0):\n",
      "    u_new = u + k * (delta*D1*A*u + f(u,v))\n",
      "    v_new = v + k * (delta*D2*A*v + g(u,v))\n",
      "    \n",
      "    return u_new, v_new"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we'll use the function one_step() to solve the problem.  It turns out that to observe the behavior we're interested in, we can take an initial condition composed of random values in the interval $[-1/2, 1/2]$.  \n",
      "\n",
      "What time step should we use?  Since the diffusion terms are stiff, we may guess that it's sufficient to consider only them when selecting a time step size.  Recall that for Euler's method and centered differences, the 2D diffusion equation requires a step size\n",
      "$$k \\le \\frac{h^2}{4\\kappa}$$\n",
      "where $\\kappa$ is the diffusion coefficient.\n",
      "\n",
      "Here is the code:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "delta=0.0021; tau1=3.5; tau2=0; alpha=0.899; beta=-0.91; gamma=-alpha;\n",
      "\n",
      "def step_size(h,delta):\n",
      "    return h**2/(5.*delta)\n",
      "\n",
      "def pattern_formation(m=10,T=1000):\n",
      "    r\"\"\"Model pattern formation by solving a reaction-diffusion PDE on a periodic\n",
      "        square domain with an m x m grid.\"\"\"\n",
      "    D1 = 0.5\n",
      "    D2 = 1.0\n",
      "    \n",
      "    # Set up the grid\n",
      "    a=-1.; b=1.\n",
      "    h=(b-a)/m;                 # Grid spacing\n",
      "    x = np.linspace(a,b,m)     # Coordinates\n",
      "    y = np.linspace(a,b,m)\n",
      "    Y,X = np.meshgrid(y,x)\n",
      "\n",
      "    # Initial data\n",
      "    u=np.random.randn(m,m)/2.;\n",
      "    v=np.random.randn(m,m)/2.;\n",
      "\n",
      "    #plt.clf(); plt.hold(False)\n",
      "    #plt.pcolormesh(x,y,u); plt.colorbar(); plt.axis('image');\n",
      "    #plt.draw()\n",
      "\n",
      "    frames = [u]\n",
      "    \n",
      "    u=u.reshape(-1)\n",
      "    v=v.reshape(-1)\n",
      "\n",
      "    A=five_pt_laplacian_sparse_periodic(m,-1.,1.)\n",
      "\n",
      "    t=0.                      # Initial time\n",
      "    k = step_size(h,delta)    # Time step size\n",
      "    N = int(round(T/k))      # Number of steps to take\n",
      "\n",
      "    \n",
      "    #Now step forward in time\n",
      "    next_plot = 0\n",
      "    for j in range(N):\n",
      "\n",
      "        u,v = one_step(u,v,k,A,delta)\n",
      "        \n",
      "        t = t+k;\n",
      "\n",
      "        #Plot every t=5 units\n",
      "        if t>next_plot:\n",
      "            next_plot = next_plot + 5\n",
      "            U=u.reshape((m,m))\n",
      "            #plt.clf()\n",
      "            #plt.pcolormesh(x,y,U)\n",
      "            #plt.colorbar()\n",
      "            #plt.axis('image')\n",
      "            #plt.title(str(t))\n",
      "            #time.sleep(0.2)\n",
      "            #clear_output()\n",
      "            #fi=plt.gcf()\n",
      "            #display(fi)  \n",
      "            #plt.clf()\n",
      "            frames.append(U)\n",
      "            \n",
      "    return X,Y,frames"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from matplotlib import animation\n",
      "import matplotlib.pyplot as plt\n",
      "from clawpack.visclaw.JSAnimation import IPython_display\n",
      "import numpy as np\n",
      "\n",
      "fig = plt.figure(figsize=[4,4])\n",
      "\n",
      "U = frames[0]\n",
      "\n",
      "# This essentially does a pcolor plot, but it returns the appropriate object\n",
      "# for use in animation.  See http://matplotlib.org/examples/pylab_examples/pcolor_demo.html.\n",
      "# Note that it's necessary to transpose the data array because of the way imshow works.\n",
      "im = plt.imshow(U.T, vmin=U.min(), vmax=U.max(),\n",
      "           extent=[x.min(), x.max(), y.min(), y.max()],\n",
      "           interpolation='nearest', origin='lower')\n",
      "\n",
      "def fplot(frame_number):\n",
      "    U = frames[frame_number]\n",
      "    im.set_data(U.T)\n",
      "    return im,"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now let's use this code to solve the problem.  Try it first with small values of $m$, then for larger values."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x,y,frames = pattern_formation(m=200,T=300)\n",
      "\n",
      "animation.FuncAnimation(fig, fplot, frames=len(frames), interval=20)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "There is a problem with our approach, of course.  For larger values of $m$, the stiff diffusion term necessitates a very small time step, so the solution advances very slowly in time, so we would need to wait hours to get to the steady state.  Perhaps we can this problem by using a method that allows a larger step size."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Implicit-explicit splitting\n",
      "If we used the backward Euler method (or another A-stable method), we could use a much larger step size, but we would need to solve a large, nonlinear system of equations at each step.  To avoid this, we will split the equation into nonlinear parts and solve each part separately, i.e. first solve\n",
      "$$\\begin{align}\n",
      "u_t & = D_1 \\nabla^2 u  \\\\\\\\\n",
      "v_t & = D_2 \\nabla^2 v\n",
      "\\end{align}$$\n",
      "and then solve\n",
      "$$\\begin{align}\n",
      "u_t & = f(u,v)  \\\\\\\\\n",
      "v_t & = g(u,v).\n",
      "\\end{align}$$\n",
      "\n",
      "This is referred to as first-order operator splitting, or Godunov splitting.  It will reduce the temporal accuracy of our solution to first-order, but we are primarily interested in the final, steady-state solution, so that's acceptable.  We can use the backward Euler method for the diffusion terms -- to deal with the stiffness -- and the forward Euler method for the reaction terms:\n",
      "$$\\begin{align}\n",
      "U^{n+1/2} & = U^n + k \\delta D_1 \\nabla^2_h U^{n+1/2}  \\\\\\\\\n",
      "V^{n+1/2} & = V^n + k \\delta D_2 \\nabla^2_h V^{n+1/2} \\\\\\\\\n",
      "U^{n+1} & = U^{n+1/2} + k f(U^{n+1/2},V^{n+1/2})  \\\\\\\\\n",
      "V^{n+1} & = V^{n+1/2} + k g(U^{n+1/2},V^{n+1/2})\n",
      "\\end{align}$$\n",
      "\n",
      "In the cell below, write a function that advances the solution one step using this approach."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def one_step(u,v,k,A,delta,D1=0.5,D2=1.0):\n",
      "    # Your code here\n",
      "    return u_new, v_new"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def step_size(h,delta):\n",
      "    return h/(10.*delta)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pattern_formation(m=150,T=300)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "How fast is this code?  Why might that be?  What could we do to make it faster?  Can you implement something that is faster than either of these approaches, based on what you have learned in the course?"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# A spectral approach\n",
      "Solving the heat equation with periodic boundary conditions is trivial using a spectral method.  Starting with the heat equation\n",
      "\n",
      "$$u_t = \\kappa u_{xx}$$\n",
      "\n",
      "we take a Fourier transform and solve the resulting ODE to obtain\n",
      "\n",
      "$$\\hat{u}(\\xi,t) = e^{-\\kappa \\xi^2 t} \\hat{\\eta}(\\xi)$$\n",
      "\n",
      "where $\\hat{u}$ is the Fourier transform of $u$ and $\\hat{\\eta}$ is the Fourier transform of the initial data.  This suggests a very simple numerical method:\n",
      "\n",
      "1. Take the FFT of $U^n$.\n",
      "2. Multiply each value $\\hat{U}^n_\\xi$ by $\\exp(-\\kappa \\xi^2 \\Delta t)$.\n",
      "3. Take the inverse FFT of $\\hat{U}$.\n",
      "\n",
      "This method has no temporal discretization error, since we were able to solve the \"semi-discretization\" exactly.  Furthermore, it is stable for any time step.  Here's an implementation of that algorithm using numpy.  Notice that we only need to take the inverse FFT when we want output."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "epsilon = 0.05\n",
      "\n",
      "# Grid\n",
      "m = 64\n",
      "x = np.arange(-m/2,m/2)*(2*np.pi/m)\n",
      "k = 1./m**2\n",
      "tmax = 30.\n",
      "\n",
      "# Initial data\n",
      "u = np.sin(x)**2 * (x<0.)\n",
      "v = np.fft.fft(u)\n",
      "\n",
      "xi = np.array([range(0,m/2) + [0] + range(-m/2+1,0)])\n",
      "eps_xi2 = epsilon * xi**2.\n",
      "E = np.exp(-k * epsilon * xi**2.)\n",
      "\n",
      "nplt = np.floor((tmax/25)/k)\n",
      "nmax = int(round(tmax/k))\n",
      "\n",
      "for n in range(1,nmax+1):\n",
      "    v = E*v\n",
      "    t = n*k\n",
      "    \n",
      "    # Plotting\n",
      "    if np.mod(n,nplt) == 0:\n",
      "        u = np.squeeze(np.real(np.fft.ifft(v)))\n",
      "        plt.clf()\n",
      "        plt.plot(x,u,linewidth=3)\n",
      "        plt.title('t='+str(t))\n",
      "        plt.xlim((-np.pi,np.pi))\n",
      "        plt.ylim((0.,1.))\n",
      "        time.sleep(0.2)\n",
      "        clear_output()\n",
      "        fi=plt.gcf()\n",
      "        display(fi)  \n",
      "        plt.clf()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "You may say that this isn't very impressive, since we can solve this problem analytically.  You're right.  The power of this method is that we can use it to solve the diffusion part of a more complicated PDE, like the reaction-diffusion system above.  At each step, we just\n",
      "\n",
      "1. Solve \\begin{align} u_t & = \\delta D_1 \\nabla^2 u \\\\ v_t & = \\delta D_2 \\nabla^2 v \\end{align} by a Fourier spectral method.\n",
      "2. Solve \\begin{align} u_t & = f(u,v) \\\\ v_t & = g(u,v) \\end{align} using an initial value ODE solver.\n",
      "\n",
      "# Exercise 15\n",
      "Implement this method in the box below.  You may find it useful to first implement and test just your 2D Fourier spectral code.  You'll want to use the function `numpy.fft.fft2()`."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Try some of the following parameter values, or others that you choose.  What kind of patterns can you generate?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "delta=0.0045; tau1=2.02; tau2=0.; alpha=2.0; beta=-0.91; gamma=-alpha;\n",
      "delta=0.0005; tau1=2.02; tau2=0.; alpha=2.0; beta=-0.91; gamma=-alpha;\n",
      "delta=0.0021; tau1=3.5; tau2=0; alpha=0.899; beta=-0.91; gamma=-alpha;\n",
      "delta=0.0045; tau1=0.02; tau2=0.2; alpha=1.9; beta=-0.85; gamma=-alpha;\n",
      "delta=0.0001; tau1=0.02; tau2=0.2; alpha=0.899; beta=-0.91; gamma=-alpha;\n",
      "delta=0.0045; tau1=0.02; tau2=0.2; alpha=1.9; beta=-0.91; gamma=-alpha;"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}