A nonlinearity lagging for the solution of nonlinear steady state reaction diffusion problems

This paper concerns with the analysis of the iterative procedure for the solution of a nonlinear reaction diffusion equation at the steady state in a two dimensional bounded domain supplemented by suitable boundary conditions. This procedure, called Lagged Diffusivity Functional Iteration (LDFI)–procedure, computes the solution by “lagging” the diffusion term. A model problem is considered and a finite difference discretization for that model problem is described. Furthermore, properties of the finite difference operator are proved. Then, sufficient conditions for the convergence of the LDFI– procedure are given. At each stage of the LDFI–procedure a weakly nonlinear algebraic system has to be solved and the simplified Newton–Arithmetic Mean (Newton–AM) method is used. This method is particularly well suited for implementation on parallel computers. Numerical studies show the efficiency, for different test functions, of the LDFI–procedure combined with the simplified Newton–AM method. Better results are obtained when in the reaction diffusion equation also a convection term is present.


Introduction
We consider a nonlinear steady state reaction diffusion equation where the diffusion coefficient and the rate of change due to a reaction depend on the solution. These two terms are denoted with σ and −g, respectively. When we use a finite difference discretization, this elliptic equation supplemented by a Dirichlet boundary condition, can be transcribed into a nonlinear system of algebraic equations. We wish to compute a solution of this system of nonlinear equations with a common iterative procedure in which the nonlinear term, corresponding to the discretization of the diffusivity σ, may be evaluated at the previous iteration (see [30]). This approach of nonlinearity lagging in the diffusivity term is denoted as Lagged Diffusivity Functional Iteration (LDFI) or Lagged Diffusivity Fixed Point iteration. In Section 2, the LDFI-procedure is stated: the nonlinear difference system is solved via a sequence of systems of weakly nonlinear difference equations where only the term corresponding to g is nonlinear. Thus, the iterates of the LDFI-procedure are the approximate solutions of the weakly nonlinear systems computed with an inner iterative solver; a criterion for acceptability of these approximate solutions is given. A stopping rule for the LDFI-procedure is also given. Since, a purpose here is to re-examine the LDFI-procedure for solving the system of nonlinear difference equations of elliptic type in the context of Parallel Computing, a simplified version of the Newton-Arithmetic Mean (Newton-AM) method ([9]) is the inner iterative solver used for the solution of the weakly nonlinear systems. This is a two-stage iterative method and is particularly suited for implementation on parallel computers ([7], [8]; see also [1], [2], [3], [32]). In Section 6 a description of the simplified Newton-AM method and a general result on the convergence are reported. The purpose of Sections 4 and 5 is to analyse the convergence of the LDFIprocedure, considering the essential features of the system of nonlinear equations generated by a finite difference discretization of a reaction diffusion model problem described in Section 3. It is important to define for this model problem the set, called B, of all grid functions defined on the discretized domain which contains the solutions of the weakly nonlinear systems and all iterates of the LDFI-procedure.
For example, the model problem studied in [21] allows to present a helpful paradigm for proving the convergence of the LDFI-procedure. In Section 4, we summarize some properties of finite difference operators defined in B that also imply the uniform monotonicity of the nonlinear mapping which defines the nonlinear difference system. In Section 5, the convergence, to a solution of the original nonlinear system, of the sequence of iterates generated with the LDFI-procedure is proved under mild and reasonable assumptions imposed on the diffusivity σ and on function g using well known standard techniques. In the section of the numerical experiments (Section 7) the behaviour of the inner-outer iterations of the procedure is examined. The effectiveness of the LDFI-procedure combined with the simplified Newton-AM method is highlighted, especially, for reaction diffusion problems where also the convection term is present.

