This is a solution to the programming assignment for the fourth lesson in the numerical-mooc course: http://openedx.seas.gwu.edu/courses/GW/MAE6286/2014_fall
The reaction diffusion equation is given by: \begin{align} \frac{\partial u}{\partial t} &= D_u \nabla ^2 u - uv^2 + F(1-u)\\ \frac{\partial v}{\partial t} &= D_v \nabla ^2 v + uv^2 - (F + k)v \end{align} where $D_u$, $D_v$, $F$, and $k$ are constants, $u(x,y,t)$ and $v(x,y,t)$ represent the concentrations of two species of particle $U$ and $V$, and $\nabla^2$ is of course the Laplacian operator: $$\nabla^2 = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}.$$
We discretize $u$ and $v$ to be $u^n_{ij}$ and $v^n_{ij}$, where $i$ represents the $y$ coordinate and $j$ represents the $x$ coordinate, both starting from the lower-left corner. $n$ represents the index of time, starting from zero. We use forward time, central space. This means we get: \begin{align} \frac{u_{ij}^{n+1} - u_{ij}^n}{\Delta t} = \frac{D_u}{\delta^2}\left[\left(u^n_{i(j+1)} - 2 u^n_{ij} + u^n_{i(j-1)}\right) + \left(u^n_{(i+1)j} - 2 u^n_{ij} + u^n_{(i-1)j}\right)\right] - u^n_{ij} (v^n_{ij})^2 + F(1-u^n_{ij})\\ \frac{v_{ij}^{n+1} - v_{ij}^n}{\Delta t} = \frac{D_v}{\delta^2}\left[\left(v^n_{i(j+1)} - 2 v^n_{ij} + v^n_{i(j-1)}\right) + \left(v^n_{(i+1)j} - 2 v^n_{ij} + v^n_{(i-1)j}\right)\right] + u_{ij}^n (v_{ij}^n)^2 - (F+k)v_{ij}^n, \end{align} where $\delta = \Delta x = \Delta y.$
Or, \begin{align} u_{ij}^{n+1} = u_{ij}^n + \Delta t\left\{\frac{D_u}{\delta^2}\left[\left(u^n_{i(j+1)} - 2 u^n_{ij} + u^n_{i(j-1)}\right) + \left(u^n_{(i+1)j} - 2 u^n_{ij} + u^n_{(i-1)j}\right)\right] - u^n_{ij} (v^n_{ij})^2 + F(1-u^n_{ij})\right\}\\ v_{ij}^{n+1} = v_{ij}^n + \Delta t\left\{\frac{D_v}{\delta^2}\left[\left(v^n_{i(j+1)} - 2 v^n_{ij} + v^n_{i(j-1)}\right) + \left(v^n_{(i+1)j} - 2 v^n_{ij} + v^n_{(i-1)j}\right)\right] + u_{ij}^n (v_{ij}^n)^2 - (F+k)v_{ij}^n\right\}. \end{align} And we define $$\Delta t = \frac{9}{40} \frac{\delta^2}{\text{max}(D_u,D_v)}.$$
We impose von-Neumann boundary conditions and demand zero flux on all boundaries. This means that: \begin{align} u^n_{i0} = u^n_{i1}\ \forall\ i\\ v^n_{i0} = v^n_{i1}\ \forall\ i\\ u^n_{0j} = u^n_{1j}\ \forall\ j\\ v^n_{0j} = v^n_{1j}\ \forall\ j\\ u^n_{i(nx-1)} = u^n_{i(nx-2)}\ \forall\ i\\ v^n_{i(nx-1)} = v^n_{i(nx-2)}\ \forall\ i\\ u^n_{(ny-1)j} = u^n_{(ny-2)j}\ \forall\ j\\ v^n_{(ny-1)j} = v^n_{(ny-2)j}\ \forall\ j \end{align}
First we import our libraries
In [3]:
import numpy
from matplotlib import pyplot
import matplotlib.cm as cm
%matplotlib inline
And the initial data that we've been given.
In [4]:
uvinitial = numpy.load('./data/uvinitial.npz')
U = uvinitial['U']
V = uvinitial['V']
In [5]:
n = 192 # number of points in x and y.
Du, Dv, F, k = 0.00016, 0.00008, 0.035, 0.065 # Bacteria 1
dh = 5./(n-1)
T = 8000
dt = .9 * dh**2 / (4*max(Du,Dv))
nt = int(T/dt)
Here we define a number of helper functions that act on the 2d grid function u or v.
In [6]:
def impose_boundary_conditions(species):
"""
Imposes the boundary conidtions (which are von-Neumann)
on the 2d grid representing the value of the species
on the domain.
INPUT:
-- species, a 2d array of floats
representing either u or v
at a given time as a function
of space.
OUTPUT:
-- out, a 2d array of
floats representing species
but with boundary conditions
imposed.
"""
out = species.copy()
out[...,0] = out[...,1]
out[...,-1] = out[...,-2]
out[0,...] = out[1,...]
out[-1,...] = out[-2,...]
return out
In [7]:
def laplacian(species,delta):
"""
Calculates the laplacian of one of the
species either u or v.
Does not calculate the boundary points.
These are set by boundary conditions.
INPUT:
-- species, a 2d array of floats
representing either u or v
at a given time as a function
of space.
-- delta, the distance in x and y
between adjacent grid points.
OUTPUT:
-- out, a 2d array representing the
Laplacian of the species.
"""
s = species # convenient shorthand
out = numpy.zeros_like(s)
out[1:-1,1:-1] = ((1./delta)**2)\
* ((s[1:-1,2:] - 2*s[1:-1,1:-1] + s[1:-1,:-2])# d/dx
+ (s[2:,1:-1] - 2*s[1:-1,1:-1] + s[:-2,1:-1])) # d/dy
return out
In [8]:
def get_rhs(u,v,delta):
"""
Calculates du/dt for use in a
time integrator.
INPUTS:
--- u, a 2d array of floats
representing u as a function
of space
--- v, a 2d array of floats
representing v as a function
of space
--- delta, distance in space.
OUTPUTS:
--- u_rhs, the right-hand-side
for u.
--- v_rhs, the right-hand-side
for v.
"""
u_rhs = Du*laplacian(u,delta) - u*v*v + F*(1-u)
v_rhs = Dv*laplacian(v,delta) + u*v*v - (F+k)*v
return u_rhs,v_rhs
In [9]:
def euler_step(u,v,dt,dh):
"""
Integrates u and v forward by one time step.
INPUTS:
--- u, the U species as a function of
space. 2D array of floats.
--- v, the V species as a function of
space. 2D array of floats.
--- dt, float. The size of the time
step
--- dh, float. The distance between
adjacent grid points in space.
OUTPUTS:
--- u_new, 2D array of floats.
u after the time step.
--- v_new, 2D array of floats.
v after the time step.
"""
u_rhs,v_rhs = get_rhs(u,v,dh)
u_new = impose_boundary_conditions(u + dt*u_rhs)
v_new = impose_boundary_conditions(v + dt*v_rhs)
return u_new,v_new
And the integration method
In [10]:
def integrate(u_initial,v_initial,dt,dh,nt):
"""
Uses forward in time, centered inspace algorithm
to integrate u and v from their initial
configurations at time t = 0 to
their final configurations at time
t = nt*dt.
INPUTS:
--- u_initial, 2D array of floats.
The initial configuration of u.
--- v_initial, 2D array of floats.
The initial configuration of v.
--- dt, the duration of a single time step.
--- dh, the grid spacing in space
--- nt, the number of time steps.
OUTPUTS:
--- u, array of floats. The solution at
the final time for the U species
--- v, like u but for the V species.
--- t, the final time.
"""
# initialize arrays
u_old = u_initial.copy()
v_old = v_initial.copy()
t = 0
# And the loop
for k in range(nt):
u,v = euler_step(u_old,
v_old,
dt,dh)
u_old = u.copy()
v_old = v.copy()
t += dt
return u,v,t
In [11]:
u,v,t = integrate(U,V,dt,dh,nt)
In [12]:
t
Out[12]:
In [15]:
nt*dt
Out[15]:
In [16]:
u[100,::40]
Out[16]:
In [ ]: