Speed of reaction-diffusion fronts in spatially heterogeneous media

Departament de Medicina, Facultat de Cie `ncies de la Salut, Universitat Internacional de Catalunya, c/ Gomera s/n, 08190-Sant Cugat del Valle `s (Barcelona), Spain Departament de Fı ́sica, Universitat de Girona, Campus Montilivi, 17071 Girona, Catalonia, Spain Department of Mathematics and Center for Biodynamics, Boston University, Boston, Massachusetts 02215, USA Department of Mathematics, UMIST-University of Manchester, Institute of Science and Technology, Manchester M60 1DQ, United Kingdom ~Received 16 September 2002; revised manuscript received 15 May 2003; published 13 October 2003 !


I. INTRODUCTION
Front propagation modeled by reaction-diffusion equations has been applied in many areas of science such as physics, biology, ecology, and chemistry ͓1,2͔.Since the pioneering works by Kolmogorov, Petrovskii, and Piskunov ͑KPP͒ ͓3͔ and Fisher ͓4͔, both in 1937, this field has been continuously growing.The basic phenomena have been described by using parabolic reaction-diffusion equations derived under the assumption that the medium in which fronts are moving is homogeneous.Although heterogeneities are always present in nature, studies of fronts in heterogeneous media have been much more recent.Some examples are porous media, random media, noisy media, ecological patchiness, etc.
Experimental studies have been developed for heterogeneous excitable media.In this context, two-dimensional front propagation in the photosensitive Belousov-Zhabotinsky modulated reaction ͓5,6͔ and patchy media ͓7͔ have been explored.Successful theoretical efforts have been also made to understand the phenomenon of front propagation in excitable media.In particular, Xin ͓8͔ has studied front solutions for reaction-diffusion equations in periodic and random media ͓8͔, and Shigesada et al. ͓9͔ have given analytical restrictions for the existence of propagating fronts in ecological patchy environments.Heterogeneous models have been also used, via computer simulations, to describe the dynamics of brain tumors ͓10͔.Moreover, diffusion coefficients depending on spatial and temporal coordinates have been recently proposed to study the formation of Alzheimer's disease senile plaques ͓11͔.Nakamura et al. ͓12͔, have studied reaction-diffusion equations when the spatially inhomogeneous reaction rate is much larger than the diffusion coefficient.Keener ͓13,14͔, Mitkov et al. ͓15͔ and Rotstein et al. ͓16͔ have studied bistable-type models in heterogeneous media in the context of calcium release waves ͓17͔.Petrovskii has analyzed the case of a spatially periodic environment ͓18͔.Many authors have recently studied the effect of the external multiplicative noise on the speed and the width of the front.These studies may be seen as a way to introduce stochastic heterogeneities in the reaction-diffusion equation through a control parameter ͓19͔.Therefore, the effect of heterogeneities on front propagation is of wide theoretical and practical interest.
Our goal in this work is to understand how deterministic heterogeneities influence the front speed of parabolic reaction-diffusion equations with a monostable reaction term, when either the diffusion coefficient or the reaction term depend on the spatial coordinate.Methods such as marginal ͑linear͒ stability ͓20͔ and variational ͓21͔ analysis have been widely used to find the asymptotic speed of a front.However, both methods do not hold, or at least they should be adequately generalized, when the reaction-diffusion equation has an explicit dependence on the spatial coordinate.Instead, we will make use of well-known analytical techniques such as singular perturbation analysis and the local speed approach, both valid for weak heterogeneities, and geometrical methods for general heterogeneities, in order to study how heterogeneities introduce corrections to the asymptotic front speed, both for pulled ͑KPP͒ and pushed ͑but monostable͒ fronts.We will also compare the analytical results and numerical simulations.
The methods we employ here have some limitations.Singular perturbative analysis may be efficiently compared to numerical results when the solution to the leading order is known and for reaction-diffusion equations with non-KPP kinetic terms.The solution for the lowest order may be found for some particular non-KPP kinetic terms but it is not known in general, although in those cases numerical solution may always be calculated.This method requires, of course, a small parameter present in the model.Therefore, it is necessary to assume that the spatial heterogeneities of the media introduce a small variation in the reaction rate and/or the diffusion coefficient ͑weak heterogeneities͒ and the characteristic length of the heterogeneities must be greater than the characteristic width of the front ͑smooth heterogeneities͒ ͓22͔.On the other hand, the geometrical method we present here, based on Hamilton-Jacobi dynamics, only holds for KPP kinetics but, contrary to the previous method, there is no need to assume either weak or smooth heterogeneities.The local speed approach is based in the assumption that for weak and smooth heterogeneities the speed of the front is given by the local value of the reaction rate U and/or the diffusion coefficient D in each spatial point, i.e., the front speed would be vӍ2ͱU(x)D(x) for KPP kinetics.Therefore, fronts with Uϭb ͑constant͒ and Dϭ f (x) would have vӍ2ͱb f (x), and the same would hold for fronts in media with Uϭ f (x) and Dϭb.But, as we shall see, in general this simple approach is not consistent either with our analytical results or with our numerical simulations.
In this paper we study first the dynamics of front motion for the following smooth heterogeneous problems where the function f satisfies f (0)ϭ f (1)ϭ0, (x,0) ϭ(x) where (x) is an initial condition that may range from the Heaviside step function ͓(x,0)ϭ1 for xϽ0 and (x,0)ϭ0 for xϾ0] to a fully developed front, D and U are the dimensionless diffusion coefficient and reaction rate, respectively, and is a small parameter.Since we expect solutions to behave like totally developed fronts, we should look at them in the asymptotic regime ͑large-space and largetime limit͒ t→t/ and x→x/.The scaling considered is equivalent to assuming that the front is totally developed independent of the way it developed from initial conditions.Equations ͑1͒ then become Consistent with the initial conditions and the existence of a front we require the solution to satisfy lim x→Ϫϱ ϭ1 and lim x→ϱ ϭ0.