The lagged diffusivity functional iteration procedure
Consider a strongly nonlinear system of algebraic equations of the form where u = (u 1 , u 2 , ..., u n ) T is a vector in R n , A(u) is a large n×n nonsingular matrix with a sparse structure, and G(u) is a continuously differentiable diagonal mapping, i.e., a nonlinear mapping whose i-th component G i is a function of only the i-th variable u i for i = 1, ..., n. s is a vector of n components independent of u. We assume that this system has a solution u * .
For solving system (1) the easiest and maybe the most common method is to lag part of the nonlinear terms in (1) generating an iterative procedure denoted as Lagged Diffusivity Functional Iteration. With this iterative procedure the nonlinear system (1) can be solved via a sequence of systems of weakly nonlinear equations where the nonlinear term is G. Specifically, given a sequence of positive numbers {ε ν } such that ε ν → 0 as ν → ∞ and an initial estimate u (0) of the solution u * of the system (1), we generate a sequence of iterates {u (ν) }, ν = 0, 1, 2, ..., with the following rule for the transition from a current iteration u (ν) to the new iterate u (ν+1) : • Find an approximate solution u (ν+1) of the nonlinear system with the criterion for acceptability of the solution Then, the LDFI-procedure is composed by a nonlinear outer iteration that generates the sequence {u (ν) } and by an inner iterative solver of the weakly nonlinear system (2). This solver must be particularly well suited for implementation on parallel computers. The termination criterion for the outer iteration is provided by the following two inequalities where τ 1 and τ 2 are prespecified tolerances ([10, Theor. 3]). For the inner iterative solver, we use the following rule ε ν+1 = 0.5ε ν , ν = 1, 2, ...

