## Explicit Methods for Reaction AdvectionDiffusion Equations

The equation whose solution we will focus on is u — D(u, x,y)q2xu + g(u, x, t)dxu + f (u, x, t), u(x, to) — «o(x) 

where t and x are the time and space variables and u = 8t«. (A generalization of  to two spatial dimensions will be considered later.) The three terms on the RHS of  are associated with diffusion, advection, and reaction, respectively. To illustrate the main issues ofnumerical solution of , we first consider its simpler version - the heat equation - where we retain only the diffusion term and take the diffusion coefficient D to be a constant:

The last two equations are called boundary conditions; such (or similar) conditions must always be supplied for a partial differential equation (PDE). Boundary conditions that fix the value of the unknown, as in , are called Dirichlet boundary conditions. Boundary condition of other types

 ( n +1 r A r J n m ) V -1 m J V m ) + 1

Figure 1 The stencil of method .

will be considered later. Note that setting the boundary conditions  to zero does not restrict the generality of our treatment. Indeed, for the more general conditions, «(0) = a(t), «(X) = b(t), one defines a new variable w = « — [(X — x)a(t) + xb(t)]/X, which satisfies  and .

We discretize tn = t0 + nr and xm = 0 + mh, where r and h are the time step and the mesh size, and n = 0, 1, 2,... and m = 0, 1, ..., X/h=M. Denote «nmx «(xm, tn) and Um to be, respectively, the exact and the numerical solutions at node (xm, tn). Replacing

one obtains the finite-difference approximation to :

Un+1 = rUl+i + (1 - 2r )un + rUl -1, un = 0, UM = 0 

where

The meaning of the notation O(rp) is explained in Numerical Methods for Local Models. Equation  is the counterpart of the explicit Euler method for ordinary differential equations (ODEs). As follows from , its accuracy is 0(t) + O(h2). (However, for an initial condition uo(x) that is either discontinuous or has a discontinuous slope, the spatial accuracy is reduced: the smaller the r in , the closer the accuracy to the maximum possible value O(h ).) The numerical solution at node (xm, tn + j) can thus be found from  and the initial condition lfm — u0(xm) if one knows the solution at nodes

±!, tn). These four nodes form a so-called stencil for scheme , as shown schematically in Figure 1. 