Lagged diffusivity fixed point iteration for solving steady-state reaction diffusion problems

The paper concerns with the computational algorithms for a steady-state reaction diffusion problem. A lagged diffusivity iterative algorithm is proposed for solving resulting system of quasilinear equations from a finite difference discretization. The convergence of the algorithm is discussed and the numerical results show the efficiency of this algorithm.


Introduction
We consider the problem of finding a solution of the system of nonlinear equations where F : ⊆ R n → R n is a continuously differentiable mapping.
We are interested in large-scale systems for which the Jacobian of F is not available or is difficult to compute.
Particularly, we consider the solution of the system of quasilinear equations

Statement of the problem
We are concerned with the numerical solution of the steady-state reaction diffusion problem by the finite difference method.
The functions σ (x, y, u), g(x, y, u) and s(x, y) are assumed to satisfy the following Smoothness Conditions for x and y in R ∪ ∂R and u in a neighbourhood of a solution of problem (2)- (3). We assume that this model problem has a solution.
(i) The functions σ (x, y, u) and g(x, y, u) are piecewise continuous in x and y and continuous in u; the source term s(x, y) is piecewise continuous in x, y. (ii) There exist two positive constants σ min and σ max such that 0 < σ min ≤ σ (x, y, u) ≤ σ max , uniformly in x, y and u.
(iii) For fixed x and y, the function σ (x, y, u) is locally Lipschitz continuous at u (uniformly in x and y), with constant > 0. (iv) For fixed x and y, the function g(x, y, u) is uniformly monotone 1 at u (uniformly in x and y), with constant c > 0 and is continuously differentiable at u.
There exist various techniques for discretizing the problem (2)- (3). The use of the truncated Taylor series to represent the derivative in Equation (2) provides more insight into the nature of the truncation error that arises when the continuous model (2) is replaced by a discrete set of finite difference equations.
Using the Taylor series approach, Equation (2) will be solved with the following standard finite difference scheme.
If the domain R is the unit square, we superimpose on R a uniform grid of points R h . More precisely, R h is defined in this manner. In the plane Oxy, we draw a system of straight lines parallel to the coordinate axes: x = x i and y = y j with mesh spacings uniform in each coordinate direction. We denote the intersection of these lines as mesh points or grid points. The coordinates of the interior mesh points (x i , y j ) of R h are given by x i = i x and y j = j y with i = 1, . . . , N x and j = 1, . . . , N y and x = 1/(N x + 1), y = 1/(N y + 1).
For simplicity, we assume x = y = h and consequently N x = N y = N. At the mesh points ofR h = R h ∪ ∂R h , the function u(x i , y j ) is to be approximated by a grid function u ij . In order to approximate partial derivatives, we shall make use of various difference quotients of grid functions. The forward, backward and centred difference quotients with respect to x and y of the grid function {u ij } at the mesh point (x i , y j ) are, respectively, , δ y u ij = 1 2 ( y u ij + ∇ y u ij ). This convenient notation was introduced by Courant et al. [3].
A particular important case is the centred second difference quotients which can be written as A centred finite difference approximation of Equation (2) at the mesh point (x i , y j ) may be written as The Dirichlet boundary condition (3) produces the following conditions: The total number of difference equations is n = N × N.
These difference equations may be written in a well-known matrix form. The above five point discretization formula which approximates the partial derivatives transforms (2) and (3) into a system of nonlinear n algebraic equations as (1) where the vector u is the approximation of the restriction of u(x, y) on R h . For natural ordering by horizontal lines of the mesh points, the irreducible n × n matrix A(u) can be partitioned in the form of a block tridiagonal matrix of order n, where each square block on the diagonal of A(u) is a diagonally dominant tridiagonal matrix of order N and each block on the subdiagonal and on the superdiagonal of A(u) is a diagonal matrix of order N. That is, we can express the matrix A(u) as Using only one index for numbering the mesh points P k , k = 1, . . . , n, where k = (j − 1)N + i, the elements a kl (u) of the kth row, k = 1, . . . , n, of the matrix A(u) are given by Now we can define the vector s ∈ R n whose component s k , k = 1, . . . , n is the value of the source term s(x, y) at the mesh point P k of R h . Analogously we define the nonlinear mapping G(u) whose kth component G k (u) is a function only of the coordinates of the mesh point P k of R h and of the approximation of u(x, y) at P k , k = 1, . . . , n.
Thus, we can write Equations (4) and (5) in the matrix form (1) This system of nonlinear equations is characterized by the following Properties (see, e.g. [12]).
(I) The system (1) has at least one solution and all the solutions belong to a well-defined closed ball in R n where ρ is independent of h and of s h .
Here, u h is the discrete l 2 (R h ) norm of grid functions {u ij }: (II) At each mesh point (x i , y j ) ofR h the backward difference quotients with respect to x and to y of a grid function {u ij } belonging to are bounded. The bound is independent of h, but depends on s h .
withc and c positive constants independent of h for which Proposition 2 Let u be a solution of the system of 1-D finite difference Equation (1). Then

Monotonicity of the mapping F
We begin with two useful lemmas which establish the discrete analogs of two continuous L 2 (R) inner products used frequently in the study of elliptic partial differential equations. and then, we have the result (7). Similarly, we obtain formula (8).

Then, we have the following expression for the discrete l 2 (R h ) inner product of the vectors A(w)u and v :
where ∇ h u ij denotes the gradient of the grid function {u ij } at the mesh point (x i , y j ) (see formulae (10) and (11)).
Proof Since the discrete l 2 (R h ) inner product of two grid functions {u ij } and {v ij } which are zero on the boundary ∂R is defined by we have We have used formulae (7) and (8) and the hypothesis that the grid functions {u ij } and {v ij } are zero at the boundary ∂R, i.e. satisfy the condition (5) on ∂R.
We define the gradient of the grid function {u ij } at (i, j) entry to be Thus, we can write Therefore, we obtain the result (9).
Remark 1 While the grid function {u ij } is defined on the whole mesh regionR h , the vector u ∈ R n represents the grid function {u ij } defined only on the interior mesh points R h . Now we want to prove that under suitable conditions the mapping (1) is uniformly monotone on the closed ball (see, Property (I)).
Let {u ij } and {v ij } be two grid functions belonging to , i.e. two solutions of the nonlinear system (1). For the condition (5), they are represented by the vectors u and v of n components.

As consequence of the identity
we have, by Equation (1), Thus, from Equation (13), the discrete l 2 (R h ) inner product Using Equation (9), we have Assumption (iv) implies that, for all grid functions u and v belonging to , there exists a positive constant c such that for all i, j = 1, . . . , N (see [16, p. 141]). The constant c is independent of h. Thus, for the discrete l 2 (R h ) inner product (see definition (10)) and from Equation (16), we have the inequality (see, definition (6)) The next step is to find a bound for Specifically, assumption (iii) says that for given grid functions u and v belonging to , there exist some w on the line between u and v and a positive constant (w ij ) such that here, is independent of h and {w ij } ∈ . Property (II) assures that there exists a constant β > 0 such that for all i, j = 1, . . . , N + 1 and all grid function {v ij } belonging to . β is independent of h. Now, we may apply the inequalities (19) and (21) and the definition (20) to estimate expression (18), then The last two terms are zero because the grid functions {u ij } and {v ij } belonging to are bounded and satisfy the boundary condition (5).
Using a well-known technical trick involving the arithmetic mean and geometric mean inequality of positive numbers where α is as yet an undetermined positive number. It now follows from Equation (14) that and from Equations (15), (17) and (22) and assumption (ii), formula (23) becomes If we set then the inequality (24) can be written as When inequality (25) shows that the mapping F(u) is uniformly monotone on . When inequality (25) holds, there exists only one solution u * of the system F(u) = 0 in . Indeed, suppose thatũ * is some other solution of F(u) = 0 in ,ũ * = u * . Then, since the mapping F(u) is uniformly monotone on and Equation (25) which is a contradiction. Hence u * is the only solution of F(u) = 0 in . Inequality (26) says that the best chances for the applicability of the algorithm to be defined in next section, are assured when the local Lipschitz constant of σ (x, y, u) is small and the constant c (in condition (iv)) on the monotonicity of g(x, y, u) is large.