A model problem
Consider a nonlinear steady state reaction diffusion convection equation of the form where ϕ = ϕ(x) is the density function at the point x of a diffusion medium Ω, σ = σ(x, ϕ) > 0 is the diffusion coefficient or diffusivity and is dependent on the solution ϕ, α = α(x) ≥ 0 is the absorption term,ṽ is the velocity vector, −g(x, ϕ) is the rate of change due to a reaction and s(x) is the source term.
In the equation (5) the convection termṽ · ∇ϕ has been taken into account; however, we will consider only convection-not dominated problems. We consider that equation (5) is subject by homogeneous Dirichlet boundary condition on the contour Γ of Ω: The functions σ(x, ϕ), α(x), s(x) and g(x, ϕ) are assumed to satisfy the following "smoothness" conditions: (i) the functions σ(x, ϕ) and g(x, ϕ) are continuously differentiable in x and continuous in ϕ; the functions α(x) and s(x) (the "source term") are continuous in x; (ii) there exist two positive constants σ min and σ max such that uniformly in x and ϕ; in addition, α(x) ≥ 0; (iii) for fixed x ∈ Ω, the function σ(x, ϕ) satisfies Lipschitz condition in ϕ with constant Λ (uniformly in x), Λ > 0; (iv) for a fixed x ∈ Ω, the function g(x, ϕ) is a uniformly monotone mapping ( [26, p. 141]) in ϕ with constant c > 0 (uniformly in x) and is continuously differentiable in ϕ.
We assume that the problem (5)-(6) has an isolated solution.
There exist various techniques for discretizing the problem (5)-(6). Using the Taylor series approach, equation (5) will be solved with the following standard finite difference scheme.
We consider Ω a rectangular domain (x ≡ (x, y) T ) with boundary Γ and we superimpose on Ω ∪ Γ a grid of points Ω h ∪ Γ h ; the set of the internal points Ω h of the grid are the mesh points (x i , y j ), for i = 1, ..., N and j = 1, ..., M , with uniform mesh size h along x and y directions respectively, i.e. x i+1 = x i + h and y j+1 = y j + h for i = 0, ..., N , j = 0, ..., M . Furthermore, at the mesh points of Ω ∪ Γ, (x i , y j ), for i = 0, ..., N + 1 and j = 0, ..., M + 1, the solution ϕ(x i , y j ) is approximated by a grid function u ij defined on Ω h ∪ Γ h and vanishing on Γ h . In order to approximate partial derivatives in (5) we shall make use of difference quotients of grid functions. The forward, backward and centered difference quotients with respect to x and to y of the grid function u ij at the mesh point (x i , y j ), are, respectively: while the centered second difference quotient with respect to x and to y can be written Providing a discretization error O(h 2 ), the finite difference approximation of (5) in (x i , y j ) is given by where the coefficients are: for i = 1, ..., N and j = 1, ..., M We denote the mesh points P k = (x i , y j ), (i = 1, ..., N, j = 1, ..., M ) and we order the points P k in lexicographic order: k = (j − 1) × N + i. We set n = N × M , and we denote the vector solution u whose components are the values of the grid function at the internal mesh points u = (u 1 , ..., u n ) T ≡ (u 11 , ..., u N 1 , u 12 , ..., u N 2 , ..., u 1M , ..., u N M ) T .
Then, formula (7) for i = 1, ..., N , j = 1, ..., M , can be written as formula where the matrix A(u) of order n has a block tridiagonal form.
The M diagonal blocks are tridiagonal matrices of order N with diagonal elements a kk (u) = D ij +D ij , the sub-and super-diagonal elements are In the following, we may consider the matrix A(u) as where A 1 (u) andÃ are the block tridiagonal matrices containing the elements {B ij , L ij , D ij , R ij , T ij } and {B ij ,L ij ,R ij ,T ij } respectively, while the matrix D is a diagonal matrix whose diagonal entries are {D ij }. Furthermore, s ∈ R n is a vector whose components are the values of the source term s(x, y) at the mesh points; the nonlinear mapping G(u) has components G k (u) = g(x i , y j , u k ), i = 1, ..., N , j = 1, ..., M and k = (j − 1)×N +i. We observe, that G k (u), the k-th component of G(u), with respect to the variable u, depends of only the k-th component u k , for k = 1, ..., n; in this case G is a diagonal mapping. We observe that the right hand side of (9) is the null vector since we have the homogeneous Dirichlet condition (6) in Γ. For grid functions {u ij } and {v ij } of this type the discrete l 2 (Ω h ) inner product and norm are defined by the formulas We say that the grid functions {u ij } defined on Ω h ∪ Γ h and vanishing on Γ h satisfy Property A if they are uniformly bounded and have uniformly bounded backward difference quotients ∇ x u ij and ∇ y u ij at each mesh point The set of all grid functions {u ij } which satisfy Property A is denoted by B. Thus, B is the set of grid functions {u ij } for which there exist some positive constants ρ and β such that The constant ρ is independent of h; also the constant β is independent of h but it depends on G(u) + s h . We assume that the system (9) has at least one solution [21], where a proof of the existence of such a solution u * of (9) in B is given; also a condition for which u * is unique in B has been obtained)

Some properties of finite difference operators
In the following we summarize some properties of finite difference operators whenṽ = 0 and α(x, y) = 0 in (5).
Here we denote Proof. We prove formula (13). We have 1 then, since v 0j = 0 for (6), the expression of the right hand side becomes and by u N +1j = 0, we have formula (13). Similarly, we obtain formula (14). ♯ We remark that from the right hand side of (13) and (14) we can swap u ij with v ij and we obtain Then, we have the following expression for the discrete l 2 (Ω h ) inner product of the vectors A(w)u and v where the n × n matrix A(w) is the one in (9), replacing u with w: Proof. Suppose the coefficients in (8), L ij , R ij , B ij and T ij , are functions of the grid function w ij , we have, Using formulae (13) and (14) and keeping into account of the inner product in l 2 (Ω h ), we have Then we have formula (15).

♯
We remark that while the grid function {u ij } is defined on the entire mesh region Ω h ∪ Γ h , the vector u ∈ R n represents the grid function {u ij } defined only on the interior mesh points Ω h , i = 1, ..., N , j = 1, ..., M . Moreover, we observe that formula (15) implies (9) and let A(w) the matrix n × n in (9) with u replaced by the vector w.
Then, if u, v and w belong to B, we have the following inequality where Λ > 0 is the Lipschitz constant of condition (iii), β > 0 is a constant for which |∇ x v ij | ≤ β and |∇ y v ij | ≤ β, and φ is an arbitrary positive number.
Proof. By using formula (15) in Lemma 2, we can write The assumption (iii) implies that, for given grid functions {u ij }, {w ij } belonging to Ω h ∪Γ h there exists a positive constant Λ such that for all i = 1, ..., N +1 and j = 1, ..., M + 1 The constant Λ is independent of h. Furthermore, Property A assures that there exists a constant β > 0 such that inequality (12) holds for all i = 1, ..., N +1 and j = 1, ..., M +1 and all grid function {v ij } belonging to B. The constant β is independent of h. Now, if we apply inequalities (18) and (12) into the expression in (17) we obtain Since the grid functions belonging to B are bounded and are equal to zero at the points of the boundary, we obtain The last expression can be written Using a well known technical trick, we have and we obtain formula (16). ♯ As consequence of Lemmas 1, 2 and 3, it is possible to prove that the mapping F (u) is uniformly monotone in B if the condition is satisfied (see [21] (and [12])). Thus, from Hadamard Theorem ( [17]), the nonlinear system (9) has a unique solution (e.g. [26, p. 143]). Note that the two hypotheses that F (u) is Lipschitz-continuous and uniformly monotone on R n are sufficient to prove that a solution of (1) exists and is unique; besides, it is possible to construct an iterative procedure that can guarantee a global convergence to the solution of (1) ( [20]).

Convergence of the LDFI-procedure
We will now investigate the solvability of the system of nonlinear difference equations (9) by applying the LDFI-procedure whenṽ = 0 and α(x, y) = 0 in (5). We will show that under the mild and reasonable restrictions (i)-(iv) imposed on the functions σ(x, ϕ) and g(x, ϕ) the problem (5)-(6) 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 ϕ. Specifically, if u (ν) is an estimate of the solution u * of (9), we will determine a new estimate of u * by solving the weakly nonlinear system (2) An approximate solution of the weakly nonlinear system (20) is computed by the simplified Newton-AM 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 ν → ∞. Here, · is the Euclidean norm. If such suitable solution u (ν+1) is found, we say that the algorithm does not break down.
Let u (0) ∈ B be arbitrary and let u (ν+1) be the solution of F ν (u) = 0 satisfying the condition (21) with F ν (u) as in (20). If all the vectors {u (ν) } belong to B and satisfy Property A with (22) instead of (12), then the sequence {u (ν) } converges to u * .
Proof. First we consider the case of α(x, y) = 0 andṽ = 0 for the problem (5)-(6). The solution u * in B of (9) satisfies the equation and the iterate u (ν+1) satisfies the equation where the discrete l 2 (Ω h ) norm of the residual F ν (u (ν+1) ) satisfies the inequality (21). Taking into account of the identity for all grid functions u, v and w belonging to B, we can write

