In [ ]:
%pylab inline
import numpy as np
import scipy.optimize
import scipy.sparse
import matplotlib.pyplot as plt
from IPython.display import clear_output
import time

The model

In this notebook we'll solve a system of reaction-diffusion PDEs in two dimensions:

\begin{align} u_t & = \delta D_1 \nabla^2 u + f(u,v) \\\\ v_t & = \delta D_2 \nabla^2 v + g(u,v) \end{align}

where $\nabla^2 u = u_{xx} + u_{yy}$ denotes the Laplacian and $f,g$ represent reaction terms.

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 \begin{align} u(-1,y,t) & = u(1,y,t) \\\\ u(x,-1,t) & = u(x,1,t) \end{align} with corresponding conditions on $v$.

Reaction equations

The reaction terms we will use are \begin{align} f(u,v) & = \alpha u (1-\tau_1 v^2) + v(1-\tau_2 u) \\\\ g(u,v) & = \beta v + \alpha \tau_1 u v^2 + u (\gamma + \tau_2 v). \end{align} 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.

In [ ]:
delta=0.0045; tau1=0.02; tau2=0.2; alpha=0.899; beta=-0.91; gamma=-alpha;

def f(u,v):
    return alpha*u*(1-tau1*v**2) + v*(1-tau2*u);

def g(u,v):
    return beta*v*(1+alpha*tau1/beta*u*v) + u*(gamma+tau2*v);

Diffusion terms

Next let's write a discretization of the diffusion terms. We'll use the 5-point Laplacian approximation:

\begin{align} \nabla^2 u(x_i,y_j) & \approx \frac{U_{i-1,j} - 2U_{ij} + U_{i+1,j}}{\Delta x} + \frac{U_{i,j-1} - 2U_{ij} + U_{i,j+1}}{\Delta y} \end{align}

For simplicity we'll assume $\Delta x = \Delta y = h$. Examine the code below carefully, noting how the periodic boundary conditions are imposed.

In [ ]:
def five_pt_laplacian_sparse_periodic(m,a,b):
    """Construct a sparse matrix that applies the 5-point laplacian discretization
       with periodic BCs on all sides."""
    # Top & bottom BCs:
    A_periodic = scipy.sparse.spdiags([e,e],[m-m**2,m**2-m],m**2,m**2).tolil()
    # Left & right BCs:
    for i in range(m):
        A_periodic[i*m,(i+1)*m-1] = 1.
        A_periodic[(i+1)*m-1,i*m] = 1.
    A = A + A_periodic
    A = A.todia()
    return A
In [ ]:
A = five_pt_laplacian_sparse_periodic(4,-1.,1.)

Forward Euler time stepping

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: $$\begin{align} U^{n+1} & = U^n + k ( \delta D_1 \nabla^2_h U^n + f(U^n,V^n)) \\\\ V^{n+1} & = V^n + k ( \delta D_2 \nabla^2_h V^n + g(U^n,V^n)). \end{align}$$ In the cell below, write a function that advances the solution one step using this approach.

In [ ]:
def one_step(u,v,k,A,delta,D1=0.5,D2=1.0):
    u_new = u + k * (delta*D1*A*u + f(u,v))
    v_new = v + k * (delta*D2*A*v + g(u,v))
    return u_new, v_new

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]$.

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 $$k \le \frac{h^2}{4\kappa}$$ where $\kappa$ is the diffusion coefficient.

Here is the code:

In [ ]:
delta=0.0021; tau1=3.5; tau2=0; alpha=0.899; beta=-0.91; gamma=-alpha;

def step_size(h,delta):
    return h**2/(5.*delta)

def pattern_formation(m=10,T=1000):
    r"""Model pattern formation by solving a reaction-diffusion PDE on a periodic
        square domain with an m x m grid."""
    D1 = 0.5
    D2 = 1.0
    # Set up the grid
    a=-1.; b=1.
    h=(b-a)/m;                 # Grid spacing
    x = np.linspace(a,b,m)     # Coordinates
    y = np.linspace(a,b,m)
    Y,X = np.meshgrid(y,x)

    # Initial data

    #plt.clf(); plt.hold(False)
    #plt.pcolormesh(x,y,u); plt.colorbar(); plt.axis('image');

    frames = [u]


    t=0.                      # Initial time
    k = step_size(h,delta)    # Time step size
    N = int(round(T/k))      # Number of steps to take

    #Now step forward in time
    next_plot = 0
    for j in range(N):

        u,v = one_step(u,v,k,A,delta)
        t = t+k;

        #Plot every t=5 units
        if t>next_plot:
            next_plot = next_plot + 5
    return X,Y,frames
