A nonlinearity lagging method for non-steady diffusion equations with nonlinear convection terms

We analyze an iterative procedure for solving nonlinear algebraic systems arising from the discretization of nonlinear, non-steady reaction-convection-diffusion equations with non-constant (and, in general, nonlinear) velocity terms. The basic idea underlying the procedure consists in lagging the diffusion and the velocity terms of the discretized system, which is thus partly linearized. After analyzing the discretized system and proving some results on the monotonicity of the operators and on the uniqueness of the solution, we prove sufficient conditions that ensure the convergence of this lagged method. We also describe the inner iteration and show how the weakly nonlinear systems arising at each lagged iteration can be solved efficiently. Finally, we analyze numerically the entire solution process by several numerical experiments.


Introduction
In many important problems and applications, it is often necessary to solve systems of nonlinear algebraic equations focus our analysis on the more general non-steady case, but we provide a few remarks on the stationary case as well.
Then, on a more operational point of view, we illustrate how to choose starting vectors and tolerances of all the iterative procedures in an efficient way. This is of fundamental importance in order to achieve a fast implementation, since we have to carry out three nested iterative procedures at each time level. Indeed, the lagging procedure transforms the nonlinear system into a sequence of weakly nonlinear systems, each of which is solved, in this paper, by the simplified Newton's method, which, in turn, requires solving a linear system at each iteration, for example by an iterative linear solver. Finally, we introduce several numerical experiments to study the behavior of the algorithm as the velocity term, the inner linear solver or other parameters of the procedure are changed.
The paper is structured as follows. In Section 2, we introduce some assumptions on the smoothness of the involved functions and describe the discretization of the differential problem, illustrating how F (u) is made.
After some initial remarks, in Section 3, we then prove some lemmas needed in the successive analysis, devoting particular attention to the analysis of terms dependent onṽ. We then characterize the existence of at least a solution to F (u) = 0. Lastly, we prove that, under some hypotheses, F (u) is monotone and that the solution of F (u) = 0 is unique.
Then, in Section 4, we describe the LDM and prove the convergence of the procedure. We also provide a few remarks on the application of the method to stationary problems.
In Section 5, we describe the solution process, describing how we solve the weakly nonlinear systems arising at each iteration of the LDM and how starting vectors and tolerances of all the iterative procedures can be chosen efficiently. The resulting algorithm is reported in Appendix.
Finally, Section 6 is devoted to the numerical experiments and Section 7 concludes this work.
Let (x, y) ∈ and letû be in a neighborhood of a solution. We assume that the following smoothness assumptions hold.

Space discretization
For simplicity, let be a 2D bounded rectangular domain. Let us superimpose to¯ a grid¯ h = h ∪ ∂ h of mesh points (x i , y j ), i = 0, . . . , N + 1, j = 0, . . . , M + 1, defined by x i+1 = x i + h, i = 0, . . . , N y j +1 = y j + h j = 0, . . . , M, where h is the mesh size along x and y (considered uniform). Space discretization is performed by approximating the space derivatives in (1) by finite difference methods. Committing a discretization error O(h 2 ) and denoting by u i,j (t) the grid function that approximates the solution u(x i , y j , t) at the mesh points (x i , y j ) of¯ h , i = 0, . . . , N + 1, j = 0, . . . , M + 1, the right-hand side of (1) can be written as x [σ (x i , y j , u i,j (t))∇ x u i,j (t)] + y [σ (x i , y j , u i,j (t))∇ y u i,j (t)] −ṽ 1 (x i , y j , u i,j (t), t)δ x u i,j (t) −ṽ 2 (x i , y j , u i,j (t), t)δ y u i,j (t) −α(x i , y j )u i,j (t) − g(x i , y j , u i,j (t)) + s(x i , y j , t), where we denote forward finite differences (in x and in y respectively) by x , y , backward finite differences by ∇ x , ∇ y and central finite differences by δ x , δ y . For compactness of notation, in the following we also denote σ (x i , y j , u i,j (t)) by σ (u i,j ),ṽ 1 (x i , y j , u i,j (t), t) byṽ 1 (u i,j , t) andṽ 2 (x i , y j , u i,j (t), t) byṽ 2 (u i,j , t).
We then write the first five terms of (3) as a matrix-vector product −A(u(t))u(t) plus a vector b(u(t)) ∈ R μ which contains terms coming from the Dirichlet boundary conditions. In particular, A(u(t)) ∈ R μ×μ is a block tridiagonal matrix. The M diagonal blocks are tridiagonal matrices of order N of elements {−L i,j ,D i,j , −R i,j }. Its lower-and upper-diagonal blocks are instead diagonal and are made of the elements {−B i,j } and {−T i,j }, respectively. Moreover, it is irreducible and, if holds, it is irreducibly diagonally dominant [20, p.23] with positive diagonal elements and non positive off-diagonal elements. In this case, A(u(t)) is, thus, a non-singular M-matrix [20, p.91]. The vector b(u(t)) can instead be easily obtained by applying the boundary conditions to the terms depending on u(t) in (3).
We thus write (3) for i = 1, . . . , N, j = 1, . . . , M in the form and (1) is replaced by the system of ordinary differential equations