Lagged diffusivity fixed point iteration
We will investigate the solvability of the system of nonlinear difference equations (1) under the assumption that the computation of the Jacobian matrix of F(u) is not convenient. For solving this system the easiest and may be the most common method is to lag part of the nonlinear term in Equation (1) (see [20]).
We will show that under weak restrictions (i)-(iv) imposed on the functions σ (x, y, u) and g(x, y, u) the problem (2)-(3) can be solved via a sequence of systems of weakly nonlinear difference equations where only g but not σ depends on the approximate solution u of u(x, y).
That is, we will apply the so-called lagged diffusivity fixed point iteration procedure.
Many numerical experiments have shown the following facts.
1. In the special case in which the diffusion term in Equation (2) dominates the reaction term, the lagged diffusivity procedure is an efficient and robust method, even if only linearly convergent, for solving problems with a highly nonlinear differential operator −div (σ ∇u). (For example, this may happen for nonlinear diffusion models in image denoising [22].) When this operator is discretized, the matrix A(u) has highly varying coefficients and Newton's method for solving the nonlinear system (1) does not work satisfactorily, in the sense that its domain of convergence is very small. 2. In the general case, it may be noted that the steeper the slope of g (in Equation (2)) and the flatter that of σ relative to that of g, the better are the chances for the applicability of the lagged diffusivity procedure.
Specifically, let u (0) ∈ be an initial approximation to the solution of u * of the system (1). The mapping represents the discretized weakly nonlinear difference system which can be solved inexactly with a Newton iterative method. This solution u (1) produces a residual F 0 (u (1) ) with bound F 0 (u (1) ) ≤ ε 1 . Here · is the Euclidean norm.
Thus, u (1) is a new approximation to u * . Now, u (1) becomes the current iterate to generate the new iterate u (2) , and so on.
In general, if u (ν) is an estimate of the solution u * of Equation (1), we will determine a new estimate of u * by solving the weakly nonlinear system An approximate solution of the weakly nonlinear system (27) is computed by a Newton iterative method in such a way that its solution u (ν+1) will be accepted if the residual F ν (u (ν+1) ) satisfies the condition where ε ν+1 is a given tolerance such that ε ν+1 → 0 for ν → ∞.
If such suitable solution u (ν+1) is found, we say that the algorithm does not break down. The special form of Equation (1)  Now, it is important to note that the iterate u (ν+1) is the solution of a weakly nonlinear reaction diffusion equation (whose diffusivity σ depends on the previous iterate) with inhomogenous term −s − F ν (u (ν+1) ).
Since the bound of the backward difference quotients depends on the inhomogeneous term (see Proposition 2), we have that there exist two constants β andβ such that instead of Equation (21), i, j = 1, . . . , N + 1.

E. Galligani
A result concerning with the convergence of the lagged diffusivity procedure is the following theorem.
Theorem 1 Let be given the nonlinear system (1) arising from the discretization of the problem (2)-(3) subject to the Smoothness Conditions (i)-(iv), and characterized by the Properties (I) and (II), which define the closed ball and the bound β (see formula (21)).

If all the vectors {u (ν) } belong to and satisfy the Property S in with the bound (29) instead of (21), then the sequence {u (ν) } converges to the unique solution u * of F(u) = 0.
Proof The unique solution u * in of Equation (1) satisfies the equation and the iterate u (ν+1) satisfies the equation where the Euclidean norm of F ν (u (ν+1) ) satisfies the inequality (28). Subtracting Equation (31) from Equation (30) and taking into account of the identity (12), we can write
Moreover, since the grid function {u (ν+1) ij } belongs to , we may assume that Thus, we obtain the inequality and iteratively Since ε ν → 0 it follows from the general Toeplitz lemma [16, p. 399] that Hence, the sequence {u (ν) } of approximate solutions converges to the solution u * of the system (1).
Let us now turn to the numerical solution of the weakly nonlinear system (27) with a Newton iterative method.
For the solution of systems of nonlinear equations of the form (1) with a Newton iterative method a common set of Standard Assumptions on F(u) is the following.
• Equation (1) has a solution u * in an open domain K of R n .
For any vector u and u (ν) belonging to K: • A(u (ν) ) is a block tridiagonal matrix of order n for any iterate u (ν) . The diagonal blocks are square (although not necessarily all of the same order) tridiagonal submatrices, and the off-diagonal blocks are diagonal submatrices. A less restrictive assumption on G(u) may be considered (see, e.g. [6]). However, with the assumption of −G(u) being diagonal, it is possible to express the rate of change due to a reaction of many reaction diffusion processes in realistic applications. (See e.g. [1,2,10,11,13,14,17,18].) As a specific example of a solver of Equation (27), we consider the Newton-arithmetic mean (AM) method and its simplified version developed in [6][7][8][9].
The Newton-AM method belongs to the general class of Newton-Iterative methods, as, for example, the well-known Newton-SOR method (see [16,Sections 7.4,10.3]).
The Newton-AM method incorporates at each stage of Newton's method, the AM method as inner iterative solver for the linear Newton equation [8]. When we apply this method to the weakly nonlinear system (27), we can reduce it to a two-stage iterative method as that described in [6,9].
The above standard assumptions on F ν (u) permit us to analyse the convergence of the Newton-AM method and of the simplified Newton-AM method [6,8,9].
In [6][7][8][9], we have also discussed the questions concerning with the global convergence of the Newton-AM method, using the line search strategy, and the monotone convergence of the simplified Newton-AM method.
These methods are well suited for implementation on parallel computers and are significantly efficient when the nonlinear differential operator in Equation (2) is general, i.e. it is composed of a diffusion term and an advection term (see [8]). The advection term may be present in a reaction diffusion equation if the reactions are taking place in a flowing fluid.