In [ ]:
from matplotlib import animation
import matplotlib.pyplot as plt
from clawpack.visclaw.JSAnimation import IPython_display
import numpy as np

fig = plt.figure(figsize=[4,4])

U = frames[0]

# This essentially does a pcolor plot, but it returns the appropriate object
# for use in animation.  See
# Note that it's necessary to transpose the data array because of the way imshow works.
im = plt.imshow(U.T, vmin=U.min(), vmax=U.max(),
           extent=[x.min(), x.max(), y.min(), y.max()],
           interpolation='nearest', origin='lower')

def fplot(frame_number):
    U = frames[frame_number]
    return im,

Now let's use this code to solve the problem. Try it first with small values of $m$, then for larger values.

In [ ]:
x,y,frames = pattern_formation(m=200,T=300)

animation.FuncAnimation(fig, fplot, frames=len(frames), interval=20)

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.

Implicit-explicit splitting

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 $$\begin{align} u_t & = D_1 \nabla^2 u \\\\ v_t & = D_2 \nabla^2 v \end{align}$$ and then solve $$\begin{align} u_t & = f(u,v) \\\\ v_t & = g(u,v). \end{align}$$

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: $$\begin{align} U^{n+1/2} & = U^n + k \delta D_1 \nabla^2_h U^{n+1/2} \\\\ V^{n+1/2} & = V^n + k \delta D_2 \nabla^2_h V^{n+1/2} \\\\ U^{n+1} & = U^{n+1/2} + k f(U^{n+1/2},V^{n+1/2}) \\\\ V^{n+1} & = V^{n+1/2} + k g(U^{n+1/2},V^{n+1/2}) \end{align}$$

In the cell below, write a function that advances the solution one step using this approach.

In [ ]:
def one_step(u,v,k,A,delta,D1=0.5,D2=1.0):
    # Your code here
    return u_new, v_new
In [ ]:
def step_size(h,delta):
    return h/(10.*delta)
In [ ]:

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?

A spectral approach

Solving the heat equation with periodic boundary conditions is trivial using a spectral method. Starting with the heat equation

$$u_t = \kappa u_{xx}$$

we take a Fourier transform and solve the resulting ODE to obtain

$$\hat{u}(\xi,t) = e^{-\kappa \xi^2 t} \hat{\eta}(\xi)$$

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:

  1. Take the FFT of $U^n$.
  2. Multiply each value $\hat{U}^n_\xi$ by $\exp(-\kappa \xi^2 \Delta t)$.
  3. Take the inverse FFT of $\hat{U}$.

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.

In [ ]:
epsilon = 0.05

# Grid
m = 64
x = np.arange(-m/2,m/2)*(2*np.pi/m)
k = 1./m**2
tmax = 30.

# Initial data
u = np.sin(x)**2 * (x<0.)
v = np.fft.fft(u)

xi = np.array([range(0,m/2) + [0] + range(-m/2+1,0)])
eps_xi2 = epsilon * xi**2.
E = np.exp(-k * epsilon * xi**2.)

nplt = np.floor((tmax/25)/k)
nmax = int(round(tmax/k))

for n in range(1,nmax+1):
    v = E*v
    t = n*k
    # Plotting
    if np.mod(n,nplt) == 0:
        u = np.squeeze(np.real(np.fft.ifft(v)))

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

  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.
  2. Solve \begin{align} u_t & = f(u,v) \\ v_t & = g(u,v) \end{align} using an initial value ODE solver.

Exercise 15

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().

In [ ]:

Try some of the following parameter values, or others that you choose. What kind of patterns can you generate?

In [ ]:
delta=0.0045; tau1=2.02; tau2=0.; alpha=2.0; beta=-0.91; gamma=-alpha;
delta=0.0005; tau1=2.02; tau2=0.; alpha=2.0; beta=-0.91; gamma=-alpha;
delta=0.0021; tau1=3.5; tau2=0; alpha=0.899; beta=-0.91; gamma=-alpha;
delta=0.0045; tau1=0.02; tau2=0.2; alpha=1.9; beta=-0.85; gamma=-alpha;
delta=0.0001; tau1=0.02; tau2=0.2; alpha=0.899; beta=-0.91; gamma=-alpha;
delta=0.0045; tau1=0.02; tau2=0.2; alpha=1.9; beta=-0.91; gamma=-alpha;