II. NONUNIFORM REACTION
We consider the problem where ␦ is the amplitude of the heterogeneities and u (x) is the reactive heterogeneity.

A. Singular perturbation analysis
This method of perturbative analysis has been already employed to study the speed of pulled fronts and it has been shown that the solvability integrals diverge ͓23͔.Therefore, we will use this method only for non-KPP kinetics.We assume ␦ϭO() ͑weak heterogeneities͒, i.e., ␦ϵ, where ϭO(1) in Eq. ͑3͒.Equation ͑3͒, together with the corresponding boundary conditions, becomes In order to study Eqs.͑4͒ we will make a nonrigorous asymptotic analysis.We assume that the domain is divided into two regions according to the space scales: a boundary layer region, whose width is O(), in which is rapidly varying; and an external region in which is almost constant, i.e., either ϭO( n 1 ) or ϭ1ϩO( n 2 ), where n 1 and n 2 are positive real numbers.
In order to solve Eq. ͑4͒ in the outer region we expand as follows: By substituting Eq. ͑5͒ into Eq.͑4͒ and collecting terms with the same powers of we get The solution of Eq. ͑6͒ is ⌽ 0 ϭ1 to the left of the boundary layer and ⌽ 0 ϭ0 to the right of the boundary layer.The solutions of Eqs.͑7͒ and ͑8͒ are ⌽ 1 ϵ0 and ⌽ 2 ϵ0, respectively.Thus, (x,t;)ϭO( 3 ) to the left of the boundary layer and (x,t;)ϭ1ϩO( 3 ) to the right of the boundary layer.Note that to the order of magnitude considered here, there is no effective difference between the homogeneous and heterogeneous cases, i.e., there is no difference in the shape of the front.
In order to study the dynamics in the interior of the boundary layer we translate Eq. ͑4͒ to the reference frame of the front, i.e., we define the new variable zϭ͓xϪS(t)͔/ where S(t) represents the position of the front.The derivatives in Eqs.͑2͒ transform according to where the dot symbol stands for the temporal derivative.We expand and S in powers of : and, in consequence, where u Ј(S 0 )ϭd u (x)/dx͉ xϭS 0 and f Ј( 0 )ϭd f ()/ d͉ ϭ 0 .Inserting Eqs.͑10͒ and ͑11͒ into Eq.͑4͒ once Eqs.
͑9͒ are taken into account, and collecting terms with equal powers of one gets O(1), O(), and O( 2 ), respectively:

