```
%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
```

    Populating the interactive namespace from numpy and matplotlib


# 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.


```
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.


```
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."""
    e=np.ones(m**2)
    e2=([1]*(m-1)+[0])*m
    e3=([0]+[1]*(m-1))*m
    h=(b-a)/(m+1)
    A=scipy.sparse.spdiags([-4*e,e2,e3,e,e],[0,-1,1,-m,m],m**2,m**2)
    # 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/=h**2
    A = A.todia()
    return A
```


```
A = five_pt_laplacian_sparse_periodic(4,-1.,1.)
plt.spy(A)
```

# 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.


```
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:


```
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
    u=np.random.randn(m,m)/2.;
    v=np.random.randn(m,m)/2.;

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

    u=u.reshape(-1)
    v=v.reshape(-1)

    A=five_pt_laplacian_sparse_periodic(m,-1.,1.)

    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
            U=u.reshape((m,m))
            plt.clf()
            plt.pcolormesh(x,y,U)
            plt.colorbar()
            plt.axis('image')
            plt.title(str(t))
            time.sleep(0.2)
            clear_output()
            fi=plt.gcf()
            display(fi)  
            plt.clf()
```

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


```
pattern_formation(m=150,T=300)
```

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.


```
def one_step(u,v,k,A,delta,D1=0.5,D2=1.0):
    # Your code here
    return u_new, v_new
```


```
def step_size(h,delta):
    return h/(10.*delta)
```


```
pattern_formation(m=150,T=300)
```

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?

Using either code, try some of the following parameter values.  What kind of patterns can you generate?


```
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;
```


```
print scipy.sparse.eye(4,4)
```


```
scipy.sparse.linalg.spsolve
```

837cea07b27cf5eebc36093a83d2358e5c072229
ipynb
Jupyter Notebook
ipython/A reaction-diffusion model for pattern formation.ipynb
bfdeandrade/Finite-Difference
97e93f394051a70b6aa2c26fa266952ae821deb2
CC-BY-2.0
2015-02-05T23:13:07.000Z
2022-01-05T03:09:04.000Z
ipython/A reaction-diffusion model for pattern formation.ipynb
MuriloHMoreira/finite-difference-course
97e93f394051a70b6aa2c26fa266952ae821deb2
CC-BY-2.0
ipython/A reaction-diffusion model for pattern formation.ipynb
MuriloHMoreira/finite-difference-course
97e93f394051a70b6aa2c26fa266952ae821deb2
CC-BY-2.0
2015-02-16T17:36:38.000Z
2021-12-11T07:13:03.000Z
Qwen/Qwen-72B
1. YES
2. YES