Time discretization
Passing to time discretization, let us introduce a time spacing t and define a series of time levels t n = n t +t 0 , n = 0, 1, . . .. Calling s n := s(t n ) and denoting by u n the approximation of u(t n ) solution of (6) at t = t n , we write the well-known θ -method (e.g., see [6]) for n = 0, 1, . . . and with 0 ≤ θ ≤ 1. In particular, in the following, we consider an implicit time discretization (hence, we consider cases where θ = 0). The use of an explicit scheme would, indeed, easily require an extremely short time step. For instance, it is easy to verify that the problems later analyzed in the numerical experiments would need a time-step as small as t = 10 −6 for the explicit method to converge. This, of course, greatly reduces efficiency. The advantages of an implicit time discretization come, however, at the cost of a nonlinearity in the discretized system, which we will later handle by the lagging procedure. Let us then write more compactly the nonlinear system arising from the discretization. Calling I the μ × μ identity matrix, we collect the known terms in a vector w = w n ∈ R μ defined as At each time level n = 0, 1, . . ., the vector u n+1 is thus given by the solution of the nonlinear algebraic system where τ := θ t.
3 Uniform monotonicity of F(u) and uniqueness of the solution of F(u) = 0 Thus, in (8), we obtained the nonlinear algebraic system F (u) = 0 that we aim to solve by the lagged diffusivity method. In the first part of this section, we assume that a solution exists and introduce some preliminary lemmas used in the following. We then better characterize the existence of at least one solution by exploiting the smoothness assumptions at the beginning of Section 2. Lastly, we study the monotonicity of the operator F (u) and the uniqueness of the solution of F (u) = 0.

Initial remarks
It is convenient to split the matrix A(u) in Finally,Ã(u) is the block tridiagonal matrix whose lower-and upper-diagonal blocks are diagonal matrices of elements We then further split A 1 (u) andÃ(u) in where A x 1 (u) andÃ x (u) are block diagonal matrices whose diagonal blocks are tridiagonal with row elements is the block tridiagonal matrix where the sub-, main and super-diagonal blocks are diagonal matrices of diagonal elements −B i,j , D y i,j and −T i,j respectively, with D y i,j = B i,j + T i,j . Finally,Ã y (u) has the same structure of A y 1 (u) and the diagonal matrices have diagonal elements −B i,j , 0 and −T i,j respectively.
It is useful to split also b(u) in a similar way. We thus write where b 1 (u) andb(u) contain, respectively, terms dependent on the diffusivity and on the velocity term. We then further split these two vectors in where b x 1 (u) andb Finally, given two grid functions u, v in¯ h , in the following we make use also of the l 2 ( h ) discrete inner product and of its associated norm · h .

Preliminary lemmas onÃ (u)
Lemma 1 Let {u i,j }, {v i,j } and {w i,j } be three grid functions defined at mesh points (x i , y j ) of a grid¯ h , i = 0, ...N + 1, j = 0, ..., M + 1, and satisfying the Dirichlet boundary conditions u i,j = u 1 (x i , y j , t) ∀(x i , y j ) ∈ ∂ h and t > 0. Then: Proof Let us split the inner product in By definition of discrete l 2 ( h ) inner product and by the form ofÃ x (u) and ofb x (u), with simple algebraic passages, we find Then, collecting terms, by definition of central finite-difference quotients, we get We find a similar relationship by proceeding analogously for the term containing A y . Using these results in (9) and collecting terms, we finally prove the lemma. (9) can be obtained also when v i,j satisfies homogeneous Dirichlet boundary conditions for Ã (u) · w −b(u), v and for Ã (u) · v, v .

Corollary 1 Expressions analogous to
Proof The proofs follow Lemma 1 without relevant modifications. Regarding Ã (u)· v, v , however, we also need to consider that, since v is null on ∂ h , we have and analogous expressions for terms inṽ 2 (u i,j )v i,j ±1 .

Lemma 2
Let {u i,j } be a grid function defined at mesh points (x i , y j ) of a grid h , i = 0, . . . , N + 1, j = 0, . . . , M + 1, such that {u i,j } satisfies the Dirichlet boundary conditions u i,j = u 1 (x i , y j , t) ∀(x i , y j ) ∈ ∂ h and t > 0. Moreover, let the backward difference quotients ∇ x u i,j and ∇ y u i,j be bounded. Then, denoting bỹ v max the bound on the velocity variable (see point 6, Section 2), Proof By δ x u i,j = (∇ x u i+1,j + ∇ x u i,j )/2 (and similarly for δ y u i,j ), we can write By the boundedness ofṽ and by the inequality ab ≤ (a 2 + b 2 )/2 with a, b real numbers, we can further evaluate Changing sign to the right-hand side of the equation and by definition of u h norm, we finally prove the lemma.

Existence of at least one solution to the discretized system
We now introduce a sufficient condition for the existence of at least one solution to the discretized system (8). We notice that we can put ourselves in this condition by a suitable choice of the discretization parameters h and τ .