͑14͒
where Lϭ‫ץ‬ zz ϩS ‫ץ0˙‬ z ϩ f ( 0 ) and L 1 ϭ‫ץ‬ zz ϩS ‫ץ0˙‬ z ϩ f Ј( 0 ).Since we assume 0 ϭ 0 (z), Eq. ͑12͒ is equivalent to the homogeneous (ϭ0) parabolic reaction-diffusion equation translated to the front reference frame (zϭxϪS ˙0t) which travels at constant speed S ˙0.We call S ˙0ϵc and therefore S 0 ϭct where we assume S(0)ϭ0.From the solvability condition for the equation at each order of the expansion we will obtain the corresponding corrections to the speed of the front.The solvability integral condition of Eq. ͑13͒ is † with null eigenvalue.The solvability integral condition for Eq.͑13͒ may be written as The integral in the numerator of Eq. ͑16͒ may be simplified by using Eq.͑12͒ and integrating by parts Finally one can obtain the first correction to the speed Note that in Eq. ͑17͒ there is no dependence on the solution of 0 but only on the function u .
The speed of the front reads, after inverting the hyperbolic scaling,

͑18͒
Before proceeding with the following order in the expansion it is necessary to solve Eq. ͑13͒.As L 1 (d 0 /dz)ϭ0 we look for a solution of the form 1 (z,t)ϭ(d 0 /dz) ϩ(d 0 /dz)zF(t) in Eq. ͑13͒, finding that After substituting S 0 ϭct, Eqs.͑17͒ and ͑19͒ into Eq.͑14͒, and applying the solvability condition ͓24͔ ͐ Ϫϱ ϱ e cz (d 0 /dz)L 1 ( 2 )dzϭ0 for Eq.͑14͒ we get and Note that in this case S ˙2 depends explicitly on the solution of 0 .In order to compute analytically the second-order correction of the speed it is necessary to have an analytical expression for the zeroth order solution 0 (z).Some exact solutions are known in the literature ͓1͔ for reaction terms of the form f ()ϭ qϩ1 (1Ϫ q ) for qу1.This source term has been applied to forest fires ͓25,26͔ and the spread of microorganisms ͓27͔.In this case, the solution for the homogeneous case takes the form It is easy to check that the integrals involved in Eqs.͑16͒ and ͑21͒ are convergent for any q.For example, for qϭ1 we have f ()ϭ 2 (1Ϫ), ␣ϭ1, and the speed of the front is given by where we have made use of Eqs.͑17͒ and ͑20͒.
For mathematical and numerical simplicity let us illustrate the above results for the case where u (x)ϭx is linear.Taking ϭ1 we have from Eq. ͑22͒ v(t)ϭ1/ͱ2ϩt/4ϩ 2 and inverting the hyperbolic scaling we obtain, for tӶO( Ϫ2 ), where the subscript SP stands for ''singular perturbation'' result.For any non-KPP f () one has The local speed approach assumes that the front position changes adiabatically in time as the front moves into a region where the characteristic parameters D and U change.For the first of Eqs.͑1͒ and for a source term f ()ϭ 2 (1Ϫ) the speed of the front is locally given by vϭͱDU(x f )/2, where x f is the position of the front.To be more specific, let us take also u (x)ϭx and ϭ1.In consequence, the dependence of speed of the front on the time is obtained by integrating the differential equation for the position of the front in dimensionless units.Taking x f (0)ϭ0 the local speed approach yields, for this case, where the subscript LA stands for ''local approach.''In Fig. 1 we compare, for different values of , the numerical results for the front speed of the first equation in Eqs.͑1͒ for U(x)ϭ1ϩ␦ u (x), u (x)ϭx, and (x,0), a Heaviside function, to the analytical solutions ͑23͒ and ͑26͒.We observe that v SP is in better agreement with numerical solutions than v LA , after an initial transient.This transient is due to two factors: it takes a certain interval of time for the front to fully develop and the asymptotic approximation is valid for u (x)ϭO(1).