Thus, we have
Using (15) and assumption (ii), we can write Assumption (iv) on g implies that, for all grid functions {u ij } and {v ij } belonging to B, there exists a positive constant c such that for all i = 1, ..., N and j = 1, ..., M . The constant c is independent of h. Thus for the discrete l 2 (Ω h ) inner product (10) we have the inequality Using Lemma 3 (see formula (19)) and taking into account of the assumption (iii) and the fact that, by Property A, the backward difference quotients |∇ x u (ν+1) ij | and |∇ y u (ν+1) ij | are bounded by inequalities (22), we can write It now follows from (23) that and from (24), (25) and (26) that where φ is a yet an undetermined positive number. Choosing Since the grid function {u (ν+1) ij } belongs to B, we may assume that Thus from (27) we have the inequality where γ = Λ 2 (β + ε ν+1 β 0 ) 2 2σ min c , and a = 2ρ/c. Now, as observed in [21], if there exists an integer ν 0 such that γ < 1 for all ν ≥ ν 0 , we can write (28) as Therefore, the sequence {u (ν) } of approximate solutions converges to the solution u * of the system (9). ♯ For sake of completeness, it easy to show (see [12]) that we have the convergence of {u (ν) } to the solution u * of the system (9) also in the cases α(x, y) > 0 andṽ = 0 for the problem (5)-(6). Indeed, since A(u) = A 1 (u) +Ã +D,Ã is the skew-symmetric part of A(u); then we have Then formula (23) becomes Furthermore setting α min = min (x,y)∈Ω α(x, y), we have then, at the left hand side of inequality (27) we have to add the term α min u * − u (ν+1) 2 h and in (28) the parameter γ becomes γ = Λ 2 (β + ε ν+1 β 0 ) 2 2σ min (α min + c) .

Solution of the weakly nonlinear system
In order to define the inner iterative solver for the nonlinear system (2) (or (20)), setting w (0) = u (ν) , the simplified-Newton method finds the solution for k = 0, 1, ..., where the matrix C ν is the Jacobian matrix of F ν evaluated at the point w (0) , i.e., C ν = F ′ ν (w (0) ) = F ′ ν (u (ν) ) and Denoting with G ′ (u) the Jacobian matrix of G(u) that has expression and taking into account the expression of C ν = A(u (ν) )+G ′ (w (0) ) = A(u (ν) )+ G ′ (u (ν) ) and the expression of F ν (w (k) ), formulae (29)- (30) are rewritten in such a way that the vector w (k+1) is the solution of the linear system for k = 0, 1, .... The system (31) is solved by the block version of the Arithmetic Mean method introduced in [29]. We remind that the matrix A(u) in (9) is a block tridiagonal matrix where the number of diagonal blocks is M .
Thus the matrix C ν has the following form: . . .
Here, {j k } denotes a sequence of positive integers. The loop over j denotes the Arithmetic Mean (AM) method. The description of the implementation and an evaluation of the effective performance of the Arithmetic Mean method on different parallel architectures are reported in the papers [29], [13], [14], [15], [16].
Letũ be a solution of the system (2). For any vector u and u (ν) belonging to an open neighbourhood K ofũ, we consider the following Standard Assumptions: • 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.
• The matrix A(u (ν) ) is irreducibly diagonally dominant and has positive diagonal entries and nonpositive off-diagonal entries for all the mesh spacings sufficiently small and for all the iterates u (ν) ∈ K.
• G(u) is a continuously differentiable diagonal mapping on R n with G ′ (u) ≥ 0 for all u ∈ K.
We report a general result on the convergence of the simplified Newton-AM method when the Standard Assumptions are satisfied. First we should define the matrix (ρ ≥ 0) and the iteration matrix and we observe that H ν = I − M −1 ν C ν . Theorem 2. Suppose the system (2) F ν (u) = 0 has a solutionũ; assume that Standard Assumptions hold for u belonging to an open neighbourhood K ofũ and that (33), i.e., are two splittings of the matrix C ν = F ′ ν (u (ν) ), u (ν) ∈ K, with the matrix H ν in (36) convergent. Then, for any j k ≥ 1, the solutionũ is an attraction point of the simplified Newton-AM iteration {w (k) } defined in (34).

