Iterative and non-iterative, full and approximate factorization methods for multidimensional reaction–diffusion equations

Applied Mathematics and Computation 174 (2006) 1586–1608

www.elsevier.com/locate/amc

Iterative and non-iterative, full and approximate factorization methods for multidimensional 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 A variety of second-order accurate, in both space and time, full and approximate factorization methods for the numerical solution of two-dimensional reaction–diffusion equations is presented. These methods may use time linearization and yield linearly implicit techniques and one-dimensional operators in each direction. It is shown that, if the factorization errors are neglected, linearly implicit approximate factorization methods provide uncoupled equations, whereas, if these errors are considered, the equations are coupled and must be solved iteratively. It is also shown that the allocation of the reaction and diffusion terms to the one-dimensional operators plays a paramount role in determining the accuracy of approximate factorization methods and preserving the symmetry of the original differential problem. Iterative, full and approximate factorization methods that do require iterations are also presented, and, for the problem considered here, these methods are shown to converge in about two iterations and provide solutions in agreement with those obtained with linearly implicit full and approximate factorization techniques.
2005 Elsevier Inc. All rights reserved.

E-mail address: [email protected] 0096-3003/$ - see front matter
2005 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2005.07.007

J.I. Ramos / Appl. Math. Comput. 174 (2006) 1586–1608

1587

Keywords: Factorization methods; Linearly implicit techniques; Iterative methods; Reaction–diffusion equations

1. Introduction Factorization methods are a branch of numerical operator-splitting techniques whereby a multidimensional operator is factorized as a sequence of lower dimensional ones, and can be broadly classified as full or approximate factorization algorithms. Full factorization methods account for the factorization errors incurred in the factorization of multidimensional operators into a sequence of one-dimensional ones which are much easier to solve than the original operator, whereas approximate factorization techniques disregard the factorization errors. The first application of approximate factorization to the numerical solution of time-dependent or evolutionary partial differential equations can be traced back to the papers of Peaceman and Rachford [1] and Douglas [2] which resulted in the development of alternating-direction implicit (ADI) techniques. More explicitly, approximate factorization was formulated by Beam and Warming [3,4] and Briley and McDonald [5] in delta and original variable form, respectively. It must be pointed out that the approximate factorization methods of Beam and Warming [3,4] and Briley and McDonald [5] are linearly implicit (approximate factorization) techniques because the nonlinear terms are linearized with respect to the previous time level and are, thus, equivalent to the first iteration of an iterative Newton–Raphson–Kantorovich (NRK) procedure, and suffer from both the linearization of the nonlinear terms and the neglect of the factorization errors. If, in approximate factorization, the linearization is performed with respect to the previous iterated values as in the NRK method, one obtains a quasilinearization technique which still does not account for the factorization errors, although these methods are iterative due to the presence of nonlinear terms. It must also be emphasized that the linearly implicit approximate factorization methods of Beam and Warming [3,4] and Briley and McDonald [5] are based on the discretization of the time variable and the (time) linearization of the resulting nonlinear elliptic equations at each time level. After the approximate factorization, these techniques provide two-point linear boundary-value problems with continuous spatial variables; if the spatial variables are discretized, then one obtains systems of linear algebraic equations. However, if the spatial variables in the original evolutionary partial differential equation are first discretized one obtains a (nonlinear) system of first-order ordinary differential equations which, upon linearization, becomes a system of linear ordinary differential equations of the form dy ¼ fðy; tÞ ¼ Ay þ gðtÞ where y and t denote dt

1588

J.I. Ramos / Appl. Math. Comput. 174 (2006) 1586–1608

the nodal values of the dependent variables, t is time, and A is a constant matrix. By decomposing the matrix A into the sum of matrices, one can use approximate factorization methods for the solution of the resulting first-order ordinary differential equations. Therefore, approximate factorization techniques can be used to solve linear and nonlinear ordinary or time-dependent partial differential equations, and result in methods such as the linear multistep approach of Warming and Beam [6], linearly implicit integration methods like Rosenbrock techniques, linearization of nonlinear techniques, iterative solution of linear systems, etc. [7,8]. Although approximate factorization methods have been applied to a variety of problems in fluid dynamics, heat transfer, combustion, biology, ecology, etc., hereon we shall be concerned with the application of these techniques to two-dimensional reaction–diffusion equations. In these applications, both iterative and non-iterative linearly implicit approximate factorization techniques have been employed, and most of them suffer from the fact that they do not account for the factorization errors. Thus, for example, the author developed families of linearly implicit approximate factorization methods for two-dimensional reaction–diffusion equations based on the different allocation of the diagonal and off-diagonal components of the Jacobian matrices of the source terms [9], one-dimensional problems [10], and two-dimensional problems [11] that employ second-order accurate discretizations of the (secondorder derivative) diffusion terms. When three-point compact finite difference formulae are employed to discretize the diffusion terms, the author had to introduce the first- and second-order spatial derivatives of the dependent variables as additional unknowns [12], and this resulted in a large increase in the bandwidth of the coefficient matrix and storage. This can, however, be avoided by appropriate rearrangement of the factorization operators [13,14]. However, in all these references, i.e., [10,11,13,14], the authors neglected the factorization errors that result upon factorizing the multidimensional operator into a sequence of one-dimensional ones and employed three-point second- or fourthorder compact finite difference formulae for the discretization of the diffusion terms. As stated previously, linearly-implicit factorization methods are based on the time discretization and linearization, and factorization of operators whereby the spatial variables or coordinates are kept continuous; therefore, these techniques are also linear RotheÕs methods which require the solution of linear second-order ordinary differential equations at each time level. The numerical solution of these equations can be achieved by means of either finite difference and finite element methods or exponential techniques [15,16] based on the piecewise linearization of these ordinary differential equations and their (piecewise) analytical integration. It is known that approximate factorization methods introduce factorization errors and may need the use of smaller time steps than those required to follow