Theorem 1 Let {u i,j } be a grid function defined at mesh points (x i , y j ) of a grid
h , i = 0, . . . , N + 1, j = 0, . . . , M + 1, such that {u i,j } satisfies the Dirichlet boundary conditions u i,j = u 1 (x i , y j , t)∀(x i , y j ) ∈ ∂ h and t > 0. Moreover, let the backward difference quotients ∇ x u i,j and ∇ y u i,j be bounded and let σ min ,ṽ max and α min be the bounds on diffusivity, velocity and absorption terms defined at the beginning of Section 2. If h ≤ 2σ min /ṽ max and then there exists at least one solution to system (8). Furthermore, all solutions belong to a ball of radius ρ with Many of these terms can be easily evaluated. Indeed, u, u = u 2 h by definition, and Regarding the other terms, we instead rely on other evaluations. Thus, by [11, Lemma 1], we easily find Therefore, combining all these evaluations, Since h ≤ 2σ min /ṽ max by hypothesis, σ min −hṽ max /2 ≥ 0 and we can further evaluate from above by canceling out the term with backward quotients. We also notice that this condition on the step h of the space discretization implies that condition (4) is satisfied as well and A(u) is then an M-matrix.
Collecting u h , we finally have This evidently means that no solution to F (u) = 0 can lie outside {u | u h ≤ ρ}.
The existence of at least one solution follows then by [8,Lemma 4.3]. Another proof, based on the invariance of the degree of a mapping under a homotopy [12, p. 156], is provided in [9, Theorem 2.1], which considers that H (u, δ), u > 0 when The condition h < 2σ min /ṽ max , often implying h small, is evidently in competition with the condition 1 − τṽ max /h + τ α min + τ c > 0, which is harder to satisfy as h gets smaller. It is however important to notice that a small step in the space grid can be compensated by an equally small step in time. Thus, in line of principle, both conditions can always be satisfied by acting on the size of space and time grids.
Finally, it is also possible to provide a bound to the backward difference quotients of a solution of (8) at each mesh point (x i , y j ) ∈¯ by proceeding similarly to [9, Lemma 2.2]. We denote this bound by β. It is then convenient to define the ball B ρ,β to which all solutions must belong. Definition 1 Let u :¯ h → R be a grid function satisfying the Dirichlet boundary conditions on ∂ h for t > t 0 and let the hypotheses of Theorem 1 be satisfied. We say that u belongs to B ρ,β if

Monotonicity analysis
Before studying the monotonicity of F (u), we need a last lemma, which relies on the bound β to the finite difference quotients of the grid functions in B ρ,β .
Lemma 3 Let u, v and w be three grid functions belonging to B ρ,β . Then, Therefore, by the triangle inequality, Next, by definition of backward and central finite-difference quotients and by boundedness of backward finite-difference quotients (12) in B ρ,β , we can write Furthermore, by Lipschitz continuity of every component of the velocity vectorṽ, we have Identical evaluations apply, respectively, to δ y v i,j and to ṽ 2 (u i,j ) −ṽ 2 (w i,j ) . Therefore, combining these results, Since, for a, b real numbers, ab ≤ (a 2 + b 2 )/2, we finally get which, by definition of discrete l 2 ( h ) norm, becomes: then F (u) is uniformly monotone in B ρ,β and the solution u * of (8) is unique.
Proof To prove that F (u) is uniformly monotone in B ρ,β , we must show that there exists a positive scalar γ satisfying To this end, let us analyze (8), adding and subtracting A(u)v and rearranging terms, we can write By the splittings of A(u) and of A(v), the right-hand side of the previous equation becomes We can make evaluations analogous to those in Theorem 1 for most of these terms, obtaining: Proceeding as in Lemma 2 and by Lemma 3, we then have, respectively Lastly, proceeding as in [11,Theorem 1], we obtain where φ is an arbitrary positive parameter. Thus, considering all previous inequalities and collecting terms, we can write Since φ is an arbitrary, positive parameter, we now choose it in a suitable way. In particular, we choose which implies We here exploit that φ is positive by the hypothesis h < 2σ min /ṽ max , which implies 2σ min − hṽ max > 0. With this choice of φ, we get h . Thus, we have obtained an equation in the form of (15) and monotonicity holds if is larger than zero. This happens when proving the first part of the theorem. The uniqueness of the solution follows then directly. Indeed, suppose that two distinct solutions u * andû to the system F (u) = 0 exist in B ρ,β . We would have As a final remark, we also notice that the monotonicity conditions mirror those in Theorem 1. In particular, if the monotonicity conditions are satisfied, the conditions of Theorem 1 are satisfied as well. Then, Theorem 2 can also be seen as a sufficient condition for the existence and the uniqueness of the solution.

The lagged diffusivity method
The basic idea of the method consists in linearizing (8) by setting up an iterative procedure where diffusivity and velocity are lagged. In order to do this, chosen a starting vector u (0) , at the (ν + 1)-th iteration, we consider the diffusivity and the velocity terms dependent on the solution u at the ν-th iteration, u (ν) . This means that, at each lagged iteration, instead of having the nonlinear terms A(u) and b(u), we have A(u (ν) ) and b(u (ν) ). A weak nonlinearity is, however, still present due to the nonlinear mapping G(u).
Hence, we compute the new iterate u (ν+1) as the solution of the weakly nonlinear system Here, System (16) can be solved approximately by an iterative method. The lagged iterate is accepted when the residual satisfies a stopping criterion where · denotes the Euclidean norm and ν is a given tolerance such that ν → 0 for ν → ∞. In the following, we refer to (18) as lagged acceptability condition.
When it is satisfied, we find u (ν+1) and the outer iteration can be restarted: the computed u (ν+1) is used to evaluate diffusivity and velocity to find u (ν+2) by the new lagged iteration, and so on. In the following, we choose the solution at the previous time level as starting vector of the lagged diffusivity method, which is we set u (0) = u n to initialize the LDM iteration at the (n + 1)-th time level. At the first time level, we instead use the initial condition. Tolerances are then defined starting from the norm of the initial residual, F (u (0) ) . Indeed, we define 1 by multiplying the initial residual by a positive constant smaller than 1, e.g., and we build a sequence of ν+1 , ν = 1, 2, . . . such that ν → 0 for ν → ∞ simply by setting The lagged diffusivity procedure is then stopped when is satisfied, where¯ is a given tolerance. As a consequence of computing ν+1 by halving 1 ν times, condition (19) is satisfied after iterations, where · denotes the ceiling function.

Convergence analysis
We now analyze the convergence of the algorithm presented in the previous subsection. In this regard, let us first better characterize the solutions of the weakly nonlinear systems arising at each lagged iteration. In particular, if we solve (16) inexactly, the approximate solution u (ν+1) solves Assuming that the hypotheses of Theorem 1 are satisfied, let us analyze It is also interesting to notice that ρ ν+1 can be expressed as . This will be useful in the analysis of convergence and it implies that ρ ≤ ρ ν+1 for all ν = 0, 1, . . .. Proceeding similarly with regard to the bound on backward difference quotients, it is possible to write β ν+1 = β + ν+1 β 0 . We then have u (ν+1) ∈ B ρ ν+1 ,β ν+1 for ν = 0, 1, . . ..

Theorem 3
Let u * ∈ B ρ,β be the solution of the nonlinear system F (u) = 0 defined in (8) with A(u) non-singular and G(u) diagonal mapping. We assume that the smoothness conditions and the hypotheses of Theorem 2 are satisfied.
Proof Let us consider F (u * ) = 0 and the lagged acceptability condition (18), satisfied by u (ν+1) : Then, let us subtract the second equation from the first one, divide by τ and cancel out constant terms, obtaining At this point, we add and subtract the same term A(u * )u (ν+1) and rearrange terms, similarly to what done at the beginning of Theorem 2. If we then take the inner product of both sides with u * − u (ν+1) , we get Next, we consider the splittings of A(u) and of A(u * ). The right-hand side of the previous equation thus becomes We can now use the previously introduced lemmas and theorems to evaluate all these terms. Thus, we find that the previous expression is larger than or equal to where φ is an arbitrary, positive parameter. Choosing φ = 1 β ν+1 (2σ min − hṽ max ) (which is positive by the hypothesis h < 2σ min /ṽ max ), the first term of the previous expression vanishes and we get We now rewrite the previous inequality by evaluating from above the term on the left-hand side. Since the iterate u (ν+1) must satisfy condition (18) and since it must belong to B ρ ν+1 ,β ν+1 (implying u * − u (ν+1) ≤ 2ρ ν+1 since ρ ≤ ρ ν+1 ), we have For compactness, let us then collect the coefficient of the first term of the right-hand side of the previous inequality in Using the relation β ν+1 = β+ ν+1 β 0 and comparing this equation with the definition of γ in the monotonicity condition, it can be noticed that γ > 0 (which is satisfied by hypothesis) implies that also ζ is positive if holds true. Assumption (20) is nonetheless plausible since the tolerance ν → 0 for ν → ∞. Therefore, inequality (20) is certainly satisfied from a certain lagging iteration. Moreover, operatively, we could grant that (20) is satisfied also at the first lagging iterations by setting 0 sufficiently small. Thus assuming ζ > 0, we can divide both sides by ζ without changing sign to the inequality. Rearranging terms, we get We can then make further considerations on ζ . Indeed, again remembering assumption (20) and γ > 0, we have Since ν → 0 for ν → ∞, we also have So, since the inequality in (22) holds in strict sense, we can assume that there exists an integer ν 0 such that Thus, definingγ proceeding iteratively from a ν 0 such thatγ < 1, by (21), we get Thus, the previous theorem proves the convergence of the method when it is applied to solving the nonlinear systems arising from a finite-difference discretization of the differential problem (1). It is nonetheless important to highlight that the lagged diffusivity method is applied to the discretized system, irrespectively of the used discretization. Thus, it can be applied, in principle, also when other discretization techniques (i.e., finite element or finite volume discretizations) are used. However, it is important to remark that the presented proofs of convergence rely on monotonicity properties which FD operators satisfy.

