In our previous article (<a href=»http://erg.biophys.msu.ru/wordpress/archives/706» alt=»discrete Laplacian»>about construction of the discrete Laplacian with numpy</a>) 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.