Numerical experiments
In this section, we consider numerical experiments of the lagged diffusivity fixed point iteration for the solution on a square domain [0, 1] × [0, 1] of the problem (2)- (3).
Different functions for the nonlinearity factors σ (x, u) and g(x, u) have been considered. We are concerned with some simple steady-state problems (2)-(3) in which reaction kinetics and diffusion are coupled.
This coupling gives rise to reaction diffusion equations which, in a one-dimensional scalar case, can look like , t) is the concentration, −g(u) represents the kinetics and σ (u) is the diffusion coefficient, here taken to be concentration dependent. (When g(u) = −ru(1 − u/k), this equation is known as the Fisher-Kolmogoroff equation [14, Vol. I, p. 400]).
More specifically, we will consider reaction diffusion problems concerning the heat transfer caused by a chemical reaction (see, e.g. [17]) and the growth of spatially distributed communities (see, e.g. [2]).
The aim of this numerical study is to make a computational verification of the effectiveness of the lagged diffusivity procedure for solving nonlinear difference systems of the type (1) in an extended domain of convergence.
The source function s(x) is chosen in order to satisfy a prespecified exact solution u * = {u * (x i , y j )} of the nonlinear system (1), i, j = 1, . . . , N. The solution u * (x, y) is chosen equal to sin(πx) sin(πy).
In the following, we list the functions σ and g and how they are referred. We observe that for g1: g > 0, g > 0 and g > 0, for any value of u; for g2: g ≥ 0 and g > 0 when u ≥ 0; for g3: g ≥ 0, g ≥ 0 and g > 0 when u ≥ 0; and then, the functions g1, g2 and g3 satisfy the standard Assumptions on g for u ≥ 0. The lagged diffusivity fixed point iteration has been implemented in a Fortran code with machine precision 2.2 × 10 −16 .
This method stops when ε ν+1 ≤ ε, with ε = 10 −4 . An experiment has been carried out with different values of the threshold ε.
The starting vector of the lagged diffusivity fixed point iteration u (0) is the null vector (u (0) = 0) or the vector whose all components are equal to 1 (u (0) = e).
In the tables, ν * indicates the iteration of the lagged diffusivity fixed point method for which condition (33) is satisfied. The number ktot is the sum of the simplified Newton method's iterations and it is included in brackets; j k indicates the number of iterations of theAM method for the solution of the linear system that occurs at each iteration of the Newton method. This number of iterations of the AM method is chosen to be fixed at each simplified Newton iteration. For the details of the AM method used here, see [19].
Here err denotes the computed relative error in l 2 (R h ) norm, and diff indicates the last difference of iterations in l 2 (R h ) norm, While, with res and res0 we indicate the residual and the initial residual in the Euclidean norm: The term 4.81(−9) indicates 4.81 × 10 −9 .
From the numerical experiments, we can draw the conclusions as below.
• We observe that, since ε ν+1 decreases, for ν increasing, as Equation (34) and the lagged diffusivity iteration method stops at the iteration ν * when the criterium for ε ν * +1 in Equation (33) is satisfied, we have where we set ε 1 = 0.1 F(u (0) ) . Then In the experiments, we obtain ν * = log 2 ε 1 ε .    • We also observe that, generally, in the experiments the outer residual F(u (ν * ) ) has the same order of ε and the error in the discrete l 2 (R h ) norm has order of hε. • When in Equation (2) σ (u) = σ 2 and g(u) = g3, the diffusion term and the reaction term are small (Degenerate Equation). If the initial iterate u (0) is zero (Table 1), we have a nonmonotone decreasing of { F(u (ν) ) } that produces a large number of iterations of the simplified Newton method; e.g. we have F(u (0) ) = 933.26 and F(u (1) ) = 8514.10. The suggestion is to change the initial vector; the same case with u (0) = e produces better results (see Table 2) and a monotone decreasing of { F(u (ν) ) }. • When in Equation (2) σ (u) = σ 1 or σ (u) = σ 3 and g(u) = g3, the diffusion term dominates the reaction term (Equation of Diffusion-type). The matrix A(u) with σ (u) = σ 3 (g(u) = g3) is more ill conditioned than A(u) with σ (u) = σ 1. When in Equation (2) σ (u) = σ 1 or σ (u) = σ 2 and g(u) = g2, the slope of σ (u) is quasi-flat and the slope of g(u) is steep. When in Equation (2) σ (u) = σ 1 or σ (u) = σ 2 and g(u) = g1, the slope of σ (u) is quasi-flat and the slope of g(u) is quasi-flat with G (u) large. When in Equation (2) σ (u) = σ 3 and g(u) = g1 or g(u) = g2, the diffusion term is comparable to the reaction term. The numerical results on these test problems (Tables 1 and 2) confirm the statements 1 and 2 at the beginning of the previous section. • In the experiment in Table 3 the function σ is chosen dependent on u and on x and y. The range of values of this function σ is larger than that of σ 1 and a comparison with Tables 1 and 2 confirms that we could expect a larger number of the simplified Newton iterations.