Remarks on the stationary case
In case we are studying a steady reaction-convection-diffusion equation, discretization has to be performed only along space. Using the finite difference schemes employed in the non-steady case, space discretization leads to an expression similar to (5). The only difference is that we do not have time dependence. Therefore, after space discretization, we do not have a system of ordinary differential equations like (6), but we get directly the nonlinear algebraic system We can then introduce the lagging iteration: with F (u) as in (24), the new iterate u (ν+1) of the lagged iteration is given by the solution of the weakly nonlinear system If A(u) is nonsingular, system (25) can be solved approximately by a convergent iterative method and the lagged iterate is accepted when the residual satisfies a stopping criterion where ν is a given tolerance such that ν → 0 for ν → ∞.
Conditions of monotonicity and convergence in stationary case can be found proceeding as done in Theorems 2 and 3.

Solution process
In this section, we analyze more in detail the solution process. In this regard, we describe inner and outer solvers and how to set their initialization and stopping criteria efficiently.
As described in Section 4.1, at each LDM iteration, we must solve the weakly nonlinear system (16), which we do approximately by an iterative procedure that is stopped when the lagged acceptability condition (18) is satisfied. We choose to solve (16) by the simplified Newton's method [12, p. 182]. We are thus left with three solution levels at each time step: 1. we use the lagged diffusivity method to linearize the nonlinear algebraic systems (8); 2. we use the simplified Newton's method to solve the weakly nonlinear algebraic system (16) arising at each lagged iteration. The method is stopped when the lagged acceptability condition (18) is satisfied; 3. we use an iterative linear solver to solve the linear system arising at each simplified Newton's iteration. In the following, the tolerance of the linear solver is chosen so that (16) is solved by an inexact Newton's method [4,5].
Having already described the lagged diffusivity method in Section 4.1, let us proceed to the analysis of the solution of the weakly nonlinear system (16).
Then, replacing F ν (u (ν+1,k) ) by its definition in (16) and canceling out opposite terms, we find that at the (k + 1)-th simplified Newton iteration, k = 0, 1, . . ., u (ν+1,k+1) is the solution of the linear system Thus, now we have to solve a linear system at each Newton iteration. This can be done approximately by an iterative method. This helps increasing the computational efficiency, especially in case of systems of large dimensions. Let us then introduce a third superscript j . At the (j + 1)-th iteration of the linear solver, j = 0, 1, . . ., we have a residual r j +1 given by where u (j +1) ≡ u (ν+1,k+1,j +1) denotes the solution of (30) at the (j + 1)-th linear iteration of the (k + 1)-th Newton iteration of the (ν + 1)-th lagged iteration. In the following, the starting vector of the linear iterative solver is given by the solution of the previous simplified Newton iteration, i.e. u (ν+1,k+1,0) = u (ν+1,k) .
Choosing a suitable stopping condition of the linear solver, we can set up an inexact Newton procedure. In this regard, we choose a forcing termσ , 0 <σ < 1, and define the toleranceˆ (k+1) of the linear solver at the (k + 1)-th Newton iteration aŝ We then stop the linear solver when the norm of r j +1 , j = 0, 1, . . ., satisfies Denoting by j end the number of iterations needed to the linear solver for satisfying the stopping condition (33), we define u (ν+1,k+1) ≡ u (ν+1,k+1,j end ) . Finally, we provide a few operative remarks on the choice of the arbitrary parameters of the procedure. The choice of the prescribed tolerance¯ is connected to the desired accuracy, as we can observe that the value of the final residual is ultimately related to it. Regarding 0 andσ , values excessively small or large may reduce efficiency, but their choice is, nonetheless, not problematic. In all our experiments, we fixed 0 =σ = 0.1.

A correction of the initialization of starting vectors
Starting vectors and stopping criteria of the used iterative procedures can be summarized as follows: The starting vectors determine also the initialization of tolerances. For instance, at the first Newton iteration of the (ν + 1)-th lagged iteration, the tolerance of the linear solver isˆ (1) =σ F ν (u (ν+1,0) ) =σ F ν (u (ν) ) . We can also easily notice that all tolerances are initialized atσ F 0 (u (0) ) =σ F (u (0) ) = 1 at the first lagged iteration.
If we do not apply any correction, however, it can happen that the initial tolerance of the linear solver,ˆ (1) , become smaller thanσ ν+1 [11,Sec. 5.2]. Being ν+1 the tolerance of the simplified Newton's method,ˆ (1) would thus be smaller than what could possibly be required by an inexact Newton's method, leading to an increase in computational cost with no foreseeable improvement to accuracy.
Then, we apply a correction on the initialization of the tolerance of the linear solver. In this way, we avoid solving too exactly the linear system and we reduce computational cost and numerical difficulties. This can be easily done by imposinĝ σ ν+1 as the minimum tolerance of the linear solver. Following [11], we thus set