B. Hamilton-Jacobi dynamics
The use of the Hamilton-Jacobi dynamics to study the front propagation is initially devoted to Freidlin ͓28͔ who treated the KPP minimal speed for slowly varying media using probabilistic ͑large deviation͒ approach but also rigorous mathematics has done by Ga ¨rtner ͓29͔ and Evans ͓30͔.However, as we will show in the last section, it is not necessary to assume either smooth or weak heterogeneities.We stress that singular perturbation analysis ͑preceding section͒ does not yield a fully analytical result for the very important KPP kinetic f ()ϭ(1Ϫ) ͓1,31͔ if one needs to go beyond first order in ␦, because the exact unperturbed solution is unknown and the solvability integrals diverge.In this section we determine the temporal evolution of the position of the reaction front for the logistic case.We replace (x/,t/) in Eq. ͑3͒ by an auxiliary field G(x,t)у0 through the exponential transformation ͑x,t͒ϭe ϪG(x,t)/ .͑27͒ We expect that (x,t) tends to a unit step function as →0.The equality G(x,t)ϭ0 determines the position of the front.Substituting Eq. ͑27͒ into Eq.͑2͒, for the KPP kinetics f ()ϭ(1Ϫ), we get, to leading order, the equation ( ϭ0) for the action functional where U(␦,x)ϭ1ϩ␦ u (x).From the analogy with the Hamilton-Jacobi equation ‫ץ‬ t G ϩH(‫ץ‬ x G,x)ϭ0, we define the Hamiltonian The solution for the action functional G(x,t) is given by G͑x,t ͒ϭ min x(sϭ0)ϭx,x(sϭt)ϭ0 ͵ and Finally, from Eq. ͑34͒ The position of the front given by G(x,t)ϭ0 is and the exact relationship for the speed after inverting the hyperbolic scaling is ͑36͒ for any ␦.For weak inhomogeneities (␦Ӷ1) one has for ϭ1, in Eq. ͑36͒, v HJ ͑ t ͒Ӎ2ϩ2t 2 ϩt 2 4 ϩO͑ 8 ͒, ͑37͒ which holds only for tӶO( Ϫ2 ).The subscript HJ stands for ''Hamilton-Jacobi'' result.
The local speed approach for the KPP kinetics yields v ϭ2ͱDU(x f ) and the differential equation for the position of the front is which after integrating under the initial condition x f (0)ϭ0 may be written as In Fig. 2 we compare Eqs.͑36͒ and ͑38͒ with the numerical solution for the first equation in Eq. ͑1͒ for different values of .After the initial transient, we observe, in general, good agreement.However, v HJ is in better agreement with numerical solutions than v LA , after the initial transient.
Another case with exact solution is u (x)ϭx 2 .In this case the Hamiltonian is the same as for the simple harmonic oscillator.Equation ͑30͒ with Eq. ͑31͒ yields FIG. 2. Comparison of the temporal evolution of the speed of fronts ͑in dimensionless units͒ between the Hamilton-Jacobi analytical result given in Eq. ͑36͒ ͑solid lines͒, the local speed approach ͑38͒ ͑dashed lines͒, and the numerical results ͑symbols͒ for different values of .This is the case u (x)ϭx and f ϭ(1Ϫ) for nonuniform reaction rate.
The energy integral and action functional are and the position of the front is given by From Eq. ͑39͒ one can see that the position as well as the speed of the front takes the infinite value just when t ϭ/4ͱ␦.However, for weak heterogeneities ͑␦Ӷ1͒ one has which does not have singularities but is valid only for t ӶO( Ϫ1 ).
If we assume weak heterogeneities ͑␦Ӷ1͒ we can approximate the speed of the front for any general (x).Details of the calculations are given in the appendix.It is important to note that Eq. ͑A7͒ is, up to ␦ order, equal to Eq. ͑18͒ for the KPP kinetic term where cϭ2 but differs for higher orders.

III. NONUNIFORM DIFFUSION
We consider now the problem where ␦ is the amplitude of the heterogeneities and D (x) is the heterogeneity function in diffusion.In this case, the local speed approach yields the same speed as for nonuniform reaction rate if D (x)ϭ u (x).As we will see below, this result is in disagreement with the singular perturbative and the Hamilton-Jacobi results.

