Linearly-implicit, approximate factorization ... - Semantic Scholar
Applied Mathematics and Computation 174 (2006) 1609–1633
www.elsevier.com/locate/amc
Linearly-implicit, approximate factorization, exponential methods for multi-dimensional reaction–diffusion equations J.I. Ramos Room I-320-D, E.T.S. Ingenieros Industriales, Universidad de Ma´laga, Plaza El Ejido, s/n 29013 Ma´laga, Spain
Abstract Linearly-implicit, approximate factorization methods are used to reduce the solution of multi-dimensional reaction–diffusion equations to a sequence of two-point boundaryvalue problems in ordinary differential equations which are solved analytically in a piecewise manner. It is shown that, depending on the allocation of the transient, diffusion and reaction terms to the second-order ordinary differential operators, several families of piecewise exponential methods can be obtained. When only diffusion is included, a quadratic analytical solution is obtained, whereas if transient and diffusion, diffusion and reaction, and transient, diffusion and reaction terms are included, piecewise exponential methods result. It is shown that exponential methods that only account for reaction and diffusion in the differential operator may result in ill-conditioned systems of algebraic equation, whereas, if diffusion, transient and reaction are included, ill-conditioning is avoided if the time step is sufficiently small. It is also shown that the most robust exponential method is that that accounts for transient and diffusion terms in the differential operator, and this method provides solutions in good accord with those
E-mail address: [email protected] 0096-3003/$ - see front matter 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2005.07.020
1610
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
of a linearly-implicit, approximate factorization technique which employs a three-point, second-order accurate discretization of the diffusion terms. It is also shown that a hybrid exponential method that solves the transient, reaction and diffusion terms at odd time levels and the transient and diffusion terms at even time levels also provides very accurate solutions. 2005 Elsevier Inc. All rights reserved. Keywords: RotheÕs methods; Approximate factorization; Linearly implicit techniques; Piecewise exponential solutions; Reaction–diffusion equations
1. Introduction Factorization methods are techniques that approximate a multi-dimensional operator by a sequence of one-dimensional ones. Since it is rarely possible to factorize exactly multi-dimensional operators, factorization techniques usually employ an iterative procedure to account for the factorization errors. If these errors are neglected, the resulting methods are referred to as approximate factorization techniques which transform a multi-dimensional nonlinear problem into a sequence of one-dimensional ones whose solution requires iterations. Iterations may be avoided by time linearization whereby the time variable in the original multi-dimensional operator is first discretized and the resulting nonlinear elliptic partial differential equation is then linearized with respect to the previous time level. Upon using approximate factorization into the resulting linear elliptic partial differential equation, one obtains a sequence of one-dimensional, linear, two-point boundary-value problems in each coordinate direction that may be discretized and written in terms of either the original dependent variables [1,2] or delta form [3,4]. These linearly-implicit, approximate factorization methods have been used with great success in the numerical solution of the Navier–Stokes and Euler equations [1–4], and nonlinear reaction–diffusion equations in one [5] and several spatial dimensions [6,7], where the diffusion terms have usually been discretized by means of second-order accurate finite difference formulae, although the techniques have also proved to be robust with three-point compact discretizations of the diffusion terms [8–10]. Since linearly-implicit, approximate factorization methods provide linear, two-point boundary-value problems in each spatial direction for reaction– diffusion equations, these problems can also be analyzed by means of locally one-dimensional techniques such as the diffusive and exponentially-fitted methods developed by the author [11] who considered one-dimensional timedependent advection–reaction–diffusion equations that upon time discretization and linearization resulted in second-order ordinary differential equations with non-constant coefficients. However, by assuming that the coefficients
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
1611
and right-hand sides of these equations are piecewise constant and integrating the resulting equations analytically, the author developed C0 and C1 exponentially-fitted methods that yield tridiagonal matrices. In C0 methods, the coefficients of the differential equation are assumed to be constant in the interval xi1 < x < xi+1 and the solution is assumed to be continuous at xi1, xi and xi+1, whereas, in C1 techniques, the coefficients of the differential equation are assumed to be constant in the interval xi1 < x < xi and the solution is assumed to be continuous and first-order differentiable at xi1 and xi. The advantage of C0 and C1 methods is that they can be used to study timedependent evolution problems such as those governed by parabolic differential equations and that these methods provide piecewise analytical solutions. By way of contrast, linearly-implicit, approximate factorization techniques that discretize the diffusion terms by either second-order or three-point, compact finite difference formulae provide only nodal values and require interpolation to determine an approximate continuous (in space) solution. On the other hand, C0 and C1 methods involve exponential terms and are more costly than standard finite difference discretizations of the diffusion terms. Recently, the author has applied C0 exponentially-fitted methods to a onedimensional, nonlinear reaction–diffusion equation [12] and found that these techniques may provide erroneous results when a Crank–Nicolson discretization is employed. In fact, in that reference, it was found that the discontinuities in the initial conditions disappeared rather slowly. Although these results are in qualitative agreement with those found by other authors who used standard finite difference discretizations for the diffusion terms [13,14], C0 methods were found to be rather sensitive to the time step employed in the calculations when the ordinary differential operators include transient, diffusion and reaction terms. For one-dimensional, one-phase Stefan problems that are characterized by the presence of a moving front, C0 techniques were also found to be rather sensitive to both the spatial step size and the time step [15], and it has been claimed that the application of C0 methods to multi-dimensional systems of reaction– diffusion equations may be plagued by serious numerical difficulties associated with the allocation of the transient, diffusion and reaction terms to the onedimensional differential operators, the spatial step size, the time step and the nonlinearity of the reaction terms [12]. It must be pointed out that the linearly-implicit, approximate factorization, exponential methods considered here (see also [12,15]) are, in fact, linear RotheÕs methods [16,17] since they are based on the discretization of the time variable while keeping continuous the spatial ones. Linear and nonlinear RotheÕs techniques have been previously used to study reaction–diffusion processes [12], one-dimensional, one-phase Stefan problems [15], two-dimensional advection–reaction–diffusion equations with boundary element methods [18], one-dimensional convection–diffusion equations with the sinc–Galerkin
1612
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
method [19], the one-dimensional nonlinear Schro¨dinger equation [20], reaction–diffusion problems with internal boundaries [21], the one-dimensional viscous Burgers equation with the Galerkin procedure [22], and the two-dimensional convection–reaction–diffusion equation [23]. Most of these studies, however, have used either standard finite difference or finite element methods for the first- and second-order spatial derivatives or locally one-dimensional methods, and considered only the convection and diffusion terms in the differential operator so that their exponential methods are able to deal with high Reynolds or Pe´clet number flows [24]. Some of these references have considered the solution of the locally one-dimensional operator to derive finite difference expressions. By way of contrast, the approach followed here considers different differential operators, provides three-point finite difference formulae upon imposing continuity requirements at three different grid points, and may consider transient, advection, diffusion and reaction terms in the differential operators. It should also emphasized that the exponential methods presented here are completely different from those derived from the discretization of the spatial coordinates, i.e., methods of lines, that yield systems of first-order, ordinary differential equations in time, which upon linearization result in systems of first-order, linear ordinary differential equations that have exponential solutions in time [25], even though the evaluation of the matrix exponential is a very hard task. The methods presented here provide piecewise exponential solutions in space. In this paper, we first (in Section 2) develop a variety of linearly-implicit, approximate factorization, exponential methods for the numerical solution of multi-dimensional reaction–diffusion equations that include as a subset the second-order standard finite difference discretization of the diffusion terms [5], and techniques that account for transient and diffusion, transient, diffusion and reaction, and reaction and diffusion terms in the differential operators. We then analyze the linear (Fourier–von Neumann) linear stability of the finite difference equations that result upon using the transient and diffusion terms in the differential operator, and comment on the difficulties that one encounters when trying to analyze the linear stability of the exponential methods that include reaction and diffusion, and transient, diffusion and reaction terms in the differential operators, in Section 3. In that section, we also point out the difficulties that are posed by using a Crank–Nicolson discretization both in terms of the linear stability and the characteristics of exponential techniques, and propose a remedy to the latter. In Section 4, the linearly-implicit, exponential, approximate factorization methods presented here are applied to one- and two-dimensional systems of nonlinear reaction–diffusion equations, and their numerical results are compared with those of linearly-implicit methods that employ second-order accurate discretizations of the diffusion terms [5,6]. Finally, a section on conclusions summarizes the main findings of this paper.
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
1613
2. Formulation The linearly-implicit, approximate factorization, exponential methods considered in this paper can be applied to one- and multi-dimensional problems. Since these methods are based on approximate factorization techniques that factorize a multi-dimensional operator into a sequence of one-dimensional ones, for the sake of conciseness, we shall only present the formulation for multi-dimensional problems. We, thus, consider the following system of twodimensional reaction–diffusion equations: oU o2 U o2 U ¼ 2 þ 2 þ SðU; t; x; yÞ; ot ox oy
ð1Þ
where S is a function in C1(RN · R · R · R ! RN), U 2 RN, x 2 R, y 2 R, and t 2 R, and the diffusion coefficients in the x- and y-directions have been set equal to unity. The independent variables t (0 6 t < 1), and x (0 6 x 6 Lx) and y (0 6 x 6 Ly) are the time and spatial coordinates, respectively, and Lx and Ly denote the domainÕs dimensions in the x- and y-directions, respectively. Although the methods presented here can be easily generalized to the case that the diffusion coefficient is a nonlinear tensor, we have adopted Eq. (1) to illustrate them because of its simplicity. Eq. (1) can be transformed into a system of nonlinear (elliptic) partial differential equations at each time step by discretizing the time variable by means of a h-method while keeping continuous the (independent) spatial variables. These nonlinear equations can, in turn, be transformed into linear ones by linearizing the nonlinear terms at tn+1 with respect to the previous time level, and the resulting equation may be expressed as Unþ1 Un o2 Unþ1 o2 U n o2 Unþ1 n n ¼h þ ð1 hÞ þ S þ hT k þ h k ox2 ox2 oy 2 þ ð1 hÞ
o2 Un þ hJn ðUnþ1 Un Þ; oy 2
ð2Þ
oF n n ðU ; t ; x; yÞ; oU
ð3Þ
where Tn ¼
oF n n ðU ; t ; x; yÞ; ot
Jn ¼
where k = tn+1 tn is the time step, Un = U(tn, x, y), and Sn = S(Un, tn, x, y). The values h = 0 and h = 1 correspond to first-order accurate (in time), explicit and implicit, respectively, methods, while h = 0.5 corresponds to a second-order accurate (in time) method. In this paper, we shall be interested on implicit techniques, and especially on 0.5 6 h 6 1.
1614
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
Eq. (3) may also be written as DU o2 DU o2 Un o2 DU o2 Un n n ¼h þ þ S þ hT k þ h þ þ hJn DU; k ox2 oy 2 ox2 oy 2
ð4Þ
where DU = Un+1 Un. Eq. (4) corresponds to a delta formulation and a coupled system of linear elliptic equations at each time level. This two-dimensional system may be reduced to sequences of one-dimensional equations by means of the following approximate factorization technique: o2 o2 I kh 2 I kdhJn I kh 2 I khJn DU ¼ RHS; ð5Þ ox oy where I denotes the unit or identity matrix, d þ ¼ 1; 2 n oU o2 Un n n RHS ¼ k þ þ S þ khT ox2 oy 2
ð6Þ ð7Þ
and the O(k2) factorization errors have been disregarded. Eq. (5) can also be written as the following system of uncoupled, ordinary differential equations which have non-constant coefficients due to the appearance of the Jacobian matrix, J, o2 n Lx ðDU Þ I kh 2 I kdhJ DU ¼ RHS; ð8Þ ox o2 Ly ðDUÞ I kh 2 I khJn DU ¼ DU ; ð9Þ oy which represent linear systems of one-dimensional operators in the x- and ydirections, respectively, and correspond to a delta formulation. These equations can also be written as o2 Lx ðU Þ I kh 2 I kdhJn U ¼ RHS; ð10Þ ox o2 n Ly ðUÞ I kh 2 I khJ Unþ1 ¼ U ; ð11Þ oy where 2 n oU o2 Un n n n n RHS ¼ k ð1 hÞ þ þ S hJ U þ khT . ox2 oy 2
ð12Þ
Note that in Eqs. (8) and (9) and Eqs. (10) and (11) the spatial coordinates have been kept continuous while time has been discretized; therefore, these equations correspond to a RotheÕs method [16].
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
1615
Since Eqs. (8)–(11) have the same form, hereon, we shall only work with only one-dimensional operator, e.g., kh
o2 U þ ðI khdHÞU ¼ R; ox2
ð13Þ
where U, H and R can be readily obtained from Eqs. (8) and (10), or Eqs. (9) and (11) if x and d in Eq. (13) are replaced for y and , respectively.
3. Numerical methods In this section, we consider a variety of methods for the numerical solution of Eq. (13). These methods are presented in an integrated fashion so that the reader can appreciate the differences in how to obtain them. 3.1. Diffusion operator method In this method, Eq. (13) is first written as kh
o2 U ¼ P R ðI khdHÞU; ox2
ð14Þ
which, upon assuming that P is constant in xi1 6 x 6 xi+1 and imposing the following conditions: Uðxi1 Þ ¼ Ui1 ;
Uðxi Þ ¼ Ui ;
Uðxiþ1 Þ ¼ Uiþ1 ;
yields the following three-point finite difference equation: kh kh kh 2 Ui1 þ 1 þ 2 2 I khdHi Ui 2 Uiþ1 ¼ Ri ; h h h
ð15Þ
ð16Þ
where h = xi+1 xi = xi xi1, which is exactly the same equation as the one that can be derived by discretizing the second-order derivative in Eq. (13) by means of second-order accurate finite differences in an equally-spaced grid [5]. Since U is a vector, Eq. (16) can be written as a system of linear algebraic equations where the matrix has a block tridiagonal structure. Hereon, we shall refer to the method represented by Eq. (16) as D. 3.2. Exponential transient–diffusion operator method In this method, Eq. (13) is first written as kh
o2 U þ U ¼ P R þ khdHU; ox2
ð17Þ
1616
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
which only includes the transient and diffusion operators in the left-hand side and, for h 5 0, can be integrated analytically in xi1 6 x 6 xi+1 by assuming that P is constant in this interval, to yield U ¼ A coshðkðx xi ÞÞ þ B sinhðkðx xi ÞÞ þ W;
ð18Þ
where A and B are constant vectors, k2 ¼ kh1 and W = Pi. Upon using Eq. (15), Eq. (18) yields the following three-point finite difference method: Ui1 2ðcoshðkhÞI þ khdð1 coshðkhÞÞHi ÞUi þ Uiþ1 ¼ 2Ri ð1 coshðkhÞÞ;
ð19Þ
which can also be written as a block tridiagonal system of linear algebraic equations. Hereon, we shall refer to the method represented by Eq. (19) as ETD. 3.3. Exponential reaction–diffusion operator method In this method, Eq. (13) is first written as the following two equations for N = 2: kh
o2 u khdH 11 u ¼ P 1 R1 u þ khdH 12 v; ox2
ð20Þ
kh
o2 v khdH 22 v ¼ P 2 R2 v þ khdH 21 u; ox2
ð21Þ
which includes the diffusion and part of the reaction terms in the left-hand side, where U = (u, v)T and the superscript T denotes transpose. For h 5 0, Eq. (20) or Eq. (21) can be integrated analytically in xi1 6 x 6 xi+1 by assuming that P is constant in this interval. For example, the solution of Eq. (20) subject to Eq. (15) yields h2 h2 uiþ1 þ ui1 2 þ ð22Þ ui þ dh2 ðH 12 Þi vi ¼ ðR1 Þi kh kh for H11 = 0, uiþ1 þ ui1 2 cosðkhÞ þ þ2
1 ð1 cosðkhÞÞ ui khdH 11
ðH 12 Þi ðR1 Þi ð1 cosðkhÞÞvi ¼ 2 ðcosðkhÞ 1Þ ðH 11 Þi khdðH 11 Þi
for dH11 > 0, where k2 = dH11, and
ð23Þ
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
uiþ1 þ ui1 2 coshðkhÞ þ þ2
1617
1 ð1 coshðkhÞÞ ui khdH 11
ðH 12 Þi ðR1 Þi ð1 coshðkhÞÞvi ¼ 2 ðcoshðkhÞ 1Þ ðH 11 Þi khdðH 11 Þi
ð24Þ
for dH11 < 0, where k2 = dH11. Eqs. (22)–(24) and the ones corresponding to v can be written as a block tridiagonal system of linear algebraic equations. Hereon, we shall refer to the method represented by Eqs. (22)–(24) as ERD. 3.4. Exponential transient–reaction–diffusion operator method In this method, Eq. (13) is first written as the following two equations for N = 2: kh
o2 u þ ð1 khdH 11 Þu ¼ P 1 R1 þ khdH 12 v; ox2
ð25Þ
kh
o2 v þ ð1 khdH 22 Þv ¼ P 2 R2 þ khdH 21 u; ox2
ð26Þ
which includes the transient term, diffusion and part of the reaction terms in the left-hand side, where U = (u, v)T and the superscript T denotes transpose. For h 5 0, Eq. (25) or Eq. (26) can be integrated analytically in xi1 6 x 6 xi+1 by assuming that P is constant in this interval. For example, the solution of Eq. (25) subject to Eq. (15) yields uiþ1 þ ui1 2ui þ h2 dðH 12 Þi vi ¼
h2 ðR1 Þi kh
ð27Þ
for Q1 = 1 khdH11 = 0, uiþ1 þ ui1 2 cosðkhÞui 2ð1 cosðkhÞÞ ¼ 2
khdðH 12 Þi vi 1 khdðH 11 Þi
ðR1 Þi ðcosðkhÞ 1Þ 1 khdðH 11 Þi
ð28Þ
for Q1 > 0, where k2 ¼ dH 11 kh1 , and uiþ1 þ ui1 2 coshðkhÞui 2 ¼ 2
khdðH 12 Þi ð1 coshðkhÞÞvi 1 khdðH 11 Þi
ðR1 Þi ðcoshðkhÞ 1Þ 1 khdðH 11 Þi
for Q1 < 0, where k2 ¼ dH 11 þ kh1 .
ð29Þ
1618
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
Eqs. (27)–(29) and the ones corresponding to v can be written as a block tridiagonal system of linear algebraic equations. Hereon, we shall refer to the method represented by Eqs. (27)–(29) as ETRD.
4. Linear stability In this section, we examine the linear stability of the methods D and ETD presented in the previous section for one-dimensional scalar problems. We, thus, consider the following partial differential equation: ou o2 u ¼ þ au; ot ox2 where a is a constant. Eq. (16) can be written taking into account Eq. (19) as kh nþ1 kh kh nþ1 2 ui1 þ 1 þ 2 2 kha uinþ1 2 uiþ1 h h h k 2k ¼ ð1 hÞ 2 ðuiþ1 þ ui1 Þ þ 1 þ ð1 hÞ ak 2 uni . h h
ð30Þ
ð31Þ
Using uni ¼ /n expðIiKhÞ where I2 = 1 in Eq. (31), one can readily obtain the following amplification factor: 2 4k 1 þ ð1 hÞ ka sin l 2 /nþ1 h ; ð32Þ ¼ 4kh /n 1 kha þ h2 sin2 l where l ¼ Kh and K is the wave number. 2 A sufficient condition for linear stability is that a 6 0 for h = 0.5 and 1. The discretization of the time variable in Eq. (30) by means of a h method yields the following ordinary differential equation: 2 n o2 unþ1 ou nþ1 n n kh þ ð1 khaÞu ¼ u þ kð1 hÞ þ au . ð33Þ ox2 ox2 If we apply to this equation the exponential techniques presented in the previous section, we may have to deal with the second-order spatial derivative in the right-hand side. For example, the use of ETD yields nþ1 nþ1 ui1 2ðcoshðkhÞ þ khað1 coshðkhÞÞunþ1 þ uiþ1 i 2 n o ui n n ¼ 2 kð1 hÞ þ au þ u i i ð1 coshðkhÞÞ; ox2
ð34Þ
which yields a tridiagonal system of linear algebraic equations, where k2 ¼ kh1 . A sufficient condition for this system to be diagonally dominant is a 6 0.
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
1619
Eq. (34) poses a problem when trying to analyze its linear stability because of the presence of the second-order spatial derivative, unless h = 1. Note that this derivative is evaluated at the previous time level and can be determined at t = 0 from the initial conditions, and, at subsequent times, from the piecewise analytical solution determined in Section 2 as follows. The analytical solution used in ETD can be written as (cf. Eq. (18)) unþ1 ¼ Ai coshðkðx xi ÞÞ þ Bi sinhðkðx xi ÞÞ þ P i ;
xi1 6 x 6 xiþ1 ; ð35Þ
where 2 n ou n P ¼ u þ kð1 hÞ þ au þ khaunþ1 ox2 n
ð36Þ
and Ai and Bi can be determined from Eqs. (15) and (35) as Ai ¼ uinþ1 P i ;
Bi ¼
nþ1 uiþ1 unþ1 iþ1 sinhðkhÞ
ð37Þ
2 nþ1
and, therefore, o oxu 2 ¼ k2 ðAi coshðkðx xi ÞÞ þ Bi sinhðkðx xi ÞÞÞ can be easily determined. Note, however, that Ai and Bi depend on the solution at the previous time level through Eqs. (36) and (37). Moreover, in ERD and ETRD, k is also a function of the previous time level and this complicates the stability analysis of these techniques. In addition, it is easily shown that ETD, ERD and ETRD result in matrices that are diagonally dominant if a is negative (cf. Eqs. (19), (24) and (29), respectively.) In view of these comments, in the next subsection, we analyze the linear stability of TD for h = 1. 4.1. Stability of the transient–diffusion method For h = 1, the use of uni ¼ /n expðIiKhÞ, where I2 = 1, in Eq. (34) yields the following amplification factor: /nþ1 1 coshðkhÞ ¼ cosð2lÞ ðcoshðkhÞ þ khað1 coshðkhÞÞ /n ¼
sinh2 b ; sin2 l þ ð1 khaÞsinh2 b
ð38Þ
where b ¼ kh2 , and this factor has an absolute value equal to or less than one if a 6 0, therefore, ETD is linearly stable according to the Fourier–von Neumann analysis presented here. It must be noted that, for kh 1, a TaylorÕs series expansion of Eq. (34) yields Eq. (31) up to O((kh)4), for h = 1.
1620
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
5. Results In this section, we present some sample results for one- and two-dimensional nonlinear reaction–diffusion equations obtained with the methods presented in this paper. 5.1. One-dimensional reaction–diffusion equations The one-dimensional reaction–diffusion system considered in this subsection can be written as ou o2 u ¼ uv; ot ox2 ov o2 v 1 ¼ þ uv v; ot ox2 2
ð39Þ ð40Þ
subject to u(0, x) = 1, v(0, x) = exp(x2), u(t, 200) = u(t, 200) = 1 and v(t, 200) = v(t, 200) = 0. This problem was solved in an equally-spaced grid consisting of NI = 800 grid points, h = 0.5 and several time steps, and some sample results are presented in Figs. 1 and 2. Fig. 1 corresponds to k = 0.05 and shows the solution obtained with D. The results presented in Fig. 1 indicate that the initial peak in v becomes two travelling waves in time that generate a widening valley in u. The errors of ETD relative to D with k = 0.01, i.e., duðkÞ ¼ uD ðk ¼ 0:01Þ uETD ðkÞ and dvðkÞ ¼ vD ðk ¼ 0:01Þ vETD ðkÞ, are shown in Figs. 2 and 3 as functions of space and time, and the time step. Figs. 2 and 3 indicate that the relative errors of both u and v increase as time is increased and are located at the inflection points of the u profile and at the relative maxima of v. The largest errors du and dv of ETD at t = 100 are less than 1 · 102 and 4 · 103, 4 · 103 and 2 · 103, 15 · 103 and 5 · 103, and 0.02 and 0.01, respectively, for k = 0.05, 0.025, 0.1 and 0.01, respectively.
Fig. 1. u (left) and v (right) as functions of x and t obtained with D, k = 0.05, NI = 800.
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
1621
Fig. 2. EU ðkÞ ¼ uD ðk ¼ 0:01Þ uETD ðkÞ (left) and EV ðkÞ ¼ vD ðk ¼ 0:01Þ vETD ðkÞ (right) as functions of x and t obtained with NI = 800. (Top: k = 0.05; bottom: k = 0.025).
On the other hand, the errors duðkÞ ¼ uD ðkÞ uETD ðkÞ and dvðkÞ ¼ vD ðkÞ vETD ðkÞ of ETD were found to decrease as the time step is decreased and, for k = 0.01 and 0.1, are, at most, 0.02 and 0.01, and 0.1 and 0.15, respectively, for u and v, respectively. In Tables 1 and 2, we show the relative errors of ETD with respect to D, as well as the relative errors of the later with respect to itself as functions of the time step. Although not shown here, the errors uD ðk ¼ 0:01Þ uD ðkÞ and vD ðk ¼ 0:01Þ vD ðkÞ exhibit the same trends as those presented in Figs. 2 and 3. The method ETRD provided nearly identical results to those shown in Figs. 2 and 3 and, for this reason, they are not presented here. Note that for jkhd(H11)ij 1, the denominators that appear in Eqs. (28) and (29) can be approximated by one, and Eq. (29) tends to Eq. (19). ERD was found to be not as robust as ETD for the grid size and time steps mentioned above. In fact, for these time steps, this method was found to incur in errors when evaluating the sine and cosine terms that appear in their finite difference equations (cf. Eq. (23)). As one can easily check, ERD may not result in diagonally dominant matrices for k = 0.05, 0.025 and 0.1, and, in some of the calculations presented here, this method was found to provide ill-conditioned matrices, especially for NI = 800. This ill-conditioning was still observed
1622
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
Fig. 3. EU ðkÞ ¼ uD ðk ¼ 0:01Þ uETD ðkÞ (left) and EV ðkÞ ¼ vD ðk ¼ 0:01Þ vETD ðkÞ (right) as functions of x and t obtained with NI = 800. (Top: k = 0.1; bottom: k = 0.01).
Table 1 Maximum errors in u of ETD with respect to D, i.e., E ¼ uD uETD , and e ¼ uD ðk ¼ 0:01Þ uD ðkÞ as functions of k Error
k = 0.05
k = 0.025
k = 0.1
k = 0.01
E e
0.15 0.15
0.06 0.06
0.3 0.3
0.02 0
Table 2 Maximum errors in v of ETD with respect to D, i.e., E ¼ vD vETD , and e ¼ vD ðk ¼ 0:01Þ vD ðkÞ as functions of k Error
k = 0.05
k = 0.025
k = 0.1
k = 0.01
E e
0.04 0.05
0.02 0.02
0.15 0.1
0.01 0
in ERD even for k = 0.01, 0.005, 0.0025 and 0.001 and is also due to the fact that, as indicated in Fig. 1, H11 = v is nil wherever v = 0 and this implies that the solution behaves as a quadratic approximation for H11 wherever v = 0 and as a trigonometric function wherever dH11 > 0.
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
1623
As will be shown below for two-dimensional reaction–diffusion equations, ETRD behaves very well and ERD is not robust due to the reasons exposed above, whereas a combination of ETRD and TD may yield good results. 5.2. Two-dimensional reaction–diffusion equations The two-dimensional reaction–diffusion system considered in this subsection can be written as ou o2 u o2 u ¼ þ u2 v; ot ox2 oy 2
ð41Þ
ov o2 v o2 v 1 ¼ þ þ u2 v v; ot ox2 oy 2 2
ð42Þ
subject to u(0, x, y) = 1, v(0, x, y) = exp((x2 + y2)), u(t, 20, x) = u(t, 20, y) = u(t,x,20) = u(t,x,20) = 1 and v(t,20,x) = v(t,20,y) = v(t,x,20) = v(t,x,20) = 0. This problem was solved in an equally-spaced grid consisting of NI = NJ = 102 grid points in the x- and y-coordinates, respectively, h = 0.5 and several time steps, and some sample results are presented in Figs. 4 and 5. Fig. 4 shows u obtained with D as a function of space at different times
Fig. 4. u as a function of x and y obtained with D, k = 0.02, NI = NJ = 102. (From left to right, top to bottom: t = 4, 6, 8, 10, 12, 14, 16, 18, 20).
1624
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
Fig. 5. v as a function of x and y obtained with D, k = 0.02, NI = NJ = 102. (From left to right, top to bottom: t = 4, 6, 8, 10, 12, 14, 16, 18, 20).
and indicates that u exhibits a valley in the middle of the computational domain which is generated by the reaction term and the initial condition for v. Both the depth and the width of this valley increase with time until about t = 20 when u reaches an almost stationary value. Fig. 5 shows that, owing to the initial conditions, v develops initially a peak in the middle of the domain. The height of this peak decreases, whereas its width increases as time increases. At about t = 8, the peak flattens and a valley is formed. The depth of this valley increases as time increases until about t = 20 when v reaches an almost stationary profile characterized by high values near the domain boundaries. In order to assess the accuracy of the linearly-implicit, approximate factorization, exponential methods presented in this paper, we have performed comparisons with D and defined the local errors as EU M ¼ uD ðx; y; tÞ uM ðx; y; tÞ and EV M ¼ vD ðx; y; tÞ vM ðx; y; tÞ, where the subscripts D and M denote the solutions obtained with D and M, respectively, and M = ETRD, ETD and RD. Figs. 6 and 7 show the errors of ETRD and ETD, respectively, as functions of space at different times for k = 0.02 and h = 0.5. The errors of ETRD are, at most, 1 · 103 at t = 4, 6, 8, 10 and 12, and, at most, 2 · 103 at t = 14, 16, 18 and 20, and occur where either juj or its gradient are largest as indicated in Fig. 6. The errors of ETD illustrated in Fig. 7 exhibit the same trends as those of ETRD and are, at most, 1 · 103 at t = 4, 6 and 8, at most, 2 · 103 at t = 10, 12 and 14, and, at most, 4 · 103 at t = 16, 18 and 20, i.e., they are larger than those of ETRD.
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
1625
Fig. 6. Errors EU ETRD ¼ du ¼ uD uETRD , where D and ETRD denote the solutions obtained with D and ETRD, respectively, k = 0.02, NI = NJ = 102. (From left to right, top to bottom: t = 4, 6, 8, 10, 12, 14, 16, 18, 20).
Although not shown here, it was observed that ERD exhibited similar problems to those found for one-dimensional systems of reaction–diffusion equations, i.e., it provided ill-conditioned matrices. On the other hand, it was also found that a combination of ETRD and ETD which uses ETRD and ETD in the x- and y-operators, respectively, at odd time levels, and in the y- and x-operators, respectively, at even time levels, respectively, provides the errors shown in Fig. 8. Hereon, we shall refer to the hybrid method just described as ETRD–TD. Fig. 8 shows that the errors of ETRD–TD are, at most, 1 · 103 at t = 4, 6, 8 and 10, at most 2 · 103 at t = 12, 14 and 16, and, at most, 4 · 103 at t = 18 and 20. Therefore, the errors of ETRD–TD are higher and smaller, respectively, than those of ETRD and ETD, respectively. Figs. 9–11 illustrate the errors in v of ETRD, ETD and ETRD–TD, respectively, as functions of space at different times for k = 0.02 and h = 0.5. The errors EV ETRD are, at most, 1 · 103 at t = 4, 14, 16, 18 and 20, and 4 · 104 at t = 6, 8, 10 and 12. The errors EV ETD and EV ETRD–TD are, at most, 1 · 103 at t = 4, 10, 12, 14, 16, 18 and 20, and 4 · 104 at t = 6 and 8. The effects of the time step on the accuracy of ETRD, ETD and ETRD–TD relative to D are presented in Tables 3 and 4 for u and v, respectively. Tables 3 and 4 indicate that the relative accuracy increases as k is decreased, ETRD is
1626
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
Fig. 7. Errors EU ETD ¼ du ¼ uD uETD , where D and ETD denote the solutions obtained with D and ETD, respectively, k = 0.02, NI = NJ = 102. (From left to right, top to bottom: t = 4, 6, 8, 10, 12, 14, 16, 18, 20).
Fig. 8. Errors EU ETRDTD ¼ du ¼ uD uETRDTD , where D and ETRD–TD denote the solutions obtained with D and ETRD–TD, respectively, k = 0.02, NI = NJ = 102. (From left to right, top to bottom: t = 4, 6, 8, 10, 12, 14, 16, 18, 20).
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
1627
Fig. 9. Errors EV ETRD ¼ dv ¼ vD vETRD , where D and ETRD denote the solutions obtained with D and ETRD, respectively, k = 0.02, NI = NJ = 102. (From left to right, top to bottom: t = 4, 6, 8, 10, 12, 14, 16, 18, 20).
more accurate than ETRD–TD which, in turn, is more accurate than ETD. This means that the accuracy of linearly-implicit, approximate factorization, exponential methods that include the reaction, transient and diffusion terms when determining the piecewise analytical solution of the one-dimensional operators is higher than that of exponential techniques that only include the transient and diffusion terms in these operators. On the other hand, one-dimensional operators that only include diffusion and reaction may result in ill-conditioned matrices. Figs. 6–11 show that the linearly-implicit, approximate factorization, exponential methods presented in this paper preserve the symmetries of the original problem (cf. Eqs. (41) and (42)). Although the results presented in this paper indicate that linearly-implicit, exponential methods have an accuracy comparable to those of linearly-implicit techniques that use standard second-order accurate finite difference formulae and, therefore, provide piecewise analytical solution and there is no need for interpolation, they may be subject to overflow for large spatial step sizes and/or small time steps because the solution of the one-dimensional operators depends in an exponential manner on the ratio phffiffik. Furthermore, when the initial conditions are not C2 and h = 0.5, linearly-implicit, exponential are not
1628
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
Fig. 10. Errors EV ETD ¼ dv ¼ vD vETD , where D and ETD denote the solutions obtained with D and ETD, respectively, k = 0.02, NI = NJ = 102. (From left to right, top to bottom: t = 4, 6, 8, 10, 12, 14, 16, 18, 20). 2
2
strictly applicable because they do require the values of ooxU2 and/or ooyU2 of the initial conditions, at the first time level (cf. Eqs. (7) and (12)). When this occurs, these second derivatives may be evaluated approximately by means of, for example, second-order accurate finite difference formulae which correspond to a quadratic (rather than exponential) approximation.
6. Conclusions Linearly-implicit, approximate factorization, exponential methods for the numerical solution of multi-dimensional reaction–diffusion equations have been developed. These techniques are based on the discretization of the time derivatives and the time linearization of the nonlinear terms, and yield linear elliptic partial differential equations at each time level. Upon factorizing the multi-dimensional elliptic operator into a sequence of one-dimensional operators while keeping the spatial coordinates continuous, two-point, linear boundary-value problems are obtained. These boundary-value problems are characterized by non-constant coefficients; however, they can be linearized (in space) in a piecewise fashion to yield linear, constant coefficients, second-order ordinary differential equations that can be solved in closed form to obtain
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
1629
Fig. 11. Errors EV ETRD–TD ¼ dv ¼ vD vETRD–TD , where D and ETRD–TD denote the solutions obtained with D and ETRD–TD, respectively, k = 0.02, NI = NJ = 102. (From left to right, top to bottom: t = 4, 6, 8, 10, 12, 14, 16, 18, 20).
piecewise analytical solutions which, upon imposing, continuity conditions at three different grid points yield block tridiagonal systems of linear algebraic equations. Depending on the allocation of the transient, diffusion and reaction terms to the second-order ordinary differential equations, several exponential techniques have been developed. If the differential operator only includes diffusion, these methods provide piecewise quadratic solutions and their finite difference equations are identical to those that result from a second-order accurate finite difference discretization in space. If transient and diffusion, transient, diffusion and reaction, and diffusion and reaction terms are included in the differential operators, these techniques provide exponential finite difference equations that depend on the diffusion, diffusion and reaction, and the reaction times, respectively. It has been show that, for the one- and two-dimensional reaction–diffusion problems considered in this paper, exponential methods that account for transient and diffusion, and transient, diffusion and reaction terms in the differential operators provide solutions of nearly the same accuracy as that of second-order (in space) accurate, linearly-implicit techniques. However, these techniques may be subject to overflow when the inverse of the Fourier number is very large, due to the exponential character of the solution. On the other
1630
Error
t=4
t=6
t=8
t = 10
t = 12
t = 14
t = 16
t = 18
t = 20
EU1(k = 0.04) EU1(k = 0.02) EU1(k = 0.01)
3
2 · 10 1 · 103 4 · 104
3
2 · 10 1 · 103 4 · 104
3
2 · 10 1 · 103 4 · 104
3
2 · 10 1 · 103 4 · 104
3
2 · 10 1 · 103 1 · 103
3
2 · 10 2 · 103 1 · 103
3
4 · 10 2 · 103 1 · 103
3
4 · 10 2 · 103 1 · 103
4 · 103 2 · 103 1 · 103
EU2(k = 0.04) EU2(k = 0.02) EU2(k = 0.01)
2 · 103 1 · 103 4 · 104
2 · 103 1 · 103 4 · 104
2 · 103 1 · 103 1 · 103
4 · 103 2 · 103 1 · 103
4 · 103 2 · 103 1 · 103
4 · 103 2 · 103 2 · 103
5 · 103 4 · 103 2 · 103
5 · 103 4 · 103 2 · 103
1 · 102 4 · 103 2 · 103
EU3(k = 0.04) EU3(k = 0.02) EU3(k = 0.01)
2 · 103 1 · 103 4 · 104
2 · 103 1 · 103 4 · 104
2 · 103 1 · 103 4 · 104
2 · 103 1 · 103 1 · 103
4 · 103 2 · 103 1 · 103
4 · 103 2 · 103 1 · 103
4 · 103 2 · 103 2 · 103
4 · 103 4 · 103 2 · 103
4 · 103 4 · 103 2 · 103
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
Table 3 Errors in u of ETRD, ETD and ETRD–TD with respect to D, i.e., EU M ¼ uD uM , where M = 1, 2 and 3 correspond to ETRD, ETD and ETRD–TD, respectively
Error
t=4
t=6
t=8
t = 10
t = 12
t = 14
t = 16
t = 18
t = 20
EV1(k = 0.04) EV1(k = 0.02) EV1(k = 0.01)
3
2 · 10 1 · 103 4 · 104
3
1 · 10 4 · 104 4 · 104
3
1 · 10 4 · 104 2 · 104
3
1 · 10 4 · 104 2 · 104
3
1 · 10 4 · 104 2 · 104
3
1 · 10 1 · 103 1 · 103
3
1 · 10 1 · 103 1 · 103
3
1 · 10 1 · 103 1 · 103
1 · 103 1 · 103 1 · 103
EV2(k = 0.04) EV2(k = 0.02) EV2(k = 0.01)
1 · 103 1 · 103 4 · 104
1 · 103 4 · 104 2 · 104
1 · 103 4 · 104 2 · 104
1 · 103 1 · 103 4 · 104
2 · 103 1 · 103 4 · 104
2 · 103 1 · 103 1 · 103
2 · 103 1 · 103 1 · 103
2 · 103 1 · 103 1 · 103
2 · 102 1 · 103 1 · 103
EV3(k = 0.04) EV3(k = 0.02) EV3(k = 0.01)
2 · 103 1 · 103 4 · 104
1 · 103 4 · 104 2 · 104
1 · 103 4 · 104 2 · 104
1 · 103 4 · 104 4 · 104
1 · 103 1 · 103 4 · 104
1 · 103 1 · 103 1 · 103
2 · 103 1 · 103 1 · 103
2 · 103 1 · 103 1 · 103
2 · 103 1 · 103 1 · 103
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
Table 4 Errors in v of ETRD, ETD and ETRD–TD with respect to D, i.e., EV M ¼ vD vM , where M = 1, 2 and 3 correspond to ETRD, ETD and ETRD–TD, respectively
1631
1632
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
hand, exponential methods that only account for the diffusion and reaction terms were found to provide ill-conditioned matrices. The exponential methods presented in this paper can be easily applied to fully implicit methods, i.e., h = 1, when these techniques employ a primitive variable form. However, if either a primitive variable formulation for h = 0.5 or a delta form for 0 < h 6 1 is employed, these exponential techniques require that the solution be continuous with continuous second-order (spatial) derivatives because of the appearance of these derivatives in the right-hand side of the ordinary differential operators. Otherwise, these techniques are not applicable in their exponential formulation, unless the second-order derivatives at the previous time level are discretized by means of second-order accurate finite difference formulae.
Acknowledgement The research reported in this paper was partially supported by Project BFM2001-1902 from the Ministerio de Ciencia y Tecnologı´a of Spain and fondos FEDER.
References [1] W.R. Briley, H. McDonald, On the structure and use of linearized block implicit methods, J. Comput. Phys. 34 (1980) 54–73. [2] W.R. Briley, H. McDonald, An overview and generalization of implicit Navier–Stokes algorithms and approximate factorization, Comput. Fluids 30 (2001) 807–828. [3] R.M. Beam, R.F. Warming, An implicit finite-difference algorithm for hyperbolic systems in conservation-law form, J. Comput. Phys. 22 (1976) 87–110. [4] R.M. Beam, R.F. Warming, An implicit factored scheme for the compressible Navier–Stokes equations, AIAA J. 16 (1978) 393–402. [5] J.I. Ramos, Linearization methods for reaction–diffusion equations: 1-D problems, Appl. Math. Comput. 88 (1997) 199–224. [6] J.I. Ramos, Linearization methods for reaction–diffusion equations: multidimensional problems, Appl. Math. Comput. 88 (1997) 225–254. [7] J.I. Ramos, Linearized factorization techniques for multidimensional reaction–diffusion equations, Appl. Math. Comput. 100 (1999) 201–222. [8] J.I. Ramos, Implicit, compact, linearized h-methods with factorization for multidimensional reaction–diffusion equations, Appl. Math. Comput. 94 (1998) 17–43. [9] Y. Gu, W. Liao, J. Zhu, An efficient high-order algorithm for solving systems of 3-D reaction– diffusion equations, J. Comput. Appl. Math. 155 (2003) 1–17. [10] W. Liao, J. Zhu, A.Q.M. Khaliq, An efficient high-order algorithm for solving systems of reaction–diffusion equations, Numer. Meth. Partial Diff. Equat. 18 (2002) 340–354. [11] J.I. Ramos, On diffusive methods and exponentially-fitted techniques, Appl. Math. Comput. 103 (1999) 69–96. [12] J.I. Ramos, Exponential methods for one-dimensional reaction–diffusion equations, Appl. Math. Comput., in press, doi:10.1016/j.amc.2004.12.003.
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
1633
[13] D. Britz, O. Østerby, J. Strutwolf, Damping of Crank–Nicolson error oscillations, Comput. Biol. Chem. 27 (2003) 253–263. [14] O. Østerby, Five ways of reducing the Crank–Nicolson oscillations, BIT 43 (2003) 811–822. [15] J.I. Ramos, Exponential numerical methods for one-dimensional, one-phase Stefan problems, Arch. Appl. Mech., in press. [16] E. Rothe, Zweidimensionale parabolische Randwertaufgaben als eindimensionaler Randwertaufgaben, Math. Ann. 102 (1930) 652–670. [17] A. Ashyralyev, P.E. Sobolevskii, Well-posedness of Parabolic Difference Equations, Birkha¨user Verlag, Basel, Switzerland, 1994. [18] N. Samec, L. Sˇkerget, Integral formulation of a diffusive–convective transport equation for racting flows, Eng. Anal. Bound. Elem. 28 (2004) 1055–1060. [19] J.L. Mueller, T.S. Shores, A new sinc–Galerkin method for convection–diffusion equations with mixed boundary conditions, Comput. Math. Appl. 47 (2004) 803–822. [20] F. Schmidt, D. Yevick, Discrete boundary conditions for Schro¨dinger–type equations, J. Comput. Phys. 134 (1997) 96–107. [21] J. Lang, High–resolution self–adaptive computations on chemical reaction–diffusion problems with internal boundaries, Chem. Eng. Sci. 51 (1996) 1055–1070. ¨ zdesß, A numerical solution of BurgersÕ equation, Appl. Math. Comput. 156 [22] E.N. Aksan, A. O (2004) 395–402. [23] T.W.H. Sheu, S.K. Wang, R.K. Lin, An implicit scheme for solving the convection–diffusion– reaction equation in two dimensions, J. Comput. Phys. 164 (2000) 123–142. [24] K.W. Morton, Numerical Solution of Convection–Diffusion Problems, Chapman & Hall, London, UK, 1996. [25] P.J. van der Houwen, B.P. Sommeijer, Approximate factorization for time-dependent partial differential equations, J. Comput. Appl. Math. 128 (2001) 447–466.
www.elsevier.com/locate/amc
Linearly-implicit, approximate factorization, exponential methods for multi-dimensional reaction–diffusion equations J.I. Ramos Room I-320-D, E.T.S. Ingenieros Industriales, Universidad de Ma´laga, Plaza El Ejido, s/n 29013 Ma´laga, Spain
Abstract Linearly-implicit, approximate factorization methods are used to reduce the solution of multi-dimensional reaction–diffusion equations to a sequence of two-point boundaryvalue problems in ordinary differential equations which are solved analytically in a piecewise manner. It is shown that, depending on the allocation of the transient, diffusion and reaction terms to the second-order ordinary differential operators, several families of piecewise exponential methods can be obtained. When only diffusion is included, a quadratic analytical solution is obtained, whereas if transient and diffusion, diffusion and reaction, and transient, diffusion and reaction terms are included, piecewise exponential methods result. It is shown that exponential methods that only account for reaction and diffusion in the differential operator may result in ill-conditioned systems of algebraic equation, whereas, if diffusion, transient and reaction are included, ill-conditioning is avoided if the time step is sufficiently small. It is also shown that the most robust exponential method is that that accounts for transient and diffusion terms in the differential operator, and this method provides solutions in good accord with those
E-mail address: [email protected] 0096-3003/$ - see front matter 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2005.07.020
1610
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
of a linearly-implicit, approximate factorization technique which employs a three-point, second-order accurate discretization of the diffusion terms. It is also shown that a hybrid exponential method that solves the transient, reaction and diffusion terms at odd time levels and the transient and diffusion terms at even time levels also provides very accurate solutions. 2005 Elsevier Inc. All rights reserved. Keywords: RotheÕs methods; Approximate factorization; Linearly implicit techniques; Piecewise exponential solutions; Reaction–diffusion equations
1. Introduction Factorization methods are techniques that approximate a multi-dimensional operator by a sequence of one-dimensional ones. Since it is rarely possible to factorize exactly multi-dimensional operators, factorization techniques usually employ an iterative procedure to account for the factorization errors. If these errors are neglected, the resulting methods are referred to as approximate factorization techniques which transform a multi-dimensional nonlinear problem into a sequence of one-dimensional ones whose solution requires iterations. Iterations may be avoided by time linearization whereby the time variable in the original multi-dimensional operator is first discretized and the resulting nonlinear elliptic partial differential equation is then linearized with respect to the previous time level. Upon using approximate factorization into the resulting linear elliptic partial differential equation, one obtains a sequence of one-dimensional, linear, two-point boundary-value problems in each coordinate direction that may be discretized and written in terms of either the original dependent variables [1,2] or delta form [3,4]. These linearly-implicit, approximate factorization methods have been used with great success in the numerical solution of the Navier–Stokes and Euler equations [1–4], and nonlinear reaction–diffusion equations in one [5] and several spatial dimensions [6,7], where the diffusion terms have usually been discretized by means of second-order accurate finite difference formulae, although the techniques have also proved to be robust with three-point compact discretizations of the diffusion terms [8–10]. Since linearly-implicit, approximate factorization methods provide linear, two-point boundary-value problems in each spatial direction for reaction– diffusion equations, these problems can also be analyzed by means of locally one-dimensional techniques such as the diffusive and exponentially-fitted methods developed by the author [11] who considered one-dimensional timedependent advection–reaction–diffusion equations that upon time discretization and linearization resulted in second-order ordinary differential equations with non-constant coefficients. However, by assuming that the coefficients
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
1611
and right-hand sides of these equations are piecewise constant and integrating the resulting equations analytically, the author developed C0 and C1 exponentially-fitted methods that yield tridiagonal matrices. In C0 methods, the coefficients of the differential equation are assumed to be constant in the interval xi1 < x < xi+1 and the solution is assumed to be continuous at xi1, xi and xi+1, whereas, in C1 techniques, the coefficients of the differential equation are assumed to be constant in the interval xi1 < x < xi and the solution is assumed to be continuous and first-order differentiable at xi1 and xi. The advantage of C0 and C1 methods is that they can be used to study timedependent evolution problems such as those governed by parabolic differential equations and that these methods provide piecewise analytical solutions. By way of contrast, linearly-implicit, approximate factorization techniques that discretize the diffusion terms by either second-order or three-point, compact finite difference formulae provide only nodal values and require interpolation to determine an approximate continuous (in space) solution. On the other hand, C0 and C1 methods involve exponential terms and are more costly than standard finite difference discretizations of the diffusion terms. Recently, the author has applied C0 exponentially-fitted methods to a onedimensional, nonlinear reaction–diffusion equation [12] and found that these techniques may provide erroneous results when a Crank–Nicolson discretization is employed. In fact, in that reference, it was found that the discontinuities in the initial conditions disappeared rather slowly. Although these results are in qualitative agreement with those found by other authors who used standard finite difference discretizations for the diffusion terms [13,14], C0 methods were found to be rather sensitive to the time step employed in the calculations when the ordinary differential operators include transient, diffusion and reaction terms. For one-dimensional, one-phase Stefan problems that are characterized by the presence of a moving front, C0 techniques were also found to be rather sensitive to both the spatial step size and the time step [15], and it has been claimed that the application of C0 methods to multi-dimensional systems of reaction– diffusion equations may be plagued by serious numerical difficulties associated with the allocation of the transient, diffusion and reaction terms to the onedimensional differential operators, the spatial step size, the time step and the nonlinearity of the reaction terms [12]. It must be pointed out that the linearly-implicit, approximate factorization, exponential methods considered here (see also [12,15]) are, in fact, linear RotheÕs methods [16,17] since they are based on the discretization of the time variable while keeping continuous the spatial ones. Linear and nonlinear RotheÕs techniques have been previously used to study reaction–diffusion processes [12], one-dimensional, one-phase Stefan problems [15], two-dimensional advection–reaction–diffusion equations with boundary element methods [18], one-dimensional convection–diffusion equations with the sinc–Galerkin
1612
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
method [19], the one-dimensional nonlinear Schro¨dinger equation [20], reaction–diffusion problems with internal boundaries [21], the one-dimensional viscous Burgers equation with the Galerkin procedure [22], and the two-dimensional convection–reaction–diffusion equation [23]. Most of these studies, however, have used either standard finite difference or finite element methods for the first- and second-order spatial derivatives or locally one-dimensional methods, and considered only the convection and diffusion terms in the differential operator so that their exponential methods are able to deal with high Reynolds or Pe´clet number flows [24]. Some of these references have considered the solution of the locally one-dimensional operator to derive finite difference expressions. By way of contrast, the approach followed here considers different differential operators, provides three-point finite difference formulae upon imposing continuity requirements at three different grid points, and may consider transient, advection, diffusion and reaction terms in the differential operators. It should also emphasized that the exponential methods presented here are completely different from those derived from the discretization of the spatial coordinates, i.e., methods of lines, that yield systems of first-order, ordinary differential equations in time, which upon linearization result in systems of first-order, linear ordinary differential equations that have exponential solutions in time [25], even though the evaluation of the matrix exponential is a very hard task. The methods presented here provide piecewise exponential solutions in space. In this paper, we first (in Section 2) develop a variety of linearly-implicit, approximate factorization, exponential methods for the numerical solution of multi-dimensional reaction–diffusion equations that include as a subset the second-order standard finite difference discretization of the diffusion terms [5], and techniques that account for transient and diffusion, transient, diffusion and reaction, and reaction and diffusion terms in the differential operators. We then analyze the linear (Fourier–von Neumann) linear stability of the finite difference equations that result upon using the transient and diffusion terms in the differential operator, and comment on the difficulties that one encounters when trying to analyze the linear stability of the exponential methods that include reaction and diffusion, and transient, diffusion and reaction terms in the differential operators, in Section 3. In that section, we also point out the difficulties that are posed by using a Crank–Nicolson discretization both in terms of the linear stability and the characteristics of exponential techniques, and propose a remedy to the latter. In Section 4, the linearly-implicit, exponential, approximate factorization methods presented here are applied to one- and two-dimensional systems of nonlinear reaction–diffusion equations, and their numerical results are compared with those of linearly-implicit methods that employ second-order accurate discretizations of the diffusion terms [5,6]. Finally, a section on conclusions summarizes the main findings of this paper.
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
1613
2. Formulation The linearly-implicit, approximate factorization, exponential methods considered in this paper can be applied to one- and multi-dimensional problems. Since these methods are based on approximate factorization techniques that factorize a multi-dimensional operator into a sequence of one-dimensional ones, for the sake of conciseness, we shall only present the formulation for multi-dimensional problems. We, thus, consider the following system of twodimensional reaction–diffusion equations: oU o2 U o2 U ¼ 2 þ 2 þ SðU; t; x; yÞ; ot ox oy
ð1Þ
where S is a function in C1(RN · R · R · R ! RN), U 2 RN, x 2 R, y 2 R, and t 2 R, and the diffusion coefficients in the x- and y-directions have been set equal to unity. The independent variables t (0 6 t < 1), and x (0 6 x 6 Lx) and y (0 6 x 6 Ly) are the time and spatial coordinates, respectively, and Lx and Ly denote the domainÕs dimensions in the x- and y-directions, respectively. Although the methods presented here can be easily generalized to the case that the diffusion coefficient is a nonlinear tensor, we have adopted Eq. (1) to illustrate them because of its simplicity. Eq. (1) can be transformed into a system of nonlinear (elliptic) partial differential equations at each time step by discretizing the time variable by means of a h-method while keeping continuous the (independent) spatial variables. These nonlinear equations can, in turn, be transformed into linear ones by linearizing the nonlinear terms at tn+1 with respect to the previous time level, and the resulting equation may be expressed as Unþ1 Un o2 Unþ1 o2 U n o2 Unþ1 n n ¼h þ ð1 hÞ þ S þ hT k þ h k ox2 ox2 oy 2 þ ð1 hÞ
o2 Un þ hJn ðUnþ1 Un Þ; oy 2
ð2Þ
oF n n ðU ; t ; x; yÞ; oU
ð3Þ
where Tn ¼
oF n n ðU ; t ; x; yÞ; ot
Jn ¼
where k = tn+1 tn is the time step, Un = U(tn, x, y), and Sn = S(Un, tn, x, y). The values h = 0 and h = 1 correspond to first-order accurate (in time), explicit and implicit, respectively, methods, while h = 0.5 corresponds to a second-order accurate (in time) method. In this paper, we shall be interested on implicit techniques, and especially on 0.5 6 h 6 1.
1614
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
Eq. (3) may also be written as DU o2 DU o2 Un o2 DU o2 Un n n ¼h þ þ S þ hT k þ h þ þ hJn DU; k ox2 oy 2 ox2 oy 2
ð4Þ
where DU = Un+1 Un. Eq. (4) corresponds to a delta formulation and a coupled system of linear elliptic equations at each time level. This two-dimensional system may be reduced to sequences of one-dimensional equations by means of the following approximate factorization technique: o2 o2 I kh 2 I kdhJn I kh 2 I khJn DU ¼ RHS; ð5Þ ox oy where I denotes the unit or identity matrix, d þ ¼ 1; 2 n oU o2 Un n n RHS ¼ k þ þ S þ khT ox2 oy 2
ð6Þ ð7Þ
and the O(k2) factorization errors have been disregarded. Eq. (5) can also be written as the following system of uncoupled, ordinary differential equations which have non-constant coefficients due to the appearance of the Jacobian matrix, J, o2 n Lx ðDU Þ I kh 2 I kdhJ DU ¼ RHS; ð8Þ ox o2 Ly ðDUÞ I kh 2 I khJn DU ¼ DU ; ð9Þ oy which represent linear systems of one-dimensional operators in the x- and ydirections, respectively, and correspond to a delta formulation. These equations can also be written as o2 Lx ðU Þ I kh 2 I kdhJn U ¼ RHS; ð10Þ ox o2 n Ly ðUÞ I kh 2 I khJ Unþ1 ¼ U ; ð11Þ oy where 2 n oU o2 Un n n n n RHS ¼ k ð1 hÞ þ þ S hJ U þ khT . ox2 oy 2
ð12Þ
Note that in Eqs. (8) and (9) and Eqs. (10) and (11) the spatial coordinates have been kept continuous while time has been discretized; therefore, these equations correspond to a RotheÕs method [16].
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
1615
Since Eqs. (8)–(11) have the same form, hereon, we shall only work with only one-dimensional operator, e.g., kh
o2 U þ ðI khdHÞU ¼ R; ox2
ð13Þ
where U, H and R can be readily obtained from Eqs. (8) and (10), or Eqs. (9) and (11) if x and d in Eq. (13) are replaced for y and , respectively.
3. Numerical methods In this section, we consider a variety of methods for the numerical solution of Eq. (13). These methods are presented in an integrated fashion so that the reader can appreciate the differences in how to obtain them. 3.1. Diffusion operator method In this method, Eq. (13) is first written as kh
o2 U ¼ P R ðI khdHÞU; ox2
ð14Þ
which, upon assuming that P is constant in xi1 6 x 6 xi+1 and imposing the following conditions: Uðxi1 Þ ¼ Ui1 ;
Uðxi Þ ¼ Ui ;
Uðxiþ1 Þ ¼ Uiþ1 ;
yields the following three-point finite difference equation: kh kh kh 2 Ui1 þ 1 þ 2 2 I khdHi Ui 2 Uiþ1 ¼ Ri ; h h h
ð15Þ
ð16Þ
where h = xi+1 xi = xi xi1, which is exactly the same equation as the one that can be derived by discretizing the second-order derivative in Eq. (13) by means of second-order accurate finite differences in an equally-spaced grid [5]. Since U is a vector, Eq. (16) can be written as a system of linear algebraic equations where the matrix has a block tridiagonal structure. Hereon, we shall refer to the method represented by Eq. (16) as D. 3.2. Exponential transient–diffusion operator method In this method, Eq. (13) is first written as kh
o2 U þ U ¼ P R þ khdHU; ox2
ð17Þ
1616
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
which only includes the transient and diffusion operators in the left-hand side and, for h 5 0, can be integrated analytically in xi1 6 x 6 xi+1 by assuming that P is constant in this interval, to yield U ¼ A coshðkðx xi ÞÞ þ B sinhðkðx xi ÞÞ þ W;
ð18Þ
where A and B are constant vectors, k2 ¼ kh1 and W = Pi. Upon using Eq. (15), Eq. (18) yields the following three-point finite difference method: Ui1 2ðcoshðkhÞI þ khdð1 coshðkhÞÞHi ÞUi þ Uiþ1 ¼ 2Ri ð1 coshðkhÞÞ;
ð19Þ
which can also be written as a block tridiagonal system of linear algebraic equations. Hereon, we shall refer to the method represented by Eq. (19) as ETD. 3.3. Exponential reaction–diffusion operator method In this method, Eq. (13) is first written as the following two equations for N = 2: kh
o2 u khdH 11 u ¼ P 1 R1 u þ khdH 12 v; ox2
ð20Þ
kh
o2 v khdH 22 v ¼ P 2 R2 v þ khdH 21 u; ox2
ð21Þ
which includes the diffusion and part of the reaction terms in the left-hand side, where U = (u, v)T and the superscript T denotes transpose. For h 5 0, Eq. (20) or Eq. (21) can be integrated analytically in xi1 6 x 6 xi+1 by assuming that P is constant in this interval. For example, the solution of Eq. (20) subject to Eq. (15) yields h2 h2 uiþ1 þ ui1 2 þ ð22Þ ui þ dh2 ðH 12 Þi vi ¼ ðR1 Þi kh kh for H11 = 0, uiþ1 þ ui1 2 cosðkhÞ þ þ2
1 ð1 cosðkhÞÞ ui khdH 11
ðH 12 Þi ðR1 Þi ð1 cosðkhÞÞvi ¼ 2 ðcosðkhÞ 1Þ ðH 11 Þi khdðH 11 Þi
for dH11 > 0, where k2 = dH11, and
ð23Þ
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
uiþ1 þ ui1 2 coshðkhÞ þ þ2
1617
1 ð1 coshðkhÞÞ ui khdH 11
ðH 12 Þi ðR1 Þi ð1 coshðkhÞÞvi ¼ 2 ðcoshðkhÞ 1Þ ðH 11 Þi khdðH 11 Þi
ð24Þ
for dH11 < 0, where k2 = dH11. Eqs. (22)–(24) and the ones corresponding to v can be written as a block tridiagonal system of linear algebraic equations. Hereon, we shall refer to the method represented by Eqs. (22)–(24) as ERD. 3.4. Exponential transient–reaction–diffusion operator method In this method, Eq. (13) is first written as the following two equations for N = 2: kh
o2 u þ ð1 khdH 11 Þu ¼ P 1 R1 þ khdH 12 v; ox2
ð25Þ
kh
o2 v þ ð1 khdH 22 Þv ¼ P 2 R2 þ khdH 21 u; ox2
ð26Þ
which includes the transient term, diffusion and part of the reaction terms in the left-hand side, where U = (u, v)T and the superscript T denotes transpose. For h 5 0, Eq. (25) or Eq. (26) can be integrated analytically in xi1 6 x 6 xi+1 by assuming that P is constant in this interval. For example, the solution of Eq. (25) subject to Eq. (15) yields uiþ1 þ ui1 2ui þ h2 dðH 12 Þi vi ¼
h2 ðR1 Þi kh
ð27Þ
for Q1 = 1 khdH11 = 0, uiþ1 þ ui1 2 cosðkhÞui 2ð1 cosðkhÞÞ ¼ 2
khdðH 12 Þi vi 1 khdðH 11 Þi
ðR1 Þi ðcosðkhÞ 1Þ 1 khdðH 11 Þi
ð28Þ
for Q1 > 0, where k2 ¼ dH 11 kh1 , and uiþ1 þ ui1 2 coshðkhÞui 2 ¼ 2
khdðH 12 Þi ð1 coshðkhÞÞvi 1 khdðH 11 Þi
ðR1 Þi ðcoshðkhÞ 1Þ 1 khdðH 11 Þi
for Q1 < 0, where k2 ¼ dH 11 þ kh1 .
ð29Þ
1618
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
Eqs. (27)–(29) and the ones corresponding to v can be written as a block tridiagonal system of linear algebraic equations. Hereon, we shall refer to the method represented by Eqs. (27)–(29) as ETRD.
4. Linear stability In this section, we examine the linear stability of the methods D and ETD presented in the previous section for one-dimensional scalar problems. We, thus, consider the following partial differential equation: ou o2 u ¼ þ au; ot ox2 where a is a constant. Eq. (16) can be written taking into account Eq. (19) as kh nþ1 kh kh nþ1 2 ui1 þ 1 þ 2 2 kha uinþ1 2 uiþ1 h h h k 2k ¼ ð1 hÞ 2 ðuiþ1 þ ui1 Þ þ 1 þ ð1 hÞ ak 2 uni . h h
ð30Þ
ð31Þ
Using uni ¼ /n expðIiKhÞ where I2 = 1 in Eq. (31), one can readily obtain the following amplification factor: 2 4k 1 þ ð1 hÞ ka sin l 2 /nþ1 h ; ð32Þ ¼ 4kh /n 1 kha þ h2 sin2 l where l ¼ Kh and K is the wave number. 2 A sufficient condition for linear stability is that a 6 0 for h = 0.5 and 1. The discretization of the time variable in Eq. (30) by means of a h method yields the following ordinary differential equation: 2 n o2 unþ1 ou nþ1 n n kh þ ð1 khaÞu ¼ u þ kð1 hÞ þ au . ð33Þ ox2 ox2 If we apply to this equation the exponential techniques presented in the previous section, we may have to deal with the second-order spatial derivative in the right-hand side. For example, the use of ETD yields nþ1 nþ1 ui1 2ðcoshðkhÞ þ khað1 coshðkhÞÞunþ1 þ uiþ1 i 2 n o ui n n ¼ 2 kð1 hÞ þ au þ u i i ð1 coshðkhÞÞ; ox2
ð34Þ
which yields a tridiagonal system of linear algebraic equations, where k2 ¼ kh1 . A sufficient condition for this system to be diagonally dominant is a 6 0.
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
1619
Eq. (34) poses a problem when trying to analyze its linear stability because of the presence of the second-order spatial derivative, unless h = 1. Note that this derivative is evaluated at the previous time level and can be determined at t = 0 from the initial conditions, and, at subsequent times, from the piecewise analytical solution determined in Section 2 as follows. The analytical solution used in ETD can be written as (cf. Eq. (18)) unþ1 ¼ Ai coshðkðx xi ÞÞ þ Bi sinhðkðx xi ÞÞ þ P i ;
xi1 6 x 6 xiþ1 ; ð35Þ
where 2 n ou n P ¼ u þ kð1 hÞ þ au þ khaunþ1 ox2 n
ð36Þ
and Ai and Bi can be determined from Eqs. (15) and (35) as Ai ¼ uinþ1 P i ;
Bi ¼
nþ1 uiþ1 unþ1 iþ1 sinhðkhÞ
ð37Þ
2 nþ1
and, therefore, o oxu 2 ¼ k2 ðAi coshðkðx xi ÞÞ þ Bi sinhðkðx xi ÞÞÞ can be easily determined. Note, however, that Ai and Bi depend on the solution at the previous time level through Eqs. (36) and (37). Moreover, in ERD and ETRD, k is also a function of the previous time level and this complicates the stability analysis of these techniques. In addition, it is easily shown that ETD, ERD and ETRD result in matrices that are diagonally dominant if a is negative (cf. Eqs. (19), (24) and (29), respectively.) In view of these comments, in the next subsection, we analyze the linear stability of TD for h = 1. 4.1. Stability of the transient–diffusion method For h = 1, the use of uni ¼ /n expðIiKhÞ, where I2 = 1, in Eq. (34) yields the following amplification factor: /nþ1 1 coshðkhÞ ¼ cosð2lÞ ðcoshðkhÞ þ khað1 coshðkhÞÞ /n ¼
sinh2 b ; sin2 l þ ð1 khaÞsinh2 b
ð38Þ
where b ¼ kh2 , and this factor has an absolute value equal to or less than one if a 6 0, therefore, ETD is linearly stable according to the Fourier–von Neumann analysis presented here. It must be noted that, for kh 1, a TaylorÕs series expansion of Eq. (34) yields Eq. (31) up to O((kh)4), for h = 1.
1620
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
5. Results In this section, we present some sample results for one- and two-dimensional nonlinear reaction–diffusion equations obtained with the methods presented in this paper. 5.1. One-dimensional reaction–diffusion equations The one-dimensional reaction–diffusion system considered in this subsection can be written as ou o2 u ¼ uv; ot ox2 ov o2 v 1 ¼ þ uv v; ot ox2 2
ð39Þ ð40Þ
subject to u(0, x) = 1, v(0, x) = exp(x2), u(t, 200) = u(t, 200) = 1 and v(t, 200) = v(t, 200) = 0. This problem was solved in an equally-spaced grid consisting of NI = 800 grid points, h = 0.5 and several time steps, and some sample results are presented in Figs. 1 and 2. Fig. 1 corresponds to k = 0.05 and shows the solution obtained with D. The results presented in Fig. 1 indicate that the initial peak in v becomes two travelling waves in time that generate a widening valley in u. The errors of ETD relative to D with k = 0.01, i.e., duðkÞ ¼ uD ðk ¼ 0:01Þ uETD ðkÞ and dvðkÞ ¼ vD ðk ¼ 0:01Þ vETD ðkÞ, are shown in Figs. 2 and 3 as functions of space and time, and the time step. Figs. 2 and 3 indicate that the relative errors of both u and v increase as time is increased and are located at the inflection points of the u profile and at the relative maxima of v. The largest errors du and dv of ETD at t = 100 are less than 1 · 102 and 4 · 103, 4 · 103 and 2 · 103, 15 · 103 and 5 · 103, and 0.02 and 0.01, respectively, for k = 0.05, 0.025, 0.1 and 0.01, respectively.
Fig. 1. u (left) and v (right) as functions of x and t obtained with D, k = 0.05, NI = 800.
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
1621
Fig. 2. EU ðkÞ ¼ uD ðk ¼ 0:01Þ uETD ðkÞ (left) and EV ðkÞ ¼ vD ðk ¼ 0:01Þ vETD ðkÞ (right) as functions of x and t obtained with NI = 800. (Top: k = 0.05; bottom: k = 0.025).
On the other hand, the errors duðkÞ ¼ uD ðkÞ uETD ðkÞ and dvðkÞ ¼ vD ðkÞ vETD ðkÞ of ETD were found to decrease as the time step is decreased and, for k = 0.01 and 0.1, are, at most, 0.02 and 0.01, and 0.1 and 0.15, respectively, for u and v, respectively. In Tables 1 and 2, we show the relative errors of ETD with respect to D, as well as the relative errors of the later with respect to itself as functions of the time step. Although not shown here, the errors uD ðk ¼ 0:01Þ uD ðkÞ and vD ðk ¼ 0:01Þ vD ðkÞ exhibit the same trends as those presented in Figs. 2 and 3. The method ETRD provided nearly identical results to those shown in Figs. 2 and 3 and, for this reason, they are not presented here. Note that for jkhd(H11)ij 1, the denominators that appear in Eqs. (28) and (29) can be approximated by one, and Eq. (29) tends to Eq. (19). ERD was found to be not as robust as ETD for the grid size and time steps mentioned above. In fact, for these time steps, this method was found to incur in errors when evaluating the sine and cosine terms that appear in their finite difference equations (cf. Eq. (23)). As one can easily check, ERD may not result in diagonally dominant matrices for k = 0.05, 0.025 and 0.1, and, in some of the calculations presented here, this method was found to provide ill-conditioned matrices, especially for NI = 800. This ill-conditioning was still observed
1622
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
Fig. 3. EU ðkÞ ¼ uD ðk ¼ 0:01Þ uETD ðkÞ (left) and EV ðkÞ ¼ vD ðk ¼ 0:01Þ vETD ðkÞ (right) as functions of x and t obtained with NI = 800. (Top: k = 0.1; bottom: k = 0.01).
Table 1 Maximum errors in u of ETD with respect to D, i.e., E ¼ uD uETD , and e ¼ uD ðk ¼ 0:01Þ uD ðkÞ as functions of k Error
k = 0.05
k = 0.025
k = 0.1
k = 0.01
E e
0.15 0.15
0.06 0.06
0.3 0.3
0.02 0
Table 2 Maximum errors in v of ETD with respect to D, i.e., E ¼ vD vETD , and e ¼ vD ðk ¼ 0:01Þ vD ðkÞ as functions of k Error
k = 0.05
k = 0.025
k = 0.1
k = 0.01
E e
0.04 0.05
0.02 0.02
0.15 0.1
0.01 0
in ERD even for k = 0.01, 0.005, 0.0025 and 0.001 and is also due to the fact that, as indicated in Fig. 1, H11 = v is nil wherever v = 0 and this implies that the solution behaves as a quadratic approximation for H11 wherever v = 0 and as a trigonometric function wherever dH11 > 0.
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
1623
As will be shown below for two-dimensional reaction–diffusion equations, ETRD behaves very well and ERD is not robust due to the reasons exposed above, whereas a combination of ETRD and TD may yield good results. 5.2. Two-dimensional reaction–diffusion equations The two-dimensional reaction–diffusion system considered in this subsection can be written as ou o2 u o2 u ¼ þ u2 v; ot ox2 oy 2
ð41Þ
ov o2 v o2 v 1 ¼ þ þ u2 v v; ot ox2 oy 2 2
ð42Þ
subject to u(0, x, y) = 1, v(0, x, y) = exp((x2 + y2)), u(t, 20, x) = u(t, 20, y) = u(t,x,20) = u(t,x,20) = 1 and v(t,20,x) = v(t,20,y) = v(t,x,20) = v(t,x,20) = 0. This problem was solved in an equally-spaced grid consisting of NI = NJ = 102 grid points in the x- and y-coordinates, respectively, h = 0.5 and several time steps, and some sample results are presented in Figs. 4 and 5. Fig. 4 shows u obtained with D as a function of space at different times
Fig. 4. u as a function of x and y obtained with D, k = 0.02, NI = NJ = 102. (From left to right, top to bottom: t = 4, 6, 8, 10, 12, 14, 16, 18, 20).
1624
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
Fig. 5. v as a function of x and y obtained with D, k = 0.02, NI = NJ = 102. (From left to right, top to bottom: t = 4, 6, 8, 10, 12, 14, 16, 18, 20).
and indicates that u exhibits a valley in the middle of the computational domain which is generated by the reaction term and the initial condition for v. Both the depth and the width of this valley increase with time until about t = 20 when u reaches an almost stationary value. Fig. 5 shows that, owing to the initial conditions, v develops initially a peak in the middle of the domain. The height of this peak decreases, whereas its width increases as time increases. At about t = 8, the peak flattens and a valley is formed. The depth of this valley increases as time increases until about t = 20 when v reaches an almost stationary profile characterized by high values near the domain boundaries. In order to assess the accuracy of the linearly-implicit, approximate factorization, exponential methods presented in this paper, we have performed comparisons with D and defined the local errors as EU M ¼ uD ðx; y; tÞ uM ðx; y; tÞ and EV M ¼ vD ðx; y; tÞ vM ðx; y; tÞ, where the subscripts D and M denote the solutions obtained with D and M, respectively, and M = ETRD, ETD and RD. Figs. 6 and 7 show the errors of ETRD and ETD, respectively, as functions of space at different times for k = 0.02 and h = 0.5. The errors of ETRD are, at most, 1 · 103 at t = 4, 6, 8, 10 and 12, and, at most, 2 · 103 at t = 14, 16, 18 and 20, and occur where either juj or its gradient are largest as indicated in Fig. 6. The errors of ETD illustrated in Fig. 7 exhibit the same trends as those of ETRD and are, at most, 1 · 103 at t = 4, 6 and 8, at most, 2 · 103 at t = 10, 12 and 14, and, at most, 4 · 103 at t = 16, 18 and 20, i.e., they are larger than those of ETRD.
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
1625
Fig. 6. Errors EU ETRD ¼ du ¼ uD uETRD , where D and ETRD denote the solutions obtained with D and ETRD, respectively, k = 0.02, NI = NJ = 102. (From left to right, top to bottom: t = 4, 6, 8, 10, 12, 14, 16, 18, 20).
Although not shown here, it was observed that ERD exhibited similar problems to those found for one-dimensional systems of reaction–diffusion equations, i.e., it provided ill-conditioned matrices. On the other hand, it was also found that a combination of ETRD and ETD which uses ETRD and ETD in the x- and y-operators, respectively, at odd time levels, and in the y- and x-operators, respectively, at even time levels, respectively, provides the errors shown in Fig. 8. Hereon, we shall refer to the hybrid method just described as ETRD–TD. Fig. 8 shows that the errors of ETRD–TD are, at most, 1 · 103 at t = 4, 6, 8 and 10, at most 2 · 103 at t = 12, 14 and 16, and, at most, 4 · 103 at t = 18 and 20. Therefore, the errors of ETRD–TD are higher and smaller, respectively, than those of ETRD and ETD, respectively. Figs. 9–11 illustrate the errors in v of ETRD, ETD and ETRD–TD, respectively, as functions of space at different times for k = 0.02 and h = 0.5. The errors EV ETRD are, at most, 1 · 103 at t = 4, 14, 16, 18 and 20, and 4 · 104 at t = 6, 8, 10 and 12. The errors EV ETD and EV ETRD–TD are, at most, 1 · 103 at t = 4, 10, 12, 14, 16, 18 and 20, and 4 · 104 at t = 6 and 8. The effects of the time step on the accuracy of ETRD, ETD and ETRD–TD relative to D are presented in Tables 3 and 4 for u and v, respectively. Tables 3 and 4 indicate that the relative accuracy increases as k is decreased, ETRD is
1626
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
Fig. 7. Errors EU ETD ¼ du ¼ uD uETD , where D and ETD denote the solutions obtained with D and ETD, respectively, k = 0.02, NI = NJ = 102. (From left to right, top to bottom: t = 4, 6, 8, 10, 12, 14, 16, 18, 20).
Fig. 8. Errors EU ETRDTD ¼ du ¼ uD uETRDTD , where D and ETRD–TD denote the solutions obtained with D and ETRD–TD, respectively, k = 0.02, NI = NJ = 102. (From left to right, top to bottom: t = 4, 6, 8, 10, 12, 14, 16, 18, 20).
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
1627
Fig. 9. Errors EV ETRD ¼ dv ¼ vD vETRD , where D and ETRD denote the solutions obtained with D and ETRD, respectively, k = 0.02, NI = NJ = 102. (From left to right, top to bottom: t = 4, 6, 8, 10, 12, 14, 16, 18, 20).
more accurate than ETRD–TD which, in turn, is more accurate than ETD. This means that the accuracy of linearly-implicit, approximate factorization, exponential methods that include the reaction, transient and diffusion terms when determining the piecewise analytical solution of the one-dimensional operators is higher than that of exponential techniques that only include the transient and diffusion terms in these operators. On the other hand, one-dimensional operators that only include diffusion and reaction may result in ill-conditioned matrices. Figs. 6–11 show that the linearly-implicit, approximate factorization, exponential methods presented in this paper preserve the symmetries of the original problem (cf. Eqs. (41) and (42)). Although the results presented in this paper indicate that linearly-implicit, exponential methods have an accuracy comparable to those of linearly-implicit techniques that use standard second-order accurate finite difference formulae and, therefore, provide piecewise analytical solution and there is no need for interpolation, they may be subject to overflow for large spatial step sizes and/or small time steps because the solution of the one-dimensional operators depends in an exponential manner on the ratio phffiffik. Furthermore, when the initial conditions are not C2 and h = 0.5, linearly-implicit, exponential are not
1628
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
Fig. 10. Errors EV ETD ¼ dv ¼ vD vETD , where D and ETD denote the solutions obtained with D and ETD, respectively, k = 0.02, NI = NJ = 102. (From left to right, top to bottom: t = 4, 6, 8, 10, 12, 14, 16, 18, 20). 2
2
strictly applicable because they do require the values of ooxU2 and/or ooyU2 of the initial conditions, at the first time level (cf. Eqs. (7) and (12)). When this occurs, these second derivatives may be evaluated approximately by means of, for example, second-order accurate finite difference formulae which correspond to a quadratic (rather than exponential) approximation.
6. Conclusions Linearly-implicit, approximate factorization, exponential methods for the numerical solution of multi-dimensional reaction–diffusion equations have been developed. These techniques are based on the discretization of the time derivatives and the time linearization of the nonlinear terms, and yield linear elliptic partial differential equations at each time level. Upon factorizing the multi-dimensional elliptic operator into a sequence of one-dimensional operators while keeping the spatial coordinates continuous, two-point, linear boundary-value problems are obtained. These boundary-value problems are characterized by non-constant coefficients; however, they can be linearized (in space) in a piecewise fashion to yield linear, constant coefficients, second-order ordinary differential equations that can be solved in closed form to obtain
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
1629
Fig. 11. Errors EV ETRD–TD ¼ dv ¼ vD vETRD–TD , where D and ETRD–TD denote the solutions obtained with D and ETRD–TD, respectively, k = 0.02, NI = NJ = 102. (From left to right, top to bottom: t = 4, 6, 8, 10, 12, 14, 16, 18, 20).
piecewise analytical solutions which, upon imposing, continuity conditions at three different grid points yield block tridiagonal systems of linear algebraic equations. Depending on the allocation of the transient, diffusion and reaction terms to the second-order ordinary differential equations, several exponential techniques have been developed. If the differential operator only includes diffusion, these methods provide piecewise quadratic solutions and their finite difference equations are identical to those that result from a second-order accurate finite difference discretization in space. If transient and diffusion, transient, diffusion and reaction, and diffusion and reaction terms are included in the differential operators, these techniques provide exponential finite difference equations that depend on the diffusion, diffusion and reaction, and the reaction times, respectively. It has been show that, for the one- and two-dimensional reaction–diffusion problems considered in this paper, exponential methods that account for transient and diffusion, and transient, diffusion and reaction terms in the differential operators provide solutions of nearly the same accuracy as that of second-order (in space) accurate, linearly-implicit techniques. However, these techniques may be subject to overflow when the inverse of the Fourier number is very large, due to the exponential character of the solution. On the other
1630
Error
t=4
t=6
t=8
t = 10
t = 12
t = 14
t = 16
t = 18
t = 20
EU1(k = 0.04) EU1(k = 0.02) EU1(k = 0.01)
3
2 · 10 1 · 103 4 · 104
3
2 · 10 1 · 103 4 · 104
3
2 · 10 1 · 103 4 · 104
3
2 · 10 1 · 103 4 · 104
3
2 · 10 1 · 103 1 · 103
3
2 · 10 2 · 103 1 · 103
3
4 · 10 2 · 103 1 · 103
3
4 · 10 2 · 103 1 · 103
4 · 103 2 · 103 1 · 103
EU2(k = 0.04) EU2(k = 0.02) EU2(k = 0.01)
2 · 103 1 · 103 4 · 104
2 · 103 1 · 103 4 · 104
2 · 103 1 · 103 1 · 103
4 · 103 2 · 103 1 · 103
4 · 103 2 · 103 1 · 103
4 · 103 2 · 103 2 · 103
5 · 103 4 · 103 2 · 103
5 · 103 4 · 103 2 · 103
1 · 102 4 · 103 2 · 103
EU3(k = 0.04) EU3(k = 0.02) EU3(k = 0.01)
2 · 103 1 · 103 4 · 104
2 · 103 1 · 103 4 · 104
2 · 103 1 · 103 4 · 104
2 · 103 1 · 103 1 · 103
4 · 103 2 · 103 1 · 103
4 · 103 2 · 103 1 · 103
4 · 103 2 · 103 2 · 103
4 · 103 4 · 103 2 · 103
4 · 103 4 · 103 2 · 103
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
Table 3 Errors in u of ETRD, ETD and ETRD–TD with respect to D, i.e., EU M ¼ uD uM , where M = 1, 2 and 3 correspond to ETRD, ETD and ETRD–TD, respectively
Error
t=4
t=6
t=8
t = 10
t = 12
t = 14
t = 16
t = 18
t = 20
EV1(k = 0.04) EV1(k = 0.02) EV1(k = 0.01)
3
2 · 10 1 · 103 4 · 104
3
1 · 10 4 · 104 4 · 104
3
1 · 10 4 · 104 2 · 104
3
1 · 10 4 · 104 2 · 104
3
1 · 10 4 · 104 2 · 104
3
1 · 10 1 · 103 1 · 103
3
1 · 10 1 · 103 1 · 103
3
1 · 10 1 · 103 1 · 103
1 · 103 1 · 103 1 · 103
EV2(k = 0.04) EV2(k = 0.02) EV2(k = 0.01)
1 · 103 1 · 103 4 · 104
1 · 103 4 · 104 2 · 104
1 · 103 4 · 104 2 · 104
1 · 103 1 · 103 4 · 104
2 · 103 1 · 103 4 · 104
2 · 103 1 · 103 1 · 103
2 · 103 1 · 103 1 · 103
2 · 103 1 · 103 1 · 103
2 · 102 1 · 103 1 · 103
EV3(k = 0.04) EV3(k = 0.02) EV3(k = 0.01)
2 · 103 1 · 103 4 · 104
1 · 103 4 · 104 2 · 104
1 · 103 4 · 104 2 · 104
1 · 103 4 · 104 4 · 104
1 · 103 1 · 103 4 · 104
1 · 103 1 · 103 1 · 103
2 · 103 1 · 103 1 · 103
2 · 103 1 · 103 1 · 103
2 · 103 1 · 103 1 · 103
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
Table 4 Errors in v of ETRD, ETD and ETRD–TD with respect to D, i.e., EV M ¼ vD vM , where M = 1, 2 and 3 correspond to ETRD, ETD and ETRD–TD, respectively
1631
1632
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
hand, exponential methods that only account for the diffusion and reaction terms were found to provide ill-conditioned matrices. The exponential methods presented in this paper can be easily applied to fully implicit methods, i.e., h = 1, when these techniques employ a primitive variable form. However, if either a primitive variable formulation for h = 0.5 or a delta form for 0 < h 6 1 is employed, these exponential techniques require that the solution be continuous with continuous second-order (spatial) derivatives because of the appearance of these derivatives in the right-hand side of the ordinary differential operators. Otherwise, these techniques are not applicable in their exponential formulation, unless the second-order derivatives at the previous time level are discretized by means of second-order accurate finite difference formulae.
Acknowledgement The research reported in this paper was partially supported by Project BFM2001-1902 from the Ministerio de Ciencia y Tecnologı´a of Spain and fondos FEDER.
References [1] W.R. Briley, H. McDonald, On the structure and use of linearized block implicit methods, J. Comput. Phys. 34 (1980) 54–73. [2] W.R. Briley, H. McDonald, An overview and generalization of implicit Navier–Stokes algorithms and approximate factorization, Comput. Fluids 30 (2001) 807–828. [3] R.M. Beam, R.F. Warming, An implicit finite-difference algorithm for hyperbolic systems in conservation-law form, J. Comput. Phys. 22 (1976) 87–110. [4] R.M. Beam, R.F. Warming, An implicit factored scheme for the compressible Navier–Stokes equations, AIAA J. 16 (1978) 393–402. [5] J.I. Ramos, Linearization methods for reaction–diffusion equations: 1-D problems, Appl. Math. Comput. 88 (1997) 199–224. [6] J.I. Ramos, Linearization methods for reaction–diffusion equations: multidimensional problems, Appl. Math. Comput. 88 (1997) 225–254. [7] J.I. Ramos, Linearized factorization techniques for multidimensional reaction–diffusion equations, Appl. Math. Comput. 100 (1999) 201–222. [8] J.I. Ramos, Implicit, compact, linearized h-methods with factorization for multidimensional reaction–diffusion equations, Appl. Math. Comput. 94 (1998) 17–43. [9] Y. Gu, W. Liao, J. Zhu, An efficient high-order algorithm for solving systems of 3-D reaction– diffusion equations, J. Comput. Appl. Math. 155 (2003) 1–17. [10] W. Liao, J. Zhu, A.Q.M. Khaliq, An efficient high-order algorithm for solving systems of reaction–diffusion equations, Numer. Meth. Partial Diff. Equat. 18 (2002) 340–354. [11] J.I. Ramos, On diffusive methods and exponentially-fitted techniques, Appl. Math. Comput. 103 (1999) 69–96. [12] J.I. Ramos, Exponential methods for one-dimensional reaction–diffusion equations, Appl. Math. Comput., in press, doi:10.1016/j.amc.2004.12.003.
J.I. Ramos / Appl. Math. Comput. 174 (2006) 1609–1633
1633
[13] D. Britz, O. Østerby, J. Strutwolf, Damping of Crank–Nicolson error oscillations, Comput. Biol. Chem. 27 (2003) 253–263. [14] O. Østerby, Five ways of reducing the Crank–Nicolson oscillations, BIT 43 (2003) 811–822. [15] J.I. Ramos, Exponential numerical methods for one-dimensional, one-phase Stefan problems, Arch. Appl. Mech., in press. [16] E. Rothe, Zweidimensionale parabolische Randwertaufgaben als eindimensionaler Randwertaufgaben, Math. Ann. 102 (1930) 652–670. [17] A. Ashyralyev, P.E. Sobolevskii, Well-posedness of Parabolic Difference Equations, Birkha¨user Verlag, Basel, Switzerland, 1994. [18] N. Samec, L. Sˇkerget, Integral formulation of a diffusive–convective transport equation for racting flows, Eng. Anal. Bound. Elem. 28 (2004) 1055–1060. [19] J.L. Mueller, T.S. Shores, A new sinc–Galerkin method for convection–diffusion equations with mixed boundary conditions, Comput. Math. Appl. 47 (2004) 803–822. [20] F. Schmidt, D. Yevick, Discrete boundary conditions for Schro¨dinger–type equations, J. Comput. Phys. 134 (1997) 96–107. [21] J. Lang, High–resolution self–adaptive computations on chemical reaction–diffusion problems with internal boundaries, Chem. Eng. Sci. 51 (1996) 1055–1070. ¨ zdesß, A numerical solution of BurgersÕ equation, Appl. Math. Comput. 156 [22] E.N. Aksan, A. O (2004) 395–402. [23] T.W.H. Sheu, S.K. Wang, R.K. Lin, An implicit scheme for solving the convection–diffusion– reaction equation in two dimensions, J. Comput. Phys. 164 (2000) 123–142. [24] K.W. Morton, Numerical Solution of Convection–Diffusion Problems, Chapman & Hall, London, UK, 1996. [25] P.J. van der Houwen, B.P. Sommeijer, Approximate factorization for time-dependent partial differential equations, J. Comput. Appl. Math. 128 (2001) 447–466.