Numerical experiments
In this section, we solve various problems by a Fortran implementation of the LDM. We are concerned with how the method behaves with a nonlinear velocity term and we consider several possible choices ofṽ: with c 1 and c 2 real constants. These choices represent different possible situations: indeed,ṽ 1 depends only on u,ṽ 2 depends on u, x and y but it is constant in time, whileṽ 3 andṽ 4 depend on t as well. All these choices evidently satisfy the initial smoothness assumptions (for u finite).
We start our analysis by verifying the effectiveness of the lagged diffusivity method. In this regard, we choose, for example, the following test problems: Problem 1 represents a situation where σ min and g are quite large, while in Problem 2 we have the opposite situation. Also the monotonicity constant of g is smaller in this latter case. Thus, Problem 2 represents a situation which is potentially more critical for the monotonicity of F (u) and for convergence.
Finally, it appears suitable to introduce also a test problem motivated by situations which could occur in real-world applications. In this regard, consider, for instance, all those cases where velocity is produced by an external force F . This is, for example, the case of the Smoluchowski diffusion equation, which describes the flow of ions dissolved in a liquid in the presence of an electric field that pulls the ions in a given direction. In this case, the velocity can be written as the quotient between the force of the field and a term ζ (called viscous drag, which accounts for the friction) and may be nonlinear (see, e.g., [1]). Conceptually, similar situations may arise also in chemical-mechanical frameworks. For instance, consider all those processes where the material which is diffusing has a larger (or smaller) viscous drag and/or causes reactions that form substances characterized by a higher friction. The velocity produced by the same force, thus, may be nonlinear with concentration. Finally, we can find similar situations also in plastic compounding, which consists in mixing additives to a molten polymer. If we add a plasticizer, the viscosity of the polymer is reduced. Thus, if the polymer is mixed by applying a fixed force F , the velocity of the fluid is larger where the concentration of the plasticizer is higher. The opposite situation may arise, on the other hand, if additives that increase viscosity (e.g., fillers) are used.
Let us then build a test problem which can constitute a simplified model of these situations. Assume, that a specie b is diffusing in a medium a. The starting concentration of b in the domain is zero and, as b diffuses, the medium is mixed, with a velocity produced by an external force F . Supposing that b has a large viscosity, the parameter accounting for the viscous drag, η, will be larger where the concentration of b is higher. For instance, assume that η increases quadratically as the concentration of b increases (in any case, it is easy to modify the problem and consider different dependence laws). Then, calling η 0 the value of η where the concentration of b is zero, the velocity may be expressed as with κ > 0 and with F 1 and F 2 representing the components of the force along the axes x and y, respectively. Notice that we have here performed some simplifying assumptions, such as that F 1 , F 2 are constant in the points of the domain. If this does not hold true (like, for instance, in case of a mechanical mixing, where we would likely have turbulences and we may stir the medium along a circumference), we can nonetheless adapt the formulation by considering F 1 and F 2 as functions of x, y.
Then, we further assume that a reaction between b and the diffusion medium occurs and that its rate is proportional to the concentration of b itself. Finally, we may consider a linear diffusivity, in accordance with the formulation of common mass diffusivities and so to analyze also a case where only the velocity term is nonlinear.
Fixing a solution (for instance, the same as that of Problem 1) and choosing some parameters, let us then solve the test problem For all problems, we choose the domain to be the square [0, 1] × [0, 1], which we discretize by a uniform grid of N × N points. We also set the initial time t 0 = 0, the final time t f = 1 and θ = 0.5. When not otherwise specified, we set N = 250 and t = 0.1.