J.I. Ramos / Appl. Math. Comput. 174 (2006) 1586–1608

1589

the time evolution of unsteady phenomena, especially in viscous flow computations when they are used as linearly implicit techniques, and/or a decrease in the convergence speed when they are used in an iterative fashion [17]. This has led to the development of iterative modified approximate factorization or iterative full factorization techniques that account for the full nonlinear terms and may account for the factorization errors [17–19]. Briley and McDonald [20] do recommend the use of noniterative approximate factorization methods for steady problems and the use of iterative methods for unsteady ones. Iterative methods that account for the factorization errors have been developed by Steinthorsson and Shih [21] for linear hyperbolic equations and by the author for nonlinear parabolic partial differential equations [9,12]. These authors eliminated the higher order derivatives that appear in the factorization and developed predictor–corrector strategies to solve iteratively the factorized equations. They also included studies on the stability of the full factorization and the convergence of their iterative techniques. In this paper, we shall refer to factorization techniques as full ones when they do account for the factorization errors, and as approximate factorization methods, otherwise. It must stated at the outset that the name full factorization derives from the fact that these methods account for the factorization methods, since the original operator can rarely be factorized in an exact manner. It should also be emphasized that full factorization methods do not provide exact solutions because the discretization of the time and space variables introduces truncation errors. As will be shown in this paper, the approximate factorization of a linearly implicit method for multidimensional reaction–diffusion equations is not unique. In fact, for a second-order temporal accuracy and two-dimensional reaction–diffusion equations, it can be easily shown (Section 3) that there is a one-parameter family of approximate factorization techniques, depending on how the reaction terms are allocated to the one-dimensional operators. This also occurs in full factorization methods, but, since these techniques account in an iterative manner for the factorization errors, they are expected to yield almost identical results. This, however, may not be the case in approximate factorization methods. The objective of this paper is severalfold. First, in Section 3, linearly implicit full factorization methods are presented based on different factorizations of the diffusion and reaction terms, and the diffusion terms. Second, several linearly implicit approximate factorization techniques that account for either the diffusion and reaction terms or the diffusion in the one-dimensional operators are presented in Section 4. Third, the accuracy and symmetry preservation of both full and approximate factorization methods is assessed in terms of local errors in both space and time, as functions of the time step in Section 5. Finally, a summary of the main findings is presented in the conclusions. It should be stressed that throughout the paper we shall be concerned with second-order

1590

J.I. Ramos / Appl. Math. Comput. 174 (2006) 1586–1608

accurate time discretizations, even though the methods are formulated in a general fashion, and that only second-order finite difference discretizations of the diffusion terms in equally-spaced grids are considered.

2. Formulation Consider the following system of two-dimensional 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, where 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 and the result can be expressed as Unþ1
Un o2 Unþ1 o2 U n o2 Unþ1 o2 Un ¼h þ ð1
hÞ 2 þ h þ ð1
hÞ 2 2 2 k ox ox oy oy þ hSnþ1 þ ð1
hÞSn ; n+1

n

ð2Þ n

n

n

n

n

