Problem set #6 PHY256¶
Reaction diffusion equations¶
The solutions of reaction–diffusion partial differential equations display a wide range of behaviours, including the formation of travelling waves and wave-like phenomena like solitons, and self-organized patterns like stripes, hexagons, or spots. Such patterns have been dubbed Turing patterns .
https://en.wikipedia.org/wiki/Turing_pattern
Reaction–diffusion systems and behavior described by reaction diffusion models are found in chemistry, biology, physics and other settings.
We consider time dependent concentrations of two chemicals $u,v$ on a periodic 2D domain $x,y$. The concentrations are $u(x,y,t)$, $v(x,y,t)$. The time dependent behavior is described with two partial differential equations \begin{align} \frac{\partial u}{\partial t} &= D_u \nabla^2 u + R_u(u,v) \\ \frac{\partial v}{\partial t} &= D_v \nabla^2 v + R_v(u,v) \end{align}
With functions $R_u()=0$ and $R_v()=0$, concentrations obey purely diffusive behavior with diffusion coefficients $D_u$, $D_v$. The diffusion rates of the two chemicals need not be the same.
The functions $R_u(u,v)$ and $R_v(u,v)$ represent chemical reactions between $u,v$ and growth and decay rates for $u,v$.
A mode that is stable in absence of diffusion may become unstable in presence of it. Instabilities induced by the presence of diffusion are called Turing Instabilities .
Brusellator Reaction Diffusion model¶
\begin{align} \frac{\partial u}{\partial t} &= D_u \nabla^2 u + \alpha - (\beta+1)u + u^2 v \\ \frac{\partial v}{\partial t} &= D_v \nabla^2 v + \beta u - u^2 v \end{align}Here $\alpha$ is a feeding rate for $u$. The parameter $\beta$ is a kill rate for $u$ that converts $u$ to $v$. The $uv^2$ term is a reaction term, producing $u$ at the expense of $v$.
Gray-Scott reaction diffusion model¶
\begin{align} \frac{\partial u}{\partial t} &= D_u \nabla^2 u - u v^2 + \alpha(1-u) \\ \frac{\partial v}{\partial t} &= D_v \nabla^2 v + u v^2 - (\alpha+\beta)v \end{align}$\alpha$ feeding rate for $u$ and a drain rate for $u,v$. $\beta$ gives a kill or drain rate for $v$. The $uv^2$ term is a reaction term, producing $v$ at the expense of $u$.
For these models $u,v > 0$.
Initial conditions:¶
The Brusselator model seems to give nice patterns with random, but low $u,v<0.05$ initial conditions.
The Gray-Scott model seems to make nice patterns with initial $u=1$, $v=0$ and some locations in the v array set to 1.
Parameters that give interesting behavior¶
Most parameters don't give interesting behavior. Here are some examples of parameters that give nice patterns with grid spacing $\Delta x=1$.
Brusselator giving nice spatial patterns: $$D_u=2, D_v=22, \alpha=5, \beta=9, dt = 0.005 $$
Brusselator giving both temporal and spatial patterns: $$D_u = 0.3, D_v = D_u/8, \alpha=2, \beta=5.4, dt=0.005$$
Gray-Scott giving nice spatial patterns: $$D_u = 0.1, D_v=D_u/2, \alpha = 0.037, \beta = 0.06, dt =1$$
Run for 10,000 steps.
Is it possible to predict when interesting behavior occurs? Yes. This is part of what Turing did with an an anlysis of (in)stability. He showed that small perturbations of certain wavelengths were likely to grow.
Numerical implementation¶
We approximate $$\frac {\partial u}{\partial t} \sim \frac{u^{n+1} - u^n}{\Delta t}$$ where $\Delta t$ is the timestep, and $u^n, u^{n+1}$ are $u$ values at consecutive times.
To update $u,v$ at timestep $n+1$ using information at timestep $n$ \begin{align} u^{n+1} &= u^n + \Delta t \left( D_u \nabla^2 u^n + R_u(u^n,v^n) \right) \\ v^{n+1} &= v^n + \Delta t \left( D_v \nabla^2 v^n + R_v(u^n,v^n) \right) \end{align}
We approximate a second derivative $$\frac{\partial^2 u}{\partial x^2} \sim \frac{u_{j-1} + u_{j+1} - 2 u_j}{\Delta x^2} $$ where $u_{j-1},u_j, u_{j+1}$ are $u$ values at consecutive grid locations in 1 dimension.
In 2 dimensions we give grid location with two indices $i,j$.
The Laplacian operator in 2-dimensions $$ \nabla = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}$$
The Laplacian at index $ij$ in the 2d grid can be computed as
$${\bf L1}: \qquad \nabla^2 u_{i,j} \approx \left( u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1} - 4 u_{i,j} \right) \frac{1}{ \Delta x^2} $$where $\Delta x$ is the grid spacing. This is a second order approximation to the Laplacian operator.
I have implemented this Laplacian in the example code below.
At time-step $n$ $${\bf L1}: \qquad \nabla^2 u_{i,j}^n \approx \left( u_{i+1,j}^n + u_{i-1,j}^n + u_{i,j+1}^n + u_{i,j-1}^n - 4 u_{i,j}^n \right) \frac{1}{ \Delta x^2} $$
The Laplacian operator is sometimes instead computed as
\begin{align} {\bf L2}: \qquad \nabla^2 u_{i,j}^n &\approx \left\{ 0.2 \left( u_{i+1,j}^n + u_{i-1,j}^n + u_{i,j+1}^n + u_{i,j-1}^n \right) - u_{i,j} \right. \\ & + 0.05 \left( u_{i+1,j+1}^n + u_{i-1,j-1}^n + u_{i-1,j+1}^n + u_{i+1,j-1}^n \right)\left. \right\} \frac{4}{ \Delta x^2} \end{align}Note these two approximations for the Laplacian do not give similar behavior unless the second one has that factor of 4. Below I ask you to write a routine to compute this second Laplacian operator!
We use periodic boundary conditions, which wraps the boundary in both x,y directions.
We can use the $\texttt{ numpy}$ routine $\texttt{ roll}$ to shift the arrays by one index vertically or horizontally or both, and with a wrap consistent with the periodic boundary conditions.
See https://docs.scipy.org/doc/numpy/reference/generated/numpy.roll.html
Above I have described a low order finite difference approximation to a reaction diffusion equation system in 2D on a grid. The recipe can also be described as a cellular automaton. See https://en.wikipedia.org/wiki/Cellular_automaton
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation, rc
# for animations
from IPython.display import HTML # requires ffmpeg
# you may need to install ffmpeg as like this on a command line: > conda install ffmpeg
%matplotlib inline
If you can't get the animations up and running, this assignment can be done without animations, by plotting many images.
# Laplacian operator evaluated in 2D
# arguments:
# dx is grid spacing, a is a 2D array
# roll shifts array by 1 index, in x or y directions
# return the computed Laplacian as 2d array
def laplacian2D(a, dx):
return (- 4*a + np.roll(a,1,axis=0) + np.roll(a,-1,axis=0)\
+ np.roll(a,+1,axis=1) + np.roll(a,-1,axis=1)) / (dx**2)
# it is faster to use a whole array routine than to loop over all indices
# returned is an array. You will write a new routine like this one for L2.
# Reaction parts of differential equations for
# different types of different Reaction diffusion equations
# arguments:
# u,v are 2D arrays,
# alpha, beta could be numbers or arrays
# these are parameters for the reaction diffusion equation
# rdtype is a string specifying the reaction diffusion model
# returns:
# 2d arrays that contain the computed non-diffusive parts of the
# reaction diffusion equation
def Reaction(u,v,alpha,beta,rdtype):
if (rdtype == 'Br'): # Brusellator model
Ru,Rv= R_Br(u,v,alpha,beta)
return Ru,Rv
if (rdtype == 'GS'): # Gray-Scott model
Ru,Rv = R_GS(u,v,alpha,beta)
return Ru,Rv
print('invalid reaction type')
# Brusselator Reaction Diffusion model
# compute the non-diffusive parts of the model
# example: Du = 2, Dv=22, alpha=5, beta=9, dt=0.005
def R_Br(u,v,alpha,beta):
Ru = alpha - (beta+1)*u + (u**2)*v
Rv = beta*u - (u**2)*v
return Ru,Rv
# The Gray-Scott Reaction Diffusion model
# alpha is feed rate for u, beta is kill or drain rate
# alpha kills v
# example: Du = 0.1, Dv=Du/2, alpha=0.037, beta=0.06, dt=1
def R_GS(u,v,alpha,beta):
# you will write a routine here!
Ru = 0*u# ? correct this
Rv = 0*u# ? correct this
return Ru, Rv
# a class for initializing and computing and showing solutions to
# a 2D reaction diffusion equation
# in this version I show both u,v images when plotted
# I also allow you to return arrays so if the animation doesn't work
# you can display images directly
class Update_RD(object):
# do this when creating a class object, axis info is passed,
# along with parameters for the reaction diffusion equation
# assumes you are passing figure and axis objects when initializing
def __init__(self, fig, ax, alpha=5.0,beta=9,Du=2,Dv=22,\
n=100,dt=0.01,dx=1,rdtype='Br',nsteps=1):
self.success = 0
self.alpha = alpha # parameter for reaction
self.beta = beta # parameter for reaction
self.Du = Du # diffusion coefficient for u
self.Dv = Dv # diffusion coefficient for v
self.n = n # dimension of 2D u,v arrays
self.rdtype = rdtype # type of reaction diffusion model, is a string
# the possible good string values are 'Br', 'GS'
self.dx = dx # horizontal grid spacing
self.dt = dt # timestep
self.nsteps = nsteps # how many timesteps to take per display
self.initialize_uv() #initialize arrays, fill with initial conditions
# axis info, two plots
self.line0, = ax[0].plot([], [], 'k-')
self.line1, = ax[1].plot([], [], 'k-')
self.fig = fig # passing this just in case we want to fuss with display
self.ax = ax
self.ax[0].set_aspect('equal') # so we get square plots
self.ax[1].set_aspect('equal')
self.ax[0].text(0,-n*0.05,'u',fontsize=15)
self.ax[1].text(0,-n*0.05,'v',fontsize=15)
# show the u,v images initial conditions in constructor
im0=self.ax[0].imshow(self.u,cmap='jet')
im1=self.ax[1].imshow(self.v,cmap='hsv')
#self.fig.colorbar(im0, ax=ax[0],fraction=0.046, pad=0.04)
#self.fig.colorbar(im1, ax=ax[1],fraction=0.046, pad=0.04)
# I wish I could get the colorbar to work!
fac_u= self.Du*self.dt/self.dx**2 # relevant for checking CFL-like condition
fac_v= self.Dv*self.dt/self.dx**2
if (fac_u>1):
print('cfl warning Du', fac_u)
if (fac_v>1):
print('cfl warning Dv', fac_v)
#### end constructor
# generate initial conditions for u,v arrays
def initialize_uv(self):
# I find that the Gray-Scott model needs these types of initial conditions
if (self.rdtype=='GS'):
self.u = np.zeros((self.n,self.n)) + 1.0 # u is 1
self.v = np.zeros((self.n,self.n)) # v is 0
for k in range(0,int(self.n/10)):
i = int(self.n*np.random.rand()) # randomly chose some locations to set
j = int(self.n*np.random.rand())
self.v[i,j]=np.random.rand() # set to number in [0,1)
else: # whereas the Brusellator seems to work with random
# but small initial conditions
self.u = 0.05*np.random.rand(self.n, self.n) # initial u data uniform distribution
self.v = 0.05*np.random.rand(self.n, self.n) # initial v data
# take a single dt timestep, update u,v arrays in place
def update(self):
du = self.Du*laplacian2D(self.u, self.dx) # compute diffusive terms
dv = self.Dv*laplacian2D(self.v, self.dx)
# compute reation terms
Ru,Rv = Reaction(self.u,self.v,self.alpha,self.beta,self.rdtype)
# update u,v arrays, Eulerian method
self.u = self.u + (du + Ru)*self.dt
self.v = self.v + (dv + Rv)*self.dt
# take nsteps, update u,v arrays in place
def update_nsteps(self):
for i in range(self.nsteps):
self.update(); # take a single step
# initialization for animation
def init(self):
lines = []
self.success = 0
self.line0.set_data([], [])
self.line1.set_data([], [])
return self.line0,self.line1,
# called each display step using animation,
# computation is only done when display is called
# this routine is designed to match our animation call!
def __call__(self, i):
if i == 0:
return self.init()
self.success += 1 #otherwise
# do nsteps times timestep for each display update
self.update_nsteps()
im0=self.ax[0].imshow(self.u,cmap='jet')
im1=self.ax[1].imshow(self.v,cmap='hsv')
#self.fig.colorbar(im0, ax=ax[0])
# it would be nice to have colorbars, but they screw up in the animation
return self.line0, self.line1,
# return u,v arrays, whatever their current values happen to be
# (for testing and used if animations don't work)
def return_uv(self):
return self.u, self.v
Problem 1. Stability¶
If $\Delta t$ is too large, the solution is unstable or extremely sensitive to $\Delta t$. Below what level for $\Delta t$ gives you consistent results for the Brusellator model with $\alpha=5, \beta=9$, $D_u = 2$, $D_v=22$? Show that you don't get similar or consistent results for similar $\Delta t$ if $\Delta t$ is too high. The numerical method becomes unstable. What condition gives stable simulations?
Problem 2. Scaling in the grid spacing¶
Show that changing grid spacing $\Delta x$ is equivalent to changing diffusion coefficients $D_u$ and $D_v$ for these reaction diffusion models. This can be done analytically.
Problem 3. Write the reaction function for the Gray-Scott model¶
Write the routine R_GS(u,v,alpha,beta) giving the reaction rates for the Gray-Scott model.
Problem 4. Writing a Laplacian operator¶
Use the $\texttt{numpy}$ routine $\texttt{roll()}$ to write a Laplacian $\nabla^2$ function that uses neighboring positions including the diagonal positions, as shown above in the markdown section where it is labelled L2. See if changing the Laplacian routine affects the morphology of patterns generated. I did not see much change in the behavior.
Problem 5. Varying parameters $\alpha, \beta$ across the grid¶
What happens if $\alpha$ and $\beta$ are functions of position in the grid? Modify/rewrite the routines to integrate a reaction diffusion equation with slowly varying $\alpha$, $\beta$ across the grid. See if you can show how the formed patterns vary as a function of $\alpha,\beta$.
I found that abrupt changes in $\alpha$ or $\beta$ gave linear features, but a periodic or sinusoidal variation in $\alpha$ and $\beta$ could make an image with very nice patterns slowly varying across it. An example is shown at the way bottom below.
Optional! You could also explore variations in the diffusion coefficients $D_u, D_v$.
# set up the Brusselator model with some nice parameters giving
# pattern formation
fig, ax = plt.subplots(1,2,figsize=(8,5)) # set up a figure with two panels
ud_Br = Update_RD(fig,ax, alpha=5.0,beta=9,Du=2,Dv=22,\
n=100,dt=0.005,dx=1,rdtype='Br',nsteps=200)
# ud_Br is the reaction diffusion object
# create this object, assign parameters
# here n is the grid length, dt is the timestep, dx is the horizontal grid spacing
# alpha,beta, Du, Dv are the parameters for the reaction diffusion model
# the string 'Br' sets model to be the Brusselator model
# the arrays have size 100x100 (sete by n)
# nsteps=200 is the number of steps between display outputs
anim_Br = animation.FuncAnimation(fig, ud_Br, frames=np.arange(50), init_func=ud_Br.init,
interval=100, blit=True)
# here 50 is the number of images shown in the animation
# this last step sets up the animation, connecting output of the integration
# to the animation.
# what is shown below is the the initialization of the array
# how finally we compute and show the animation
#HTML(anim_Br.to_html5_video())
HTML(anim_Br.to_jshtml()) # I like this type of display even better!
# image on the left is u, that on the right is v