Verification of the method and analysis of the linear solvers
Let us solve the test problems considering all the forms ofṽ introduced in (35) with c 1 = c 2 = 1, α = 0 and different inner linear solvers: the Arithmetic Mean (AM) method [17] with ω = 1, the BiConjugate Gradient-stabilized method with parameter l (BiCGstab(l)) [19] and the GMRES method [18]. The implementation of the GMRES follows [7, p. 45], which uses Givens rotations and re-orthogonalization. We choose to consider these solvers because, beingṽ = 0, the coefficient matrix of the linear systems arising at each Newton iteration cannot be symmetric, thus preventing the use of a conjugate gradient method. Depending on the problem, the non-symmetry can then be more or less significant, and it is thus interesting to compare the AM method, which is known to be well-suited for solving strongly non-symmetric systems [16], and Krylov solvers such as BiCGstab(l) and GMRES. Moreover, since we here aim especially at analyzing the behavior of the algorithm, we do not apply any preconditioner. Preconditioning techniques can, nonetheless, be easily employed if higher efficiency is required.
We stop the lagging iteration when ν+1 ≤¯ = 10 −4 , while starting vectors and stopping criteria of all the iterative procedures are chosen as described in Section 5. We also set a maximum number of Newton iterations k max = 500 for each lagging iteration and a maximum number of linear iterations j max = 10,000 for each Newton iteration.
The results are reported in Tables 1, 2 and 3 and are referred to the last time level, except for t tot , which denotes the whole time required to compute the solution from t 0 to t f . By ν end , k end and j end we denote, respectively, the total number of lagged, Newton and linear solver iterations at the last time level. Finally, res 0 denotes the Euclidean norm of the initial residual, res denotes the Euclidean norm of the final residual and err h and err 2 denote the global error in l 2 ( h ) and in relative Euclidean norm respectively. We are able to compute the correct solution in all cases. Indeed, irrespective of the problem and of the used inner linear solver, the algorithm always converges and global errors never exceed 10 −3 . The choice of the linear solver is nonetheless important, as we can see looking at j end and at the total time t tot . Indeed, we see that the AM method requires many more iterations than the BiCGstab(l). Time t tot is higher as well: for instance, in Problem 1 when we use the BiCGstab(1) as linear solver, the LDM computes the solution in less than one fifth of the time required when the AM is used. The difference is yet more remarkable for Problem 3. This had to be expected: indeed, setting c 1 = c 2 = 1,ṽ remains quite small. Thus, the Jacobian is weakly non-symmetric and the AM method is thus slower than the other analyzed linear solvers.
In the BiCGstab(l) method, we also notice that large values of l lead to a reduction of the number of linear iterations but not to a reduction of the computational time, which rather tends to increase for l > 2. This is due to the increased computational cost of each iteration of the BiCG-stab(l) method. Since the BiCGstab(1) (which is equivalent to the Bi-CGSTAB in [22]) can efficiently solve all the analyzed problems, increasing l is thus not needed. In the following, we thus focus on the AM and on the BiCG-stab(1) methods.
Let us then see what happens as c 1 and c 2 are increased, leading to a more nonsymmetric Jacobian matrix. We restrict our analysis to choices of c 1 and c 2 which satisfy the condition h < 2σ min /ṽ max within the discretization used in the above experiments. This condition, moreover, ensures that the AM method converges (see [16]) since it implies that A(u) is an M-matrix.
For instance, let us consider Problem 1 (which is more suitable for this analysis, since σ min is larger) andṽ =ṽ 4 . Taking into account that bothṽ 1 andṽ 2 attain larger values when t is small, it is easy to verify that h < 2σ min /ṽ max is satisfied if c 1 < 200 and c 2 < 96 at t = 0. We thus consider these two values as thresholds for the choice  of c 1 and c 2 . In Table 4, we report the results computed for different choices of c 1 and c 2 . We notice that, as expected, the AM method gets faster and faster as c 1 and c 2 increase, which is when the Jacobian gets more non-symmetric. Even more, beingṽ variable and dependent on t, the non-symmetry of the Jacobian vary at different time levels. As mentioned above, withṽ =ṽ 4 the Jacobian gets less non-symmetric at t increases, since t damps the values ofṽ. So, we expect the AM method to become less efficient as t increases. Eventually, when t is large, the BiCGstab(1) should become competitive also when c 1 and c 2 are large. In order to confirm this, let us then see what happens for different choices of c 1 and c 2 . The results are reported in Tables 5,  6 and 7 and represented in Figs. 1, 2, and 3.
In Table 5 and in Fig. 1, we consider a weakly non-symmetric case, where we set c 1 = c 2 = 1. In this case, the Jacobian is always almost symmetric and the BiCG-stab is evidently faster at any time level.
In Table 6 and in Fig. 2, we instead set c 1 = 50 and c 2 = 95. Non-symmetry is thus much more relevant, and we see that at the beginning the algorithm using the AM as linear solver is twice as fast as the one using the BiCG-stab(1). However, the difference gets smaller as t increases, until, for t > 0.8, the BiCG-stab(1) prevails. The total time required for computing the solution from t = 0 to t = 1 is however smaller when using the AM method, as already seen in Table 4. Choosing t = 1 as final time, the AM thus prevails.  This is even more true for the case c 1 = 195 and c 2 = 95, whose results are reported in Table 7 and in Fig. 3. Here, the AM method tends to be faster also at terminal time levels.

Effect of discretization, α and initialization of linear solver and tolerances
We complete the analysis of the algorithm by summarizing some results that show what happens varying the discretization, changing the initialization and/or stopping criteria of the linear solvers and when α is different from zero.
In the following, when not otherwise specified, we consider Problem 1 withṽ = v 4 , c 1 = c 2 = 50, α = 0 and we again set N = 250 and t = 0.1. With this choice, the Jacobian presents a significant non-symmetry, which, however, is not so big as to make the AM method more efficient than the BiCGstab(l), as we can see in Table 4. The following results have thus been obtained using the BiCGstab(1) as inner linear solver.
Let us start by considering the liner solver itself and see what happens if we modify initializations and stopping criteria. The results are reported in Table 8. In the second column of the table, we declare which initialization of the stopping criterion of the linear solver is used. No bound onˆ (k+1) denotes that we are using (32) to define the tolerance of the linear solver at the first Newton iteration, whileˆ (k+1) ≥σ ν+1 indicates that we are using the tolerance in (34). Finally,ˆ (k+1) ≥σ¯ denotes an intermediate case, where the bound on the tolerance of the linear solver at the first Newton iteration is given in a non-dynamical way by using¯ .
Global errors do not change in all analyzed cases: indeed, some differences could be spotted only if we considered two more decimal digits. However, different choices of initializations and stopping criteria deeply influence the efficiency of the algorithm, here represented by the total number on linear iterations required to compute  the solution at the last time level, j end . Indeed, since we use always the same linear solver, j end is a better indicator than, e.g., t tot , which is less reproducible. We notice that j end is ten times larger when the starting vector of the linear solver is initialized by zero instead of by the solution at the previous Newton iteration, u (ν+1,k) . This is consistent with what happens for constantṽ and also with what we expected: indeed, the choice u (ν+1,k+1,0) = u (ν+1,k) was made in order to choose a starting vector closer to the solution.
Analogously, we notice that we achieve a consistent reduction of the computational cost with no loss in accuracy if we apply the initialization of the tolerance of the linear solver described in Section 5.2. Indeed, j end is reduced by about 20% when u (ν+1,k+1,0) = 0 and by almost 50% when u (ν+1,k+1,0) = u (ν+1,k) , with no increase of global errors.
Passing to the analysis of the discretization, let us set u (ν+1,k+1,0) = u (ν+1,k) and (k+1) ≥ σ ν+1 and let us change the number N of points in which the space domain is discretized. In this regard, we also note that, since c 1 = c 2 = 50, the condition h < 2σ min /ṽ max holds only for N > 104. We get the results in Table 9, where EOC denotes the experimental order of convergence. We notice that global errors decrease as N increases. On the other hand, also the initial residual res 0 increases, leading to higher number of lagged iterations and, thus, of linear iterations, with j end that roughly doubles when N is doubled. The increase in computational times is, then, steeper, since increasing N means increasing also the order of all the matrices involved, thus making the entire solution process more complicated.
Passing to time discretization, Table 10 reports the results obtained as t is varied. Also in this case, global errors decrease as t is reduced. Moreover also res 0 decreases with t, implying that also the times required for computing the solution at each time level get smaller (e.g., see t last , which refers to the time needed for computing the solution at the last time level). However, a smaller t implies that we also need more time levels in order to compute the solution from t 0 to t f . Thus, the total time t tot increases as t is reduced. However, we cannot arbitrarily increase the size of the time step, since the conditions in Theorems 1 and 2 are eventually not satisfied for excessively large time steps (if h is not increased as well).  For completeness, we finally see what happens when α is not zero. The results are reported in Table 11.
We see that α does not have a significant effect on the results, as we expected: indeed, it only acts on the diagonal elements, strengthening (due to its sign) the diagonal dominance of A(u). The first three rows of Table 11 show that, however, a larger α can have the effect of reducing the complexity of the problem (since j end decreases as α increases) and, possibly, also the global errors. This happens also when α is variable: e.g. α = 10/(10 −3 + x + y) 2 leads to α ∈ [10 4 , 4.998] in the domain we are considering and j end is smaller than in the cases where α is smaller. However, when α is variable it is harder to evaluate the effect of α itself, since its value can vary significantly in different parts of the domain.