where k = t
t is the time step, U = U(t , x, y), and S = S(U , t , 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 secondorder accurate (in time) method. In this paper, we shall be interested on implicit techniques, and especially on 0.5 6 h 6 1. Instead of using a single parameter h to define the implicitness of the reaction and diffusion processes in Eq. (2), it is possible to use three parameters, 0 6 hx 6 1, 0 6 hy 6 1 and 0 6 hR 6 1 for the diffusion terms in the x- and y-directions and the reaction terms, respectively. However, the resulting Eq. (2) would only be O(k2) accurate for hx = hy = hR = 0.5. Eq. (2) is, in general, nonlinear due to the nonlinearity of S; however, if the nonlinear term Sn+1 is approximated by means of its Taylor polynomial of first degree around (Un, tn, x, y), Eq. (2) becomes the following linear elliptic partial differential equation at each time step

J.I. Ramos / Appl. Math. Comput. 174 (2006) 1586–1608

1591

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

ð3Þ

where Tn ¼

oF n n ðU ; t ; x; yÞ; ot

Jn ¼

oF n n ðU ; t ; x; yÞ. oU

ð4Þ

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

ð5Þ

where DU = Un+1
Un. Eq. (5) corresponds to a delta formulation, and both Eqs. (3) and (5) correspond to large systems of linear algebraic equations and couple all the dependent variables at (tn+1, x, y). This two-dimensional system may be reduced to sequences of one-dimensional equations by means of the factorization techniques presented in the next sections. It must be noted that the factorization of Eq. (3) or (5) is not unique, but, as will be shown below, some factorization methods may be more convenient than others when one accounts for the factorization errors.

3. Linearly implicit full factorization methods In this section, we present some factorization methods for Eq. (5). 3.1. First linearly implicit full factorization method The two-dimensional differential operator which appears in (5) may be written in the following full factorization form:

 o2 o2 I
kh 2 I
kdhJn ð6Þ I
kh 2 I
k
hJn DU ¼ RHS þ EAF ; ox oy where I denotes the unit or identity matrix, d þ
¼ 1; RHS ¼ k


2 n  oU o2 Un n n þ þ S þ khT ox2 oy 2

ð7Þ ð8Þ

1592

J.I. Ramos / Appl. Math. Comput. 174 (2006) 1586–1608

and EAF denotes the factorization errors, i.e.,
4  2 o o2 n n o n 2 I þ
J þ dJ I þ d
ðJ Þ EAF ¼ k 2 h2 DU. ox2 oy 2 ox2 oy 2 Eq. (6) can also be written as
 o2 Lx ðDU Þ  I
kh 2 I
kdhJn DU ¼ RHS þ EAF ; ox
Ly ðDUÞ 

 o2 n I
kh 2 I
k
hJ DU ¼ DU ; oy

ð9Þ

ð10Þ

ð11Þ

which represent linear systems of one-dimensional operators in the x- and ydirections, respectively. However, since EAF depends on DU (cf. Eq. (9)), Eqs. (10) and (11) are coupled, unless the approximation factorization errors which are O(k2) are neglected in Eq. (10). When the approximation factorization errors are disregarded, one obtains approximate factorization techniques. Eq. (10) contains a mixed fourth-order derivative which was not present in Eq. (1) (cf. Eq. (9)); this mixed derivative may be eliminated as follows. From Eq. (11), one can easily obtain DG 

o2 DU 1 ¼ ½ðI
k
hJn ÞDU
DU ; 2 oy kh

ð12Þ

which substituted into Eq. (10) yields o2 Lx ðDU Þ ¼ RHS þ kh 2 ðDU
DU Þ ox
 o2 2 þ k 2 h2 dJn 2 I þ d
ðJn Þ DU; oy

ð13Þ

which only contains second-order spatial derivatives. In order to solve Eqs. (11) and (13) and thus account for the approximate factorization errors, one can use an iterative predictor–corrector strategy [9,12]. Note, however, that, due to the coupling between Eqs. (11) and (13), this iterative technique requires the solution of these two equations. Note also that the factorization method presented here represents a one-parameter family of full factorization techniques since it depends on d (cf. Eq. (7)). Note that the one-dimensional operators of Eqs. (11) and (13) include both diffusion and reaction. For the sake of convenience, we shall refer to the linearly implicit full factorization method presented in this section with h = 0.5 as LF1. If the approximate factorization errors in this method are neglected, the one-dimensional operators in the x- and y-directions are uncoupled and no iterative procedure is required for its solution.

J.I. Ramos / Appl. Math. Comput. 174 (2006) 1586–1608

1593

3.2. Second factorization method The full factorization method presented here is analogous to the one described above except that the nonlinear terms Sn+1 are not initially linearized, i.e., the time discretization of Eq. (1) is first written as DU o2 DU o2 Un o2 DU o2 Un ¼h þ þ ð1
hÞSn þ hSnþ1 þ h þ ; 2 2 k ox oy 2 ox oy 2 which can be written in full factorization form as

 o2 o2 I
kh 2 I I
kh 2 I DU ¼ RHS þ EAF ; ox oy

ð14Þ

ð15Þ

or
 o2 Lx ðDU Þ  I
kh 2 I DU ¼ RHS þ EAF ; ox 




ð16Þ

 o2 I DU ¼ DU ; oy 2

ð17Þ


2 n  oU o2 Un n nþ1 þ þ ð1
hÞS þ hS ox2 oy 2

ð18Þ

Ly ðDUÞ 

I
kh

where RHS ¼ k

and EAF denotes the approximate factorization errors, i.e., EAF ¼ k 2 h2

o4 DU. ox2 oy 2

ð19Þ

Although the approximate factorization errors of Eq. (16) are much simpler than those of the technique presented in the previous section, these terms involve a mixed fourth-order spatial derivative, and Eqs. (16) and (17) must be solved iteratively due to the presence of Sn+1 in Eq. (18) and the factorization errors in Eq. (16). The mixed fourth-order derivative terms can be eliminated as follows. From Eq. (17), one can easily obtain DG 

o2 DU 1 ¼ ½DU
DU ; oy 2 kh

ð20Þ

which substituted into Eq. (16) yields Lx ðDU Þ ¼ RHS þ kh

o2 ðDU
DU Þ; ox2

which only contains second-order spatial derivatives.

ð21Þ

1594

J.I. Ramos / Appl. Math. Comput. 174 (2006) 1586–1608

Linearization of Sn+1, i.e., Sn+1 = Sn + k Tn + JnDU + O(k2), transforms Eq. (21) into Lx ðDU Þ ¼ RHS þ khJn DU þ kh

o2 ðDU
DU Þ; ox2

ð22Þ

where


2 n  oU o2 Un n n RHS ¼ k þ þ S þ khT ð23Þ ox2 oy 2 and the solution of the coupled Eqs. (17) and (22) provides the value of DU. The approximation method presented in this section requires iterations because Eqs. (17) and (22) are coupled even when the approximate factorization errors, i.e., the third term in the right-hand side of Eq. (22), are neglected due to the presence of the Jacobian matrix in this equation. Note that the onedimensional operators of Eqs. (17) and (22) include only diffusion. For the sake of convenience, we shall refer to the linearly implicit full factorization method presented in this section with h = 0.5, as LF2. 4. Linearly implicit approximate factorization methods In this section, we present some linearly implicit approximate factorization methods for Eq. (5). 4.1. First linearly implicit approximate factorization method This method can be easily obtained from the first full factorization technique presented above by neglecting the approximate factorization errors and can be written as
 o2  n Lx ðDU Þ  I
kh 2 I
kdhJ DU ¼ RHS; ð24Þ ox
Ly ðDUÞ 

 o2 n I
kh 2 I
k
hJ DU ¼ DU ; oy

ð25Þ

where RHS ¼ k


2 n  oU o2 Un n n þ þ S þ khT . ox2 oy 2

ð26Þ

It must be pointed out that Eqs. (24) and (25) are uncoupled and can be solved sequentially. These equations also represent a one-parameter family of methods depending on d, i.e., the allocation of the reaction terms to the x- and y-operators. Moreover, Eqs. (24) and (25) include the transient, diffusion and reaction terms in the Lx and Ly one-dimensional operators.

J.I. Ramos / Appl. Math. Comput. 174 (2006) 1586–1608

1595

Hereon, we shall refer to the approximate factorization methods presented in this section with d = 0.5 as LA1. 4.2. Second linearly implicit approximate factorization method This method can be easily obtained from the second full factorization technique presented above by neglecting the approximate factorization errors and can be written as
 o2 Lx ðDU Þ  I
kh 2 I DU ¼ RHS þ hJn DU; ð27Þ ox


 o2 I
kh 2 I DU ¼ DU ; oy

ð28Þ


2 n  oU o2 Un n n þ þ S þ khT . ox2 oy 2

ð29Þ

Ly ðDUÞ  where RHS ¼ k

Note that Eqs. (27) and (28) are linearly coupled and, therefore, an iterative procedure must be used to solve them. In this paper, we have used the predictor–corrector strategy developed by the author [9,12]. 4.3. Third linearly implicit approximate factorization method Eq. (2) can be written as
 o2 o2 I
kh 2 I
kh 2 I Unþ1 ¼ RHS; ox oy

ð30Þ

where
2 n  oU o2 U n RHS ¼ U þ kð1
hÞ þ þ kð1
hÞSn þ khSnþ1 . ox2 oy 2 n

ð31Þ

Eq. (30) can be written in factorization form as Lx ½Ly ðUnþ1 Þ ¼ RHS þ E1 ;

ð32Þ

Lx ½Ly ðUnþ1 Þ ¼ M x ½M y ðUn Þ þ E1 þ E2 þ kð1
hÞSn þ khSnþ1 ;

ð33Þ

or

where
 o2 Lx V  I
kh 2 I V; ox

ð34Þ

1596

J.I. Ramos / Appl. Math. Comput. 174 (2006) 1586–1608


 o2 Ly V  I
kh 2 I V; oy
M xV 

Ly V 

I þ kð1
hÞ

 o2 I V; ox2


 o2 I þ kð1
hÞ 2 I V; oy

ð35Þ

ð36Þ

ð37Þ


 o2 o2 Unþ1 E1 ¼ k h 2 ; ox oy 2

ð38Þ


 o2 o2 U n E2 ¼
k ð1
hÞ 2 . ox oy 2

ð39Þ

2 2

2

2

If the O(k2) terms in Eq. (33) are neglected, one can easily show that this equation can be written in approximate factorization form (with accuracy of O(k2)) as Lx ðU Þ ¼ M x ½M y ðUn Þ þ kð1
hÞSn ;

ð40Þ

Ly ðUnþ1 Þ ¼ U þ khSnþ1

ð41Þ

and the last equation must be solved iteratively due to the nonlinearity of the source term S. The need for iterations can be avoided completely by linearizing the source term with respect to the previous time level and writing Eq. (41) as ðLy
khJn ÞUnþ1 ¼ U þ khðSn
Jn Un Þ;

ð42Þ

The linearly implicit approximate factorization method corresponding to Eqs. (40) and (42) with h = 0.5 is hereafter referred to as LA2. Note that the x-operator involves diffusion, whereas the y-one involves both diffusion and reaction.

5. Results In this section, we present some sample results corresponding to N = 2,
20 6 x 6 20,
20 6 y 6 20, U = (u, v)T, where T denotes transpose, the following source terms: S u ¼
u6 v;

S v ¼ u6 v
0:5v;

ð43Þ

initial conditions uðx; y; 0Þ ¼ 1;

vðx; y; 0Þ ¼ expð
ðx2 þ y 2 ÞÞ

ð44Þ

J.I. Ramos / Appl. Math. Comput. 174 (2006) 1586–1608

1597

and Dirichlet boundary conditions uoD = 1 and voD = 0, where oD denotes the boundary of the domain D = [
20, 20] · [
20, 20]. Figs. 1 and 2 show the profiles of u and v, respectively, at selected times, as functions of x and y and correspond to LA1, equally-spaced grids, second-order accurate discretizations of the diffusion terms, h = 0.5, k = 0.02 and NI = NJ = 102, where NI and NJ denote the number of grid points in the x- and y-directions, respectively. Note that in these and the following figures only the results corresponding to even number of grid points are drawn. Fig. 1 shows that u is being consumed in the center of the computational domain as a consequence of the increase of v there. Both the depth and the width of the valley of the u profile increase as time increases. Fig. 2 shows that the amplitude of the initial peak of the v profile decreases as time increases. This peak flattens and becomes a valley; both the peak and valley propagate from the center to the boundaries of the domain. In order to assess the numerical errors of the LA1, LA2 and LF2, we have compared the results of these methods with those obtained with LF1. Remember that LF1 and LF2 are linearly implicit full factorization techniques which

Fig. 1. u as a function of x and y obtained with LA1, k = 0.02, NI = NJ = 102 (from left to right, top to bottom: t = 4, 6, 8, 10, 12, 14, 16, 18, 20).

1598

J.I. Ramos / Appl. Math. Comput. 174 (2006) 1586–1608

Fig. 2. v as a function of x and y obtained with LA1, k = 0.02, NI = NJ = 102 (from left to right, top to bottom: t = 4, 6, 8, 10, 12, 14, 16, 18, 20).

were solved iteratively by means of the predictor–corrector method developed by the author [9,12] using the following convergence criterion: NI X NJ X ðDUkþ1 ÞT  DUkþ1
ðDUk ÞT  DUk 6 10
16 ; NI  NJ i¼1 j¼1

ð45Þ

where the superscript k denotes the kth iteration within the time step, whereas LA1 and LA2 are linearly implicit approximate factorization techniques that uncouple the solution of the x-operator from that of the y-operator. The error of the full and approximate factorization methods presented in this paper have been evaluated as duM  EUM = uM(x, y, t)
uLF1(x, y, t) and dvM  EVM = vLF1(x, y, t)
vM(x, y, t), where the subscript LF1 denotes the solution obtained with LF1 and M, those obtained with LA1, LA2 and LF2. Figs. 3–8 illustrate the errors for k = 0.02 and NI = NJ = 102 and indicate that the largest errors of LA1 are, at most, 5 · 10
3 at t = 6 and, at most, 1 · 10
3 for t > 6, for the u component (cf. Fig. 3). On the other hand, the errors of LA2 are, at most, 5 · 10
6 at t = 6, 2 · 10
6 at t = 8, 10, 16, 18 and 20, and, at most, 5 · 10
7 at t = 12 and 14, for the u component (cf. Fig. 4); therefore, the errors of LA2 are not monotonic functions of time. The errors of LF2 are, at most, 1 · 10
3 for t P 6 and, therefore, are almost not functions of time. Fig. 6 indicates that the errors in v predicted by LA1 are, at most, 2 · 10
3, at t = 6, and, at most, 1 · 10
3 for t > 6. On the other hand, the errors in v of

J.I. Ramos / Appl. Math. Comput. 174 (2006) 1586–1608

1599

Fig. 3. Errors in u, i.e., EU1 = du = uLF1
uLA1, where LF1 and LA1 denote the solutions obtained with LF1 and LA1, respectively, k = 0.02, NI = NJ = 102 (from left to right, top to bottom: t = 4, 6, 8, 10, 12, 14, 16, 18, 20).

LA2 are, at most, 5 · 10
6, at t = 6, at most, 5 · 10
7 for t = 8, 16 and 18, and, at most, 2 · 10
7 for t = 12 and 14 as shown in Fig. 7, whereas the errors of LF2 are, at most, 1 · 10
4, at t = 6, 8 and 10, at most, 1 · 10
4, at t = 6, 8 and 10, at most, 2 · 10
4 for t = 12, 14 and 16, and, at most, 1 · 10
3, at t = 18 and 20 as illustrated in Fig. 8. Figs. 4 and 7 indicates that the errors of LA2 are not symmetric with respect to the x- and y-axes, whereas those of LA1 and LF2 are symmetric as illustrated in Figs. 3 and 6 and Figs. 5 and 8, respectively. Since for the problem considered in this paper, the solution should be symmetric with respect to the x- and y-axes, it can be concluded that LA2 does not provide symmetric solutions, and this drawback is associated with Eqs. (40) and (42) that show that the x-operator includes diffusion, whereas the y-one includes both diffusion and reaction. However, as the results presented in Figs. 4 and 7 indicate this method provides more accurate answers with respect to LF1 than LA1 and LF2. The results presented in Figs. 3–8 indicate that the accuracy of LF2 is O(k2) compared to that LF1 as one should have expected, for these two methods are second-order accurate (in time) full factorization techniques. Figs. 3–8 also indicate that, for v, LA2 is more accurate than LF2 by, at least, two orders of magnitude, and LF2 is more accurate than LA1, despite the fact that LA1 and LA2 are linearly implicit approximate factorization techniques of O(k2)

1600

J.I. Ramos / Appl. Math. Comput. 174 (2006) 1586–1608

Fig. 4. Errors in u, i.e., EU2 = du = uLF1
uLA2, where LF1 and LA2 denote the solutions obtained with LF1 and LA2, respectively, k = 0.02, NI = NJ = 102 (from left to right, top to bottom: t = 4, 6, 8, 10, 12, 14, 16, 18, 20).

accuracy in time. For u, the results presented in Figs. 3–5 indicate that LF2 is approximately five times more accurate than LA1 but, at least two orders of magnitude less accurate than LA2. Therefore, when solving systems of two-dimensional reaction–diffusion equations by means of (iterative) linearly implicit full factorization and approximate factorization techniques, one should assess both the accuracy and the symmetry of the results for the different components of the vector U as functions of space and time, and, as the results presented above indicate, the accuracy of approximate factorization depends on the allocation of the source terms to either the x- or y-operators, even when these techniques have the same order of temporal accuracy. This point is again emphasized in Tables 1 and 2 which show the errors of the factorization methods presented in this paper as functions of time and k. It must be noted that, unless stated otherwise, the errors of the methods considered in this paper were nil for t 6 4. Tables 1 and 2 show that the errors of LF2, LA1 and LA2 with respect to LF1 are not monotonic functions of the time step and time, and that the relative errors of LF2 are smaller than those of LA1, whereas those of LA2 are the smallest for the three time steps considered in this paper. Although the numerical errors of LA1, LA2, LF1 and LF2 decrease in a quadratic fashion as the time step is decreased, the results presented in Tables 1

J.I. Ramos / Appl. Math. Comput. 174 (2006) 1586–1608

1601

Fig. 5. Errors in u, i.e., EU3 = du = uLF1
uLF2, where LF1 and LF2 denote the solutions obtained with LF1 and LF2, respectively, k = 0.02, NI = NJ = 102 (from left to right, top to bottom: t = 4, 6, 8, 10, 12, 14, 16, 18, 20).

and 2 indicate that the difference between the solutions obtained with these methods is, in most cases, a non-monotonic function of the time step and time. Moreover, these tables also show that, if the solution obtained with LF1 is assumed to be exact, the results obtained with LA2 are more accurate than those of LA1 and LF2 for k = 0.04, 0.02 and 0.01, despite the fact that this method does not account for the factorization errors. Since both LF1 and LF2 are linearly implicit fully factorized methods which are second-order accurate in time and, therefore, their solutions can be written as ULF1 = Ue + CLF1k2 and ULF2 = Ue + CLF2k2, respectively, then ULF1
ULF2 = C4k2, where the subscript e denotes the exact solution, C4 = CLF1
CLF2, and the results presented in Figs. 3–8 and Tables 1 and 2 show that C4 is a function of both t and k. An analogous comment applies to LA1 and LA2. This implies that the measurement of the accuracy based on the order of magnitude of the time step, i.e., O(k2), may be misleading because it ignores the fact that the local error is, in fact, Ck2, and C may depend on time, time step, spatial step size and higher-order spatial and temporal derivatives, and, most importantly, may be very large. Methods for which C is a constant independent of time and time step are known as uniformly convergent.

1602

J.I. Ramos / Appl. Math. Comput. 174 (2006) 1586–1608

Fig. 6. Errors in v, i.e., EV1 = dv = vLF1
vLA1, where LF1 and LA1 denote the solutions obtained with LF1 and LA1, respectively, k = 0.02, NI = NJ = 102 (from left to right, top to bottom: t = 4, 6, 8, 10, 12, 14, 16, 18, 20).

The validity of the comment made in the previous paragraph has also been verified with the following iterative full and approximate factorization methods which do not use time linearization. 5.1. First iterative full factorization method In this technique, the nonlinear Eq. (2) was first written as

 o2 o2 Lx ½Ly ðUkþ1 Þ  I
kh 2 I
kdhJk I
kh 2 I
k
hJk Ukþ1 ox oy ¼ RHS þ EAF ;

ð46Þ

d þ
¼ 1;

ð47Þ


RHS ¼ kð1
hÞ

o2 Un o2 Un þ þ Sn ox2 oy 2



þ khðSk
Jk Uk Þ

and EAF denotes the approximate factorization errors, i.e.,
4  2 o o2 k 2 2 k o k 2 I þ
2 J þ dJ I þ d
ðJ Þ Ukþ1 . EAF ¼ k h ox2 oy 2 ox oy 2

ð48Þ

ð49Þ

J.I. Ramos / Appl. Math. Comput. 174 (2006) 1586–1608

1603

Fig. 7. Errors in v, i.e., EV2 = dv = vLF1
vLA2, where LF1 and LA2 denote the solutions obtained with LF1 and LA2, respectively, k = 0.02, NI = NJ = 102 (from left to right, top to bottom: t = 4, 6, 8, 10, 12, 14, 16, 18, 20).

Eq. (46) can be written as the following sequence of one-dimensional operators
 o2 ð50Þ Lx ðU Þ  I
kh 2 I
kdhJk U ¼ RHS þ EAF ; ox Ly ðU

kþ1


 o2 k Þ  I
kh 2 I
k
hJ Ukþ1 ¼ U . oy

ð51Þ

The mixed fourth-order spatial derivative in Eq. (50) can be eliminated in an analogous manner to that shown previously, so that this equation can be expressed as  o2  Lx ðU Þ ¼ RHS þ kh 2 Ukþ1
U ox
 2 2 2 k o k 2 þ k h dJ I þ d
ðJ Þ Ukþ1 oy 2

ð52Þ

and Eqs. (51) and (52) can be solved iteratively by means of the predictor–corrector strategy developed by the author [9,12] until an analogous convergence criterion to the one presented before is satisfied.

1604

J.I. Ramos / Appl. Math. Comput. 174 (2006) 1586–1608

Fig. 8. Errors in v, i.e., EV3 = dv = vLF1
vLF2, where LF1 and LF2 denote the solutions obtained with LF1 and LF2, respectively, k = 0.02, NI = NJ = 102 (from left to right, top to bottom: t = 4, 6, 8, 10, 12, 14, 16, 18, 20).

Table 1 Errors in u of LA1, LA2 and LF2 with respect to LF1, i.e., EUi = uLF1
uM, where M = 1, 2 and 3 correspond to LA1, LA2 and LF2, respectively Error

t=6

EU1 (k = 0.04) EU1 (k = 0.02) EU1 (k = 0.01)

5 · 10
3 5 · 10
3 2 · 10
3 2 · 10
3 2 · 10
3 2 · 10
3 2 · 10
3 2 · 10
3 5 · 10
3 1 · 10
3 1 · 10
3 1 · 10
3 1 · 10
3 1 · 10
3 1 · 10
3 1 · 10
3 2 · 10
3 1 · 10
3 1 · 10
3 1 · 10
3 1 · 10
3 1 · 10
3 1 · 10
3 1 · 10
3

t=8

t = 10

t = 12

t = 14

t = 16

t = 18

t = 20

EU2 (k = 0.04) EU2 (k = 0.02) EU2 (k = 0.01)

1 · 10
5 5 · 10
6 5 · 10
6 2 · 10
6 2 · 10
6 5 · 10
6 5 · 10
6 5 · 10
6 5 · 10
6 2 · 10
6 2 · 10
6 5 · 10
7 5 · 10
7 2 · 10
6 2 · 10
6 2 · 10
6 2 · 10
6 5 · 10
7 2 · 10
7 1 · 10
7 1 · 10
7 2 · 10
7 2 · 10
7 5 · 10
7

EU3 (k = 0.04) EU3 (k = 0.02) EU3 (k = 0.01)

2 · 10
3 1 · 10
3 1 · 10
3 1 · 10
3 1 · 10
3 1 · 10
3 1 · 10
3 2 · 10
3 1 · 10
3 1 · 10
3 1 · 10
3 1 · 10
3 1 · 10
3 1 · 10
3 1 · 10
3 1 · 10
3 1 · 10
3 2 · 10
4 2 · 10
4 2 · 10
4 2 · 10
4 1 · 10
3 1 · 10
3 1 · 10
3

The main inconvenience of the first iterative full factorization method just described is that the operators Lx and Ly contain the Jacobian matrix that changes from one iteration to the next one. This is avoided in the second iterative full factorization method presented next.

J.I. Ramos / Appl. Math. Comput. 174 (2006) 1586–1608

1605

Table 2 Errors in v of LA1, LA2 and LF2 with respect to LF1, i.e., EVi = vLF1
vM, where M = 1, 2 and 3 correspond to LA1, LA2 and LF2, respectively Error

t=6

t=8
3

t = 10
3

t = 12


3


3

t = 14
3

t = 16
3

t = 18

t = 20


3

EV1 (k = 0.04) EV1 (k = 0.02) EV1 (k = 0.01)

5 · 10 2 · 10 1 · 10 1 · 10 1 · 10 1 · 10 1 · 10 1 · 10
3 2 · 10
3 1 · 10
3 1 · 10
3 1 · 10
3 1 · 10
3 1 · 10
3 1 · 10
3 1 · 10
3 1 · 10
3 1 · 10
3 2 · 10
4 2 · 10
4 2 · 10
4 1 · 10
3 1 · 10
3 1 · 10
3

EV2 (k = 0.04) EV2 (k = 0.02) EV2 (k = 0.01)

1 · 10
5 2 · 10
6 2 · 10
6 2 · 10
6 2 · 10
6 2 · 10
6 2 · 10
6 2 · 10
6 5 · 10
6 5 · 10
7 2 · 10
7 2 · 10
7 2 · 10
7 5 · 10
7 5 · 10
7 5 · 10
7 2 · 10
6 2 · 10
7 5 · 10
8 5 · 10
8 5 · 10
8 1 · 10
7 1 · 10
7 1 · 10
7

EV3 (k = 0.04) EV3 (k = 0.02) EV3 (k = 0.01)

2 · 10
4 2 · 10
4 2 · 10
4 1 · 10
3 1 · 10
3 1 · 10
3 1 · 10
3 2 · 10
3 1 · 10
4 1 · 10
4 1 · 10
4 2 · 10
4 2 · 10
4 2 · 10
4 1 · 10
3 1 · 10
3 5 · 10
5 5 · 10
5 5 · 10
5 1 · 10
4 1 · 10
4 1 · 10
4 2 · 10
4 2 · 10
4

5.2. Second iterative full factorization method In this technique, the nonlinear Eq. (2) is first written as

 o2 o2 Lx ½Ly ðUkþ1 Þ  I
kh 2 I I
kh 2 I Ukþ1 ¼ RHS þ EAF ; ox oy d þ
¼ 1;

ð53Þ ð54Þ




o2 Un o2 Un RHS ¼ kð1
hÞ þ þ Sn ox2 oy 2



   þ kh Sk þ Jk Ukþ1
Uk

ð55Þ

and EAF denotes the approximate factorization errors, i.e., EAF ¼ k 2 h2

o4 Ukþ1 . ox2 oy 2

ð56Þ

Eq. (53) can be written as the following sequence of one-dimensional operators
 o2 Lx ðU Þ  I
kh 2 I U ¼ RHS þ EAF ; ð57Þ ox Ly ðU

kþ1


 o2 Þ  I
kh 2 I Ukþ1 ¼ U . oy

ð58Þ

The mixed fourth-order spatial derivative in Eq. (57) can be eliminated in an analogous manner to that shown previously so that this equation can be expressed as Lx ðU Þ ¼ RHS þ kh

 o2  kþ1 U
U 2 ox

ð59Þ

1606

J.I. Ramos / Appl. Math. Comput. 174 (2006) 1586–1608

and Eqs. (58) and (59) can be solved iteratively by means of the predictor–corrector strategy developed by the author [9,12] until an analogous convergence criterion to the one presented before is satisfied. The operators in Eqs. (58) and (59) do not depend on the Jacobian matrix, i.e., they are constant, and can be inverted and stored in the iterative procedure. The first and second iterative approximate factorization methods considered in this paper can be easily obtained from the first and second iterative full approximation techniques described above, respectively, by simply neglecting the factorization errors, i.e., the second term in the right-hand side of Eqs. (52) and (59), respectively. Calculations performed with the first and second iterative full approximate factorization methods show second-order temporal accuracy for h = 0.5 and similar trends to those of LF1 and LF2, thus indicating that, for the mildly nonlinear system of reaction–diffusion equations considered in this paper, there are few differences between the results corresponding to iterative and linearly-implicit full factorization methods. Analogous results were obtained with the first and second iterative approximate factorization methods. Therefore, for the problem considered here, one does not need to iterate the reaction terms, and linearly-implicit full and approximate factorization methods yield nearly the same results as iterative full and approximate factorization methods, respectively. However, the CPU times of iterative factorization methods are larger than those of linearly-implicit factorization techniques as illustrated in Table 3, where the first and second iterative full factorization methods are referred to as IF1 and IF2, respectively, while the corresponding iterative approximate factorization algorithms are referred to as IA1 and IA2, respectively, when h = 0.5. The CPU times presented in Table 3 refer to equally-spaced grids with NI = NJ = 102 and different time steps, and show that linearly-implicit full factorization methods require about twice as much CPU time than the corresponding approximate factorization ones. The CPU times of IA1, IA2, IF1 and IF2 are smaller than twice those of LA1, LA2, LF1 and LF2, respectively, thus indicating that iterative factorization methods converge on average in two iterations with the convergence criterion employed in this paper.

Table 3 CPU times (seconds) on a HP J5000 at 440 MHz Method

LA1

LA2

LF1

LF2

IF1

IF2

IA1

IA2

k = 0.04 k = 0.02 k = 0.01

74.4 148.1 295.5

68.4 135.9 272.1

165.9 314.6 615.6

131 255.2 503.4

119.4 254.9 318.3

107.6 229.5 309.2

306.7 576.3 637.5

214.1 493.1 526.9

J.I. Ramos / Appl. Math. Comput. 174 (2006) 1586–1608

1607

6. Conclusions Full and approximate factorization methods for the numerical solution of multidimensional reaction–diffusion equations have been developed. Full factorization methods account for the factorization errors and may be linear and nonlinear. In linear full factorization methods, the original reaction–diffusion equation is linearized with respect to the previous time level and these methods provide linear systems of algebraic equations at each time step. However, owing to the inclusion of the factorization errors, an iterative procedure is required to determine the solution. Nonlinear full factorization methods require iterations due to both the factorization errors and the nonlinearity of the original operator. Approximate factorization techniques do not account for the factorization errors, have the same order of temporal and spatial accuracy as full factorization methods, and may be linear and nonlinear. Linearly implicit approximate factorization methods yield uncoupled linear equations methods in each direction, whereas nonlinear ones require iterations due to the nonlinearity of the source terms. In this paper, it has been shown that the factorization of a multidmensional reaction–diffusion operator is not unique, and several full and approximate, linear and nonlinear, factorization methods have been presented. It has also been shown that, both linearly implicit and nonlinear, approximate factorization errors may not preserve the symmetry of the original differential problem due to the allocation of the reaction and diffusion terms to the one-dimensional operators, even though the errors of these techniques may be small. It has also been shown that iterative full and approximate factorization methods do converge in about two iterations for the two-dimensional system of reaction–diffusion equations considered in this paper, and provide almost the same solution as linearly implicit full and approximate factorization techniques. It has also been shown that the formal order of accuracy may be misleading because the ‘‘constant’’ that multiplies this order may not be constant, and may be very large. Acknowledgements 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] D.W. Peaceman, H.H. Rachford Jr., The numerical solution of parabolic and elliptic partial differential equations, J. Soc. Ind. Appl. Math. 3 (1955) 28–41.

