In our previous article about construction of the discrete Laplacian with numpy we discuss the discrete Laplacian calculation for numerical solving of diffusion PDE. Here we generalize this technique for the  two-dimensional space. We describe the complete algorithm of solving of reaction-diffusion equations with implementation of Crank-Nicolson scheme for diffusion part together with Runge-Kutta scheme for reactional part.

Two-dimensional discrete Laplacian

2D Laplacian operator can be described with matrix N2xN2, where N is a grid spacing of a square reactor. Note that we consider a square reactor, but the same technique can be easily modified for a rectangular reactor MxN.

Consider the most rough discretization 4×4 of the reactor with periodic boundary conditions. Let us fix the point (2,2) «o» and denote with «i» all points that contribute to the second space derivative in the point «o»:

   x i x x
   i o i x
   x i x x
   x x x x

One may say that point (1,1) also contributes to second derivative in (2,2) point, but here we consider only the contribution of  left/right and top/bottom. This approximation is called a «5-point stensil»:

\frac{\partial^2 f}{\partial x^2} = \frac{4f(x,y)-f(x-h,y)-f(x+h,y)-f(x,y-h)-f(x,y+h)}{h^2} (1)

The discrete Laplacian for the 4×4 reactor has the form:

-4  1  0  1  0  0  0  0  0  0  0  0  1  0  0  0 
 1 -4  1  0  1  0  0  0  0  0  0  0  0  1  0  0 
 0  1 -4  1  0  1  0  0  0  0  0  0  0  0  1  0 
 1  0  1 -4  1  0  1  0  0  0  0  0  0  0  0  1 
 0  1  0  1 -4  1  0  1  0  0  0  0  0  0  0  0 
 0  0  1  0  1 -4  1  0  1  0  0  0  0  0  0  0 
 0  0  0  1  0  1 -4  1  0  1  0  0  0  0  0  0 
 0  0  0  0  1  0  1 -4  0  0  1  0  0  0  0  0 
 0  0  0  0  0  1  0  1 -4  1  0  1  0  0  0  0 
 0  0  0  0  0  0  1  0  1 -4  1  0  1  0  0  0 
 0  0  0  0  0  0  0  1  0  1 -4  1  0  1  0  0 
 0  0  0  0  0  0  0  0  1  0  1 -4  1  0  1  0 
 1  0  0  0  0  0  0  0  0  1  0  1 -4  1  0  1 
 0  1  0  0  0  0  0  0  0  0  1  0  1 -4  1  0 
 0  0  1  0  0  0  0  0  0  0  0  1  0  1 -4  1 
 0  0  0  1  0  0  0  0  0  0  0  0  1  0  1 -4

It can be easily verified that this matrix calculates the second derivative according to the equation (1). «+1» and «-1» diagonals correspond to the «left point» and «right point» contributions to the second derivative. «+4» and «-4» diagonals correspond to the  «upper point» and «bottom point» contributions. «+12» and «-12» diagonals give the periodicity. Let us generalize this result for any N using the following function:

def gen_laplacian_pbc(N):
   data1 = [-4]*N**2
   data2 = [1]*N**2
   L = sp.spdiags( [ data2, data2, data2, data1, data2, data2, data2 ], \
     [ -N*(N-1), -N, -1, 0, 1, +N, +N*(N-1) ], \
     N**2, N**2 )
   return L.tocsc()

Diffusion equation in 2D space

The algorithm developed for the 1D space can be slightly modified for 2D calculations. To apply the Laplacian we should linearize the matrix of function values:

v_lin = v.ravel()

For visualization, this linearized vector should be transformed to the initial state:

v = v_lin.reshape((num,num))
imshow(v)

Finally, we can write the algorithm solving the Fick equation:

u = initcond(num) # some function defining initial conditions
L = gen_laplacian_pbc(num)
dt = 0.005
v = np.copy(u).ravel()
i = 0

print 'Lets draw!'
figure()
imshow(u)
show(block=False)

while not Stop:
   i += 1
   v += dt*D*L.dot(v)

   if i % 1000 == 0:
      print '.',
      sys.stdout.flush()
      imshow(v.reshape((num,num)))
      draw()

show(block=True)

Implicit scheme for solving the diffusion equation

To render the calculations more stable we can use the implicit scheme. We calculate the space derivative using simultaneously vector values from current and next steps:

  \frac{V_{2} - V_{1}}{\tau} = \frac{D}{2}(\Delta V_{1} + \Delta V_{2})

That is equivalent to:

  \left(\frac{1}{\tau}\mathbb{E}-\frac{D}{2}\Delta\right)V_{2}=\frac{V_{1}}{\tau}+\frac{D}{2}\Delta V_{1}

Thus we obtained a standard problem of linear equation calculation. The matrix Δ is sparse, and thus, matrix \left(1/\tau\mathbb{E}-D\Delta/2\right) is also sparse. Moreover, for 1D case this matrix is band tridiagonal matrix and the numerical algorithm has complexity of O(N) where N is a number of cells (space discretization).

By L± we denote matrix \left(1/\tau\mathbb{E}\pm D\Delta/2\right); then we can rewrite the problem as:

L_{-}\cdot v=L_{+}\cdot v (2)

The fastest way to solve this problem is to (i) build preliminary matrices; (ii) apply the L+ part; (iii) solve the linear system with L- matrix:

ECSC = sp.identity(num, format='csc')
L = gen_laplacian_pbc(num)
LM = ECSC/dt - D/2*L
LP = ECSC/dt + D/2*L

As we have to solve the same linear problem many-many times we can optimize the algorithm by using the preliminary LU decomposition implemented in method splu() of sparse matrix object. This method returns an object that has method solve():

LMlu = lg.splu(LMo)
v_s = LMlu.solve(v)

Finally, the code for PDE solution using implicit scheme remains very simple:

u = initcond(num)
v = np.copy(u)
while True:
i += 1
v = LM.solve(LP.dot(v))

Implicit scheme for solving a reaction-diffusion equation

Let us demonstrate how the discrete Laplacian is used in explicit scheme of solving the reaction-diffusion system

\frac{\partial\vec{v}}{\partial t}=R(\vec{v})+D\Delta\vec{v}

First several steps of the explicit algorithm are the following:

  1. Let V_0 be starting conditions.
  2. Use the R mapping for finding V_{R1}=V_0+\tau R(V_0).
  3. Applying the discrete Lapplacian for finding the diffusion part: V_{D1}=V_{R1}+\tau D\Delta V_{R1}.
  4. Repeat step 2 computing V_{R2} from V_{D1}.
  5. Repeat step 3 computing V_{D2} from V_{R2}.
  6. Etcetera..

This is the starting point for a numerical scheme for improving both reaction and diffusion part. For improving of the reactional part  we calculate the change at step 2 with Runge-Kutta technique.