Discretization of the convective term and convection-dominated problems
Up to now, we have considered the solution of systems arising from space discretizations where the convective term is discretized by central finite differences. This choice is convenient in the analysis of the algorithm and it also makes so that the discretization of the convective term has the same order of accuracy as the discretization of the diffusion term (i.e., O(h 2 )). However, as remarked in Section 2.1, it also has the shortcoming that A(u) is an M-matrix only if h satisfies (4). Evidently, satisfying this condition requires smaller and smaller values of h as the problem gets more convection-dominated, with important consequences on the convergence theorems. Moreover, in convection-dominated problems, it is especially desirable that A(u) be an M-matrix so to be able to use the AM method, which, as seen in the previous subsections, is particularly well-suited to solve systems with non-symmetric matrices.  Considering other discretizations of the convective term, it is nonetheless easy to notice that A(u) is an M-matrix for any choice of h if we instead employ the upwind discretization, at the cost, however, of reducing the order of accuracy to O(h). Thus, we here consider this alternative discretization in order to show numerically its behavior and to provide an useful framework for convection-dominated problems. Further information on this case (referred to the case of constant velocity term) can be also found in [10].
Let us then consider, for instance, Problem 1 withṽ =ṽ 4 . In Table 12, we report the results obtained with different discretizations as the value of c 1 and c 2 increases. Observe that, in all analyzed cases, condition (4) is not respected, as c 1 and c 2 exceed the threshold values reported in Section 6.1 (which are 200 and 96, respectively).
As we expected, the algorithms employing central finite differences on the convective term incur in difficulties as convection becomes more dominant. The AM method, in particular, immediately fails. Algorithms using BiCGstab(l) method work better (especially when l is quite large), but they ultimately fail as convection becomes more dominant. Moreover, also when they work, they are less efficient than the corresponding algorithm employing the upwind discretization. On the other hand, if the discretization is performed by an upwind scheme, the AM method always converges without presenting any problem (on the contrary, the algorithm gets faster with larger values ofṽ). This was to be expected, since the Jacobian is now an M-matrix, as required by the convergence of the AM method.

Conclusions
We have presented an iterative procedure which solves nonlinear steady and nonsteady reaction-convection-diffusion equations with a variable velocity term. We have discretized the partial differential problems and proved that the LDM applied to the resulting nonlinear algebraic systems converges when some assumptions are satisfied. In this context, we have also studied the uniform monotonicity of the finitedifference operator, which is crucial for the uniqueness of the solution and for the convergence of the LDM itself.
We have then described the solution of the weakly nonlinear algebraic systems generated by the LDM and provided some details over the implementation of the entire procedure.
Finally, we have provided several numerical experiments showing the behavior of the LDM in a variety of situations. We have showed that the LDM can successfully compute the solution of all the considered test problems. We have also pointed out that the choice of the linear inner solver affects the efficiency of the method: indeed, the fastest linear solver depends on the non-symmetry of the coefficient matrix of the linear system, which, in turn, is determined by the velocity termṽ. We have then showed that efficiency is affected also by other factors, including the refinement of the discretization and the choice of starting vectors and stopping criteria.