1608

J.I. Ramos / Appl. Math. Comput. 174 (2006) 1586–1608

[2] J. Douglas Jr., On the numerical integration of uxx + uyy = ut, J. Soc. Ind. Appl. Math. 3 (1955) 42–65. [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] W.R. Briley, H. McDonald, On the structure and use of linearized block implicit methods, J. Comput. Phys. 34 (1980) 54–73. [6] R.F. Warming, R.M. Beam, An extension of A-stability to alternating direction implicit methods, BIT 19 (1979) 395–417. [7] P.J. van der Houwen, B.P. Sommeijer, Approximate factorization for time-dependent partial differential equations, J. Comput. Appl. Math. 128 (2001) 447–466. [8] M.A. Botchev, J.G. Verwer, A new approximate matrix factorization for implicit time integration in air pollution modeling, J. Comput. Appl. Math. 157 (2003) 309–327. [9] J.I. Ramos, Linearized factorization techniques for multidimensional reaction–diffusion equations, Appl. Math. Comput. 100 (1999) 201–222. [10] J.I. Ramos, Linearization methods for reaction–diffusion equations: 1-D problems, Appl. Math. Comput. 88 (1997) 199–224. [11] J.I. Ramos, Linearization methods for reaction–diffusion equations: multidimensional problems, Appl. Math. Comput. 88 (1997) 225–254. [12] J.I. Ramos, Implicit, compact, linearized h-methods with factorization for multidimensional reaction–diffusion equations, Appl. Math. Comput. 94 (1998) 17–43. [13] 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. [14] W. Liao, J. Zhu, A.Q.M. Khaliq, An efficient high-order algorithm for solving systems of reaction–diffusion equations, Numer. Methods Partial Differen. Equat. 18 (2002) 340–354. [15] J.I. Ramos, On diffusive methods and exponentially-fitted techniques, Appl. Math. Comput. 103 (1999) 69–96. [16] J.I. Ramos, Exponential methods for one-dimensional reaction–diffusion equations, Appl. Math. Comput., in press, doi:10.1016/j.amc.2004.12.003. [17] R.W. MacCormack, Iterative modified approximate factorization, Comput. Fluids 30 (2001) 917–925. [18] E. Giladi, J.B. Keller, Iterative solution of elliptic problems by approximate factorization, J. Comput. Appl. Math. 85 (1997) 287–313. [19] C. Eichler-Liebenow, P.J. van der Houwen, B.P. Sommeijer, Analysis of approximate factorization in iteration methods, Appl. Numer. Anal. 28 (1998) 245–258. [20] W.R. Briley, H. McDonald, An overview and generalization of implicit Navier–Stokes algorithms and approximate factorization, Comput. Fluids 30 (2001) 807–828. [21] E. Steinthorsson, T.I.-P. Shih, Methods for reducing approximate-factorization errors in twoand three-factored schemes, SIAM J. Sci. Comput. 14 (1993) 1214–1236.