A. Singular perturbation analysis
As in the preceding section, we assume ␦ϵ in Eq. ͑41͒ where ϭO(1) and we assume the existence of an outer and a boundary layer region.Equation ͑41͒, together with the corresponding boundary conditions, becomes In order to solve Eq. ͑42͒ in the outer region we use expansion ͑5͒.We can easily see that the effect of the het-erogeneity in the coefficient of diffusion does not affect the solution at least until O( 3 ).Thus (x,t;)ϭO( 3 ) to the left of the boundary layer and (x,t;)ϭ1ϩO( 3 ) to the right of the boundary layer.
In order to study the dynamics of Eq. ͑42͒ inside the boundary layer we substitute Eqs.͑9͒-͑11͒ into Eq.͑42͒ and collect terms with equal powers of .We obtain Eq. ͑12͒ for the lowest order, and for first and second orders, respectively.From the solvability conditions and assuming ϭ 0 (z) one has S ˙0ϭc constant and The speed of the front is finally given by It is very interesting to note that, in this case, up to second order in the speed correction does not depend on the solution of 0 .For D (x)ϭx, ϭ1, and general non-KPP kinetic term one has, for tӶO( Ϫ2 ), after inverting the hyperbolic scaling.In Fig. 3 we compare Eq. ͑44͒ for f ϭ 2 (1Ϫ) and Eq.͑26͒ with the numerical solution of the second equation in Eqs.͑1͒ with D(x)ϭ1 ϩ␦ D (x) and D (x)ϭx for different values of .In this case v LA is in slightly better agreement with numerical solutions than v SP , contrary to the previous case.

B. Hamilton-Jacobi dynamics
The Hamilton-Jacobi equation for problem ͑41͒, to leading order ͑ϭ0͒, with a KPP kinetic term is The Hamilton equations are The equation for x(s) is the corresponding integral energy is and the action functional is Let us now be more specific for the two cases where Eq.
͑46͒ has exact solution.First we take D (x)ϭx.From Eqs.
͑31͒ and ͑46͒ one has and from Eq. ͑48͒ one gets Therefore the speed is given, in an exact form, by after inverting the hyperbolic scaling.Note first of all that this result is equal to that obtained from the local speed approach.As in Sec.II, note that the result obtained in this section is essentially a leading order approximation while the result obtained in Sec.III A is an O() approximation.
In Fig. 4 we compare Eqs.͑49͒ ͑taking ␦ϭ) and ͑38͒ with the numerical result of the second equation in Eqs.͑1͒.
In this case the agreement between both analytical methods and the numerical results is not so good as for nonuniform reaction rate.
From the local speed approach the speed for a given dependence of the reaction rate on the spatial coordinate and constant diffusion coefficient is equal to the speed when the reaction rate is constant and the diffusion coefficient depends on the spatial coordinate with the same functional dependence as the above reaction rate.We have checked numerically for some values of that the speed for the linear dependence of the reaction rate is not equal, although very similar, to the speed for the linear dependence of the diffusion coefficient.
From the singular perturbative analysis and the Hamilton-Jacobi methods we can conclude that for the problems ‫ץ‬ t ϭ‫ץ‬ xx ϩ͓1ϩ␦ u (x)͔ f () and ϩ␦ D (x)͔‫ץ‬ x ͖ϩf() with non-KPP f (), such that f (0) ϭ f (1)ϭ0, where ␦, Ӷ1 ͑weak and smooth heterogene-ities͒ and u,D (x) is a continuous and derivable function, the speed of the front has the formal solution v͑ t ͒ϭcϩ 1 2 c u,D ͑ ct ͒␦ϩO͑ ␦ 2 ͒, for tӶO( Ϫ1 ) where c is the asymptotic ͑constant͒ speed for the homogeneous problem ‫ץ‬ t ϭ‫ץ‬ xx ϩ f ().If is an increasing/decreasing function of space, the front is accelerated/decelerated.

IV. FRACTAL MEDIA
In this section we illustrate the advantages of using the Hamilton-Jacobi method for dealing with heterogeneous media.In particular, we get an exact expression for the front speed propagation in fractal media.The reaction-diffusion process in a fractal may be described by the equation for the probability density of O'Shaughnessy and Procaccia ͓32͔ coupled to a KPP kinetic term where d is the dimension of the fractal, is an index which is 0 for the classical normal situation ͑Euclidean media͒, and D is a kind of diffusion coefficient.After taking into account the hyperbolic scaling r→r/ and t→t/ and the field G (r,t)ϭϪln(r,t), one has from Eq. ͑50͒, The first and third terms in the right-hand side of Eq. ͑51͒ have the same order of magnitude and in the asymptotic limit (→0) both terms may be neglected in front of the second term and in consequence, the Hamilton-Jacobi for the front propagation in a fractal media is where G(r,t)ϭlim →0 G (r,t).Proceeding as in the above sections one has and the exact expression for the speed of the front, once the hyperbolic scaling is inverted, is , ͑52͒ which describes a decelerated front.In Fig. 5 we compare the exact result ͑52͒ with the results of the numerical solutions of Eq. ͑50͒ for different values of .We stress that in this case of physical interest, in contrast to what would happen if using the local approach, we have not assumed weak or smooth heterogeneities because we have applied the Hamilton-Jacobi method.