Numerical studies
In this section we consider a numerical experimentation of the LDFI method for the solution on a rectangular domain of the model problem (5) with homogeneous Dirichlet boundary condition (6). Different functions for the nonlinearity factors σ(x, ϕ) and g(x, ϕ) and for α(x) have been considered. The source function s(x) is chosen in order to satisfy a prespecified exact solution u * = ϕ(x j , y j ) of the nonlinear system (1), i = 1, ..., N , j = 1, ..., M ; different choices for ϕ(x, y) are examined. In the following we list the involved functions and how they are referred. The functions σ are dependent on ϕ and are: σ1 : σ(ϕ) = 0.5 + 0.5ϕ, σ2 : σ(ϕ) = 0.02 + 0.5ϕ 2 , σ3 : σ(ϕ) = 1/(0.02 + 0.5ϕ).
Set Set The LDFI-procedure has been implemented in a Fortran code with machine precision 2.2 × 10 −16 .
In the experiments, we consider as stopping criterium for LDFI-procedure the satisfaction of both the inequalities (4) The approximate solution computed, at each iteration of LDFI-procedure, by the simplified Newton method satisfies the stopping rule with ε 1 = 0.1 F (u (0) ) and ε ν+1 = min{0.5 ε ν , ε}, ν = 1, 2, .... The threshold ε is chosen 10 −5 , 10 −3 or 10 −2 . In Tables 3-7, ε is chosen equal to 10 −5 . The starting vector of the LDFI-procedure u (0) is the vector whose all components are equal to 1. In all the experiments we have N = M . In the tables, it indicates the number of iterations of the LDFI-procedure. The number ktot, the sum of the simplified Newton method's iterations, is expressed in brackets.
Here err denotes the computed relative error in the Euclidean norm, i.e.
with res and res0 we indicate the residual and the initial residual in the Euclidean norm: and diff indicates the last difference of iterations The term 7.60(−10) indicates 7.60 × 10 −10 . We indicate with j k the number of iterations of the Arithmetic Mean method for the solution of the system (31) ) + s, that occurs at each iteration k of the simplified Newton method. At each iteration j, j = 1, ..., j k , of the Arithmetic Mean method, M − 1 independent 2 × 2 block linear systems of order 2N have to be solved; the block Gaussian elimination method is used and, the parameter ρ in the AM method is chosen equal to zero (see [29]).

Conclusions
In this paper we have analysed the LDFI-procedure combined with the simplified Newton-AM method for the solution of finite difference nonlinear systems. For the convergence of the LDFI-procedure we have considered: • a model problem where smoothness conditions on the functions involved in the equation of the model are assumed; • the grid functions, i.e. discrete approximations of the solution of the model problem, satisfy Property A (i.e., they are uniformly bounded and have uniformly bounded backward difference quotients); • the solution of the weakly nonlinear difference system that occurs at each iteration of the LDFI-procedure, is solved inexactly by a convergent iterative solver; • the theorem of convergence of the LDFI-procedure is proved with standard techniques.
From the numerical experiments the following conclusions can be drawn: • from Tables 1 and 2, we can observe that when the values of the function σ(ϕ) increase (σ3 has larger values than σ1 and σ2 for ϕ ∈ [0, 1]), then the total number of the simplified Newton iterations increases; • from Tables 1, 2 and 3, we observe that when the values of the function g(ϕ) are rapidly increasing for 0 ≤ ϕ ≤ 1 or the values of the function α(x, y) are large, then the diagonal of the matrix C ν becomes more dominant; it implies a reduction of the total number of the simplified Newton iterations; • from Tables 1 and 2, we remark that there is no appreciable reduction of the total number of the simplified Newton iterations when we solve the weakly nonlinear system with a looser accuracy than that imposed on the LDFI-procedure. Thus, in the strategy of choice of the parameters in criteria (4) and (21), these experiments suggest that the parameter τ 2 must have approximately the same value of the threshold ε.
• from Tables 5 and 7, we remark that the LDFI-procedure combined with the simplified Newton-AM method gives better results (in terms either of total number of simplified Newton iterations or of the number of LDFI-procedure iterations) when the coefficients ofṽ increase, i.e., when the deviation from asymmetry 2 of the matrix C ν increases. This is a peculiar feature of the Arithmetic Mean linear solver ( [29]), especially when it is implemented as inner solver in a two iteration levels procedure, such as the Newton-AM method ([9], [11]).