V. CONCLUSIONS
We have studied how deterministic heterogeneities influence the front speed of parabolic reaction-diffusion equations where the reaction term and/or the diffusion coefficient depend on the spatial coordinates.We have derived analytic expressions for the speed of fronts that are valid for initially fully developed fronts or for more general initial conditions in the asymptotic limit.The singular perturbative analysis and the geometrical method of Hamilton-Jacobi have been employed to find the speed of the fronts propagating in deterministic heterogeneous media.The singular perturbative analysis has been used when spatial heterogeneities are weak (␦Ӷ1) and smooth (Ӷ1) and may be successfully applied only for fronts with non-KPP kinetics ͑pushed fronts͒.The expressions obtained for the speed are power series of , where secular terms appear and in consequence are not uniformly valid in time.However, for the simplest case of linear heterogeneities these expressions have been compared to numerical solutions exhibiting a good agreement.
The Hamilton-Jacobi method we used here only holds for fronts with KPP kinetics ͑pulled fronts͒.However, this method allows us to work without assuming either smooth or weak heterogeneities.We have compared the results for the simplest case of linear heterogeneities with numerical solutions and a good agreement is found again.Exact solution for the speed of fronts traveling in fractal media is obtained and compared to numerical solutions.We have found an excellent agreement and it has been shown both analytically and numerically that the front is decelerated.
Finally, the local speed approach has been compared with the above analytical methods and numerical solutions.For nonuniform diffusion coefficient this approach slightly im-FIG.5. Comparison of the temporal evolution of the speed of fronts ͑in dimensionless units͒ derived from the Hamilton-Jacobi method ͑52͒ ͑dashed lines͒ with numerical solutions of Eq. ͑50͒ ͑symbols͒ for fractal media.For ϭ0 one recovers the Fisher result.We have taken DϭUϭdϭ1.
proves the singular perturbative method and yields the same result found with the Hamilton-Jacobi method.However, for nonuniform reaction rate both singular perturbative and Hamilton-Jacobi are in better agreement with numerical results than the local speed approach.

APPENDIX: WEAK HETEROGENEITIES IN HAMILTON-JACOBI DYNAMICS
In this appendix we develop the calculations to obtain the speed of fronts for weak heterogeneities (␦Ӷ1) for both nonuniform reaction and nonuniform diffusion.

͑A5͒
The position xϭx(t) of the front comes from G(x,t)ϭ0 which has to be solved, by using the expansion x(t)ϭx 0 (t) ϩ␦x 1 (t)ϩ␦ 2 x 2 (t)ϩ O(␦ 3 ), to obtain and the speed is once the hyperbolic scaling is inverted.

Nonuniform diffusion coefficient
Proceeding analogously as in the preceding section the perturbed solution to Eq. ͑46͒ is given by x 0 ͑ s ͒ϭxϪxs/t, for 0ϽsϽt, and the energy integral and the action functional are Finally, from the temporal derivative of the position of the front given by Gϭ0 one gets

FIG. 1 .
FIG. 1.Comparison of the temporal evolution of the speed of fronts ͑in dimensionless units͒ between the singular perturbation analytical result given in Eq. ͑23͒ ͑solid lines͒, the local speed approach ͑26͒ ͑dashed lines͒, and the numerical results ͑symbols͒ for different values of .This is the case u (x)ϭx and f ϭ 2 (1 Ϫ) for nonuniform reaction rate.

FIG. 3 .
FIG.3.Comparison of the temporal evolution of the speed of fronts ͑in dimensionless units͒ between the singular perturbed analytical result given in Eq. ͑44͒ for cϭ1/ͱ2 ͑solid lines͒, the local speed approach ͑26͒ ͑dashed lines͒, and the numerical results ͑sym-bols͒ for different values of .This is the case D (x)ϭx and f ϭ 2 (1Ϫ) for nonuniform diffusion coefficient.