On the Lagged Diffusivity Method for the Solution of Nonlinear Finite Difference Systems

Abstract: In this paper, we extend the analysis of the Lagged Diffusivity Method for nonlinear, non-steady reaction-convection-diffusion equations. In particular, we describe how the method can be used to solve the systems arising from different discretization schemes, recalling some results on the convergence of the method itself. Moreover, we also analyze the behavior of the method in case of problems presenting boundary layers or blow-up solutions.

If we discretize Equation (1), at each time level the differential problem can be approximated by a nonlinear, algebraic system.However, computing the solution of the discretized system is generally complicated.For instance, if we use a Newton's method, we incur in the problem of computing the Jacobian matrix, which is usually difficult and computationally onerous.Therefore, it is suitable to first linearize the discretized system by setting up an iterative procedure which "lags" the diffusivity.In this way, we obtain the so-called Lagged Diffusivity Method (LDM).The analysis of this kind of strategies can be traced back to [1].More recently, the LDM applied to steady and non-steady diffusion equations with reaction and/or convection terms has been studied in [2] and following works.
The nonlinear algebraic system is thus replaced by a series of systems of weakly nonlinear algebraic equations, whose Jacobian matrix is much easier to compute.Therefore, these systems can be easily solved by a Newton's method.Lastly, at each Newton iteration we solve a system of linear algebraic equations (which we again do iteratively, for enhanced efficiency).This paper addresses some topics which have not been considered in earlier works.In particular, we • compare different discretizations.In particular, we consider what happens when the space is discretized by an upwind method or by central finite differences.We also describe different time discretizations, considering the θ-method and the IMEX schemes; study numerically what happens in some particularly critical cases, where some smoothness assumptions are not satisfied.This is the case of boundary layer and blow-up solutions.
Regarding the discretization, this is relevant since the discretization can modify the properties of the discretized systems.This can be exploited to improve the solution procedure: for instance, we notice that the upwind method can significantly reduce the computational times when ṽ is large (at the cost of a worse order of approximation).
The analysis of blow-up solutions is instead interesting for two reasons.First, we evaluate the behavior of the method in critical situations where the LDM is not analytically proven to converge.Indeed, the convergence of the method requires some smoothness assumptions on σ, α, g and s which are not satisfied in proximity of blow-ups and boundary layers.Second, we analyze numerically the problem, while several papers regarding blow-up solutions deal with the analytical problem of determining if and when the blow-up occurs (e.g., if T is finite or not) and they are often referred to 1D differential problems with Neumann boundary conditions.Regarding analyses involving multidimensional domains and/or Dirichlet boundary conditions, it is instead worth citing [3][4][5], which are concerned with problems similar to the ones we analyze in the following.In particular, criteria for determining whether the blow-up occurs in a finite time or not and for estimating the blow-up time itself are given by using the Hopf's maximum principle (e.g., see [6] ( §2.3) or [7]).Other interesting papers on blow-up include [8,9], while [10] is a comprehensive review on blow-up in parabolic equations.
In Section 2 we analyze the discretization.In particular, we describe how the discretized system changes when space is discretized by central finite differences or by an up-wind scheme and when time is integrated by a θ−method or by IMEX strategies.
In Section 3 we then summarize the method and consider how it is affected by the discretization itself.We also report some known results on the monotonicity of the finite-difference (FD) operators and on the convergence of the method.
In Section 4 we provide some details on the implementation, an algorithm and a short summary of the solution method.
Lastly, Section 5 is devoted to the numerical experiments.First, we show numerically the differences which arise when the discretization is changed.Then, we study problems characterized by blow-up and by boundary layer solutions in order to analyze the behavior of the LDM in critical cases, when the smoothness assumptions are not satisfied.

Space Discretization
For simplicity, let us assume Ω to be a 2D bounded rectangular domain.Let us superimpose a grid Ω h = Ω h ∪ ∂Ω h to Ω.For simplicity, we use uniform mesh size h along x and y.Thus, Ω h is made of the mesh points (x i , y j ), i = 1, ..., N, j = 1, ..., M with: Using this grid, u i,j (t) : Ω h → R denotes the grid function which 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.This grid function satisfies the initial condition for (x i , y j ) ∈ Ω h and t = t 0 and the Dirichlet boundary conditions for (x i , y j ) ∈ ∂Ω h , t > t 0 .
We can now approximate the derivatives of Equation ( 1) by finite difference methods.In this regard, several choices are possible.In particular, it is interesting to analyze the space discretizations obtained by central finite differences and by the upwind scheme.Indeed, this choice influences the properties of the matrix of the discretized system, as we will see in the following.Regarding the upwind scheme, we consider, for instance, the case of positive ṽ.If ṽ < 0 the backward difference quotients which approximate the first-order derivatives are simply replaced by forward difference quotients.
Thus, denoting the forward finite differences (in x and in y respectively) by ∆ x and ∆ y , the backward finite differences by ∇ x , ∇ y and the central finite differences by δ x , δ y , we have:

•
Discretization by central finite differences (error O(h 2 )) • Discretization by upwind scheme (error O(h)), case ṽ > 0 In both cases, we can write the discretization in a more convenient way by defining Central FD Upwind Then, setting Li,j = L i,j + Li,j (analogously for Bi,j , Ri,j , Ti,j and Di,j ), the right hand-side of Equations ( 2) and (3) can be written formally in the same way as Bi,j u i,j−1 (t) + Li,j u i−1,j (t) − Di,j u i,j (t) + Ri,j u i+1,j (t) + Ti,j u i,j+1 (t) − g(x i , y j , u i,j (t)) + s(x i , y j , t). (5)   In the following, we denote σ(x i , y j , u i,j ) by σ(u i,j ) for compactness of notation.Finally, we use matrix notation to write the entire system of ordinary differential equations arising from the space discretization.In this regard, we use the row lexicographic ordering for mesh points and define P k = (x i , y j ) with k = (j − 1)N + i, for i = 1, ..., N, j = 1, ..., M. We can then write the vector u(t) containing all the u i,j (t) at internal mesh points.Thus, the space discretization leads to Again, we notice that formally this expression is valid both when the discretization is performed by central FD and when the upwind scheme is used.The following differences are however present: • when we use the upwind scheme, the order of the approximation is O(h); when we use central FD, the order of the approximation is O(h 2 ); • when we use the upwind scheme, A(u(t)) is irreducibly diagonally dominant [11] (p. 23) irrespective of the choice of h.Having also positive diagonal elements and non positive off-diagonal elements, it is always a non-singular M-matrix [11] (p. 91).With central FD, A(u(t)) is irreducibly diagonally dominant (and a non-singular M-matrix) only if h satisfies a condition on the step-size.In particular, by Equation (4) if is easy to find that this condition corresponds to The vector b(u(t)) in Equation ( 6) derives from the Dirichlet boundary conditions u(x i , y j , t) = u 1 (x i , y j , t) at points (x i , y j ) ∈ ∂Ω h .Considering the Dirichlet BCs in Equation (1) and using, again, the row-lexicographic order, its components are where the coefficients L, B, R, T depend on u at the inner point which is neighbor to the considered point of ∂Ω h .At corners, where the same inner point is neighbor to two boundary points, we then have two contributions: e.g., b , y j , u k ), with i = 1, ..., N, j = 1, ..., M and k = (j − 1)N + i.It is, then, easy to notice that G(u) is a diagonal mapping [12] (p.11).

Time Discretization
We now consider time discretization.Let us introduce a time spacing ∆t, which defines a series of time levels: the n-th time level t n is defined by t n = n∆t + t 0 , n = 0, 1, ....The time derivatives of Equation ( 6) can be approximated by several different schemes.We could use, for example, the well-known θ−method (e.g., see [13]), which is completely implicit when θ = 0.An alternative, as noted in the Introduction, is given by IMEX methods, which have been employed to solve diffusion equations when convection or reaction terms can be treated explicitly.
For instance, let us apply the θ-method.Denoting by u n the approximation of u(t n ) (with u(t n ) solution of Equation ( 6) at t = t n ) and by s n ≡ s(t n ), we get: for n = 0, 1, ... and with 0 ≤ θ ≤ 1.By simple algebraic manipulations, we find that at each time level n = 0, 1, ... the vector u n+1 is given by the solution of the nonlinear algebraic system where τ = ∆tθ, I is the µ × µ identity matrix and w = w n ∈ R µ is a vector containing the known terms, defined as If we use an IMEX method, we similarly get a nonlinear algebraic system to solve at each time level.The difference is that some terms are treated explicitly.Thus, if, for instance, the reaction term is treated explicitly, at the time level n + 1 the mapping G(u) is evaluated at time n; G(u) then becomes a known term which appears solely in the known vector w.
Analogously, if the velocity is not linear, it could be useful to treat explicitly the convection term.In this case, the part of A(u) depending on ṽ would appear solely in the known vector w.

The Lagged Diffusivity Method
We now introduce the lagged diffusivity method to solve the discretized system, referring, for instance, to Equation (9).
Solving Equation ( 9) by a Newton's method would be onerous, since the Jacobian matrix of F(u) is generally hard to compute.The lagged diffusivity method, on the other hand, linearizes the discretized system by iteratively "lagging" the diffusivity term.At each iteration we must thus solve a weakly nonlinear algebraic system, whose Jacobian matrix is significantly easier to compute.
Let us study more in detail the procedure applied to Equation (9).Lagging the diffusivity means that the new iterate u (ν+1) of the (ν + 1)-th lagging iteration is the solution of the weakly nonlinear system When A(u) is a nonsingular M-matrix, I + τA(u) is evidently nonsingular ∀u ∈ R µ .Thus, we now have to solve Equation ( 10) at each lagging iteration.We can do this by a Newton's method (or, in general, by an iterative method) which is stopped when the residual satisfies a stopping criterion where • denotes the Euclidean norm and ν is a given tolerance such that ν → 0 for ν → ∞.When the inner solver converges, we find u (ν+1) , which is used, in turn, to lag the diffusivity at the (ν + 2)-th lagged iteration.Proceeding in the same way, we thus find u (ν+2) and so on for the following iterations.

Effect of the Discretization
In the previous paragraphs, we applied the lagged diffusivity method to Equation ( 9), which was obtained by discretizing time using a θ−method.The time discretization, however, affects some properties of the lagging iteration.For instance, if we can use an IMEX method and treat G(u) explicitly, the lagged diffusivity method completely linearizes the discretized system.
Analogously, an IMEX scheme can simplify the analysis of cases with non-constant velocity terms: if we can treat non-constant terms ṽ(u) explicitly, all the convection terms are evaluated at the previous time level and are, thus, known.Both monotonicity and convergence analyses can then be traced back to the results obtained for constant ṽ.
Figure 1 summarizes all these cases, representing which functions are unknown at each time level and at each lagging iteration.
Finally, although formally the system does not change, also the space discretization has an effect on the solution procedure.Indeed, we noticed earlier that A(u) is always an M-matrix when the upwind discretization is used, while we have a condition on the step-size h when central finite differences are used.If A(u) is an M-matrix, also the coefficient matrices of the linear systems arising from applying a Newton's method to solve Equation ( 10) are M-matrices.This is important since some linear solvers converge only when the coefficient matrix is an M-matrix.We can thus always use these solvers when we use the upwind discretization, while we can do so only when condition Equation ( 7) is satisfied when we use central FD.This is of great importance especially when ṽ is large: indeed, in this case the coefficient matrix of the linear systems is highly non-symmetric and we would like to use splitting methods such as the Arithmetic Mean (AM) [14], which are particularly efficient in this case but require that the coefficient matrix is an M-matrix.Thus, also considering that Equation (7) implies that h must be really small when ṽ is large, we may want to use an upwind discretization in order to use these linear solvers.This is further studied in the numerical experiments.Here A 1 (u) represents the part of A dependent on σ(u) and Ã(u) represents the part of A dependent on ṽ(u).

Monotonicity of the Finite-Difference Operator and Convergence Analysis
Let us assume that the following smoothness properties hold true: 1.
the functions σ and g are continuous in u and continuously differentiable in (x, y).Moreover, the functions α and s are continuous in their variables; 2.
for fixed (x, y) ∈ Ω, the function g is uniformly monotone in u (uniformly in x, y) with constant c > 0 [15] (p.114); moreover, it is continuously differentiable in u.
Under these smoothness assumptions, the monotonicity of the finite different operator and the convergence of the algorithm are assured for grid functions u i,j belonging to a set B ρ,β , where ρ is a bound on the discrete and β is a bound over the absolute value of the backward difference quotients of u i,j .
Thus, considering for instance the discretization by central finite differences and time integration by the θ−method, the following theorems hold true Theorem 1.Let u * ∈ B ρ,β be a solution of the nonlinear system F(u) = 0.If (13) then F(u) is uniformly monotone in B ρ,β and the solution u * of Equation ( 9) is unique.
Theorem 2. Let u * ∈ B ρ,β be the solution of the nonlinear system F(u) = 0 defined in Equation (9) with A(u) non-singular M-matrix and G(u) diagonal mapping.We assume the smoothness Properties 1-4 to be satisfied and that condition Equation (13) on the monotonicity of F(u) is satisfied.
The proofs of these theorems and more information on B ρ,β can be found in [16].In the literature it is also possible to find proofs of the convergence of the procedure for less general cases (e.g., see [2] for steady state reaction diffusion problems).Proceeding analogously, similar results can be easily found also for the other discretizations.

Solution Procedure
In this section, we present some insights on the implementation that is used in the numerical experiments.This implementation follows the description in [2], with the difference that the inner systems are here solved by an inexact, simplified Newton iteration and that we use the correction on the initialization of tolerances presented in [16] and here summarized in Section 4.2.Then, we also briefly summarize the systems arising from the discretization and from the LDM and provide an algorithm which describes the entire procedure.

Solution of the Inner Systems
We solve the lagged system of Equation ( 10) by a simplified Newton method [12] (p.182).At each Newton's iteration we then have to solve a system of linear equations.For achieving the maximum efficiency, we do this approximately by an iterative linear solver, whose tolerance is chosen so that the solution procedure of Equation (10) makes up an inexact (simplified) Newton method [17].
In particular, denoting the Newton's iteration by a second superscript k = 0, 1, ..., the simplified Newton iteration at the (ν + 1)-th lagged iteration consists in computing the solution ∆u (k+1) = u (ν+1,k+1) − u (ν+1,k) of the linear system where F ν (u (ν+1,0) ) is the Jacobian matrix of F ν (u) evaluated at u (ν+1,0) , which is where G (u (ν+1,0) ) is the Jacobian matrix of G(u) evaluated at u (ν+1,0) .From Equation (15), we can appreciate the advantage of the lagging iteration.Indeed, since the LDM linearizes A(u), it is easy to compute the Jacobian matrix of F ν (u), since the only nonlinear term is the diagonal mapping G(u) .Finally, the linear system arising at each Newton iteration is solved by an iterative solver.We henceforth refer to the inner linear iteration by a third superscript j.

Starting Vectors and Stopping Criteria
Starting vectors and tolerances of the iterative procedures are chosen as in the diagram in Figure 2, where η ∈ (0, 1) and ν is initialized by 1 = ˜ 0 F(u (0) ) , with ˜ positive constant smaller than 1.

Lagged diffusivity method
Simplified Newton method Iterative linear solver The choice of the starting vectors is made so to start as close as possible to the solution of the system.Moreover, it also allows to write conveniently the linear system arising at each Newton iteration (e.g., see [16]).The stopping criteria, on the other hand, descend directly by the description of the LDM and by the choice of an inexact Newton iteration to solve the inner systems.

Starting vector
The starting vectors determine also the initialization of the tolerances: for example, at the first Newton iteration of the (ν + 1)-th lagged iteration, the tolerance of the linear solver is ) .We can also easily notice that all the tolerances are initialized at η F 0 (u (0) ) = η F(u (0) ) = 1 at the first lagged iteration.
However, if we do not apply any correction, the initial tolerance of the linear solver, ˆ 1 , may become smaller than the minimum possible tolerance of the simplified Newton method, η ν+1 .Thus, ˆ 1 would be smaller than what could be required by an inexact Newton method, leading to an increase in computational cost with no foreseeable improvement in accuracy.
To avoid this situation, we simply impose η ν+1 as the minimum acceptable tolerance of the linear solver: This does not lead to a worsening of the accuracy of the lagged procedure (indeed, η ν+1 is the tolerance of the simplified Newton method), while significantly reducing the computational cost.

Summary of the Solution Method
The diagram in Figure 3 summarizes the equations and the systems that we used, illustrating how the partial differential equation is first transformed in a system of ordinary differential equations and then in a series of nonlinear algebraic systems, which are afterwards progressively linearized by the iterative methods.Finally, we can write the algorithm summarizing the entire procedure as follows.

9:
Compute residual r j+1 of the linear system Compute Newton residual F ν (u (ν+1,k+1) ) 14: if F ν (u (ν+1,k+1) ) ≤ ν+1 then break 15: Update linear solver tol.: ˆ k+1 = η F ν (u (ν+1,k+1) ) Update vectors and matrices: find F ν+1 (u (ν+1) ) 19: end for 23: n = n + 1 24: end for Here, break means that we exit the current loop and go back to the outer one.On the other hand, return means that we exit from the entire procedure and the algorithm returns the final result.If the smoothness assumptions in Section 3.2 are satisfied (and if the inner linear solver converges), then the convergence of Algorithm 1 is ensured by Theorem 2 or by analogous theorems, depending on the chosen discretization.

Numerical Experiments
We consider a square domain Ω = [0, 1] × [0, 1], which we discretize by a uniform grid of N × N inner points.We build test problems from known solutions u * (specified in the following) by computing an appropriate source term.When not otherwise specified, time discretization is performed by a θ−method with ∆t = 0.1 and θ = 0.5 and we compute the solution from a starting time t 0 = 0 to a final time t f = 1.
In the test problems we choose α = 0, σ = 0.4 + 0.5u and g = 100e 0.5u , which has the form of the radiation model.We instead consider several different choices of u * and of ṽ.
The linear systems arising at each Newton iteration are solved considering Krylov and splitting methods.In particular, we use the Preconditioned Conjugate Gradient (PCG) method [18] when the coefficient matrix of the linear system is symmetric and positive definite (i.e.ṽ = 0).Here, the used preconditioner is the diagonal matrix M whose elements m i,i are given by the l 2 -norm of the i-th row of the coefficient matrix of the linear system.We instead compare the AM method [14] and the BiConjugate Gradient-stabilized method with parameter l (BiCGstab(l)) [19] when ṽ = 0.
We choose = 10 −4 ; the other stopping criteria are determined as in Section 4. The maximum number of simplified inexact Newton iterations is k max = 500 and the maximum number of linear iterations is j max = 10, 000.In the following subsections, ν * , k * and j * denote, respectively, the total number of lagged, Newton and linear solver iterations required for computing the solution at the last time level.By res 0 and res we instead denote, respectively, the initial residual and the last residual of the linear system at the last time level, while err h and rel.err 2 denote the global error (in l 2 (Ω h ) and in relative Euclidean norm respectively).

Effect of the Discretization Scheme
Let us analyze the effect of the space discretization.In Table 1 we report the results obtained for different choices of ṽ with u * = u * 1 := (1 + x − y) 3 t.We set N = 250; with this choice, Equation ( 7) is satisfied for ṽ1 , ṽ2 < 200.In Table 1, we see that the two discretization schemes are identical for ṽ = 0, also in terms of global errors and number of iterations.Indeed, in this case, the first-order derivative in the differential equation is multiplied by a null term.Thus, A(u) is symmetric positive definite and the two discretizations are equivalent.
In the other cases, both methods always work properly when ṽ respects the spacing condition, although global errors are larger when we use the upwind discretization.As mentioned above, this is to be attributed to the discretization of the first-order derivative by forward or backward finite difference quotients, leading to a larger discretization error.
Lastly, when the spacing condition is not respected anymore, we can have problems in the solution of the linear system.Indeed, the convergence of the AM method requires that the coefficient matrix is an M-matrix.Although the AM still works for values of ṽ for which this is not satisfied (e.g., ṽ = (300, 300) T ), we soon have problems for larger values of the velocity term.For example, ṽ = (500, 500) T (and any other larger velocity term) leads the AM to stall.Also the BiCGstab(1) presents some problems and it already stalls for ṽ = (300, 300) T .Contrary to the AM, however, the BiCGstab(1) does not explicitly require that the coefficient matrix is an M-matrix.Its stalling is, however, not surprising.Indeed, the BiCGstab(1) (which is equivalent to the BiCG-STAB [20]) is known for stalling especially when the coefficient matrix is real but its eigenvalues are complex with large imaginary part.This could certainly be the case, since, irrespective of the chosen discretization, the coefficient matrices are non-symmetric when ṽ = 0. Complex eigenvalues are, thus, possible.
Then, it is interesting to see what happens increasing l in the BiCGstab(l) method.Indeed, l has been introduced exactly to overcome these difficulties.Let us then consider yet larger values of ṽ and analyze what happens using BiCGstab(2) and BiCGstab(4) as well.The results are reported in Table 2.As we expected, the AM and the BiCGstab(1) are not able to solve the linear systems arising when ṽ is large and central FD are used.Indeed, they stall and the maximum number of allowed linear iterations is reached.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 coefficient matrix is now an M-matrix, as required by the convergence of the AM method.The BiCGstab(1), instead, keeps stalling for large values of ṽ, as we could expect from what we said in the previous paragraph.However, we can see that this can be overcome by increasing the parameter l of the BiCGstab(l): the algorithms using BiCGstab (2) or BiCGstab(4), indeed, converge in all cases.We can also see the big difference between the global errors obtained when discretizing by central finite differences or by an upwind scheme.Indeed, they are in the order of 10 −6 for central finite differences and in the order of 10 −3 for the upwind scheme.However, it is also to be noted that the upwind scheme is much faster, especially when we use the AM as inner solver.Thus, we may opt for central finite differences with, for instance, BiCGstab(2) as linear solver if we require smaller errors, while we may choose an upwind method with AM if we require lower computational cost.

Blow-Up
We now apply the LDM to problems with blow-up and boundary layer solutions.In particular, we consider the solutions in Figure 4, where T denotes the blow-up time.
Blow-up zones Plot of u * 4 at t T u * 5 = 1 sin((T − t)π) + sin(πx) + sin(πy) , Blow-up zones Plot of u * 5 at t T We set T = 1 for u * = u * 2 and u * = u * 3 ; instead, u * 4 and u * 5 require T < 1 by the periodicity of the trigonometric functions.In these two cases we therefore set T = 0.8.
With no loss of generality, the velocity vector ṽ is here set equal to zero and the system matrix A(u) is a Stieltjes matrix.We thus use the PCG solver for solving the linear systems arising at each Newton iteration and discretize space by central finite differences.
In the following results, max u denotes the maximum of the computed solution vector at the last computed time level, while max u * denotes the maximum of the analytical solution u * at the last computed time level in the point where max u occurs.Analogously, we define min u and min u * .

Preliminary Results and Analysis
Let us first analyze how the computed solution evolves as the blow-up time approaches.To this end, let us start from t 0 = 0 and see what happens at different times, until reaching T; at first, we here choose ∆t = 0.1 and N = 100.In Figure 5 we report what we get considering e.g., u * 2 (example of solutions with blow-up on one corner), u * 4 (example of solutions with blow-up on one or more edges) and u * 5 (example of solutions with blow-up on more corners).In Table 3 we report some information on errors and data of the iterative procedure.Of course, in all these cases, the iterative procedure stops when we are still quite far from the blow-up time: indeed, using ∆t = 0.1, the last computable time step is at t = 0.9 when T = 1 and t = 0.7 when T = 0.8.Indeed, the blow-up occurs at t = 0.9 + ∆t and at t = 0.7 + ∆t respectively.We therefore need a smaller time step in order to get closer to the singularity.However, we notice that the computed solutions behave as the analytical ones (as proven by the errors in Table 3) and that, as the blow-up time approaches, they get more and more similar to the ones in Figure 4.
It is also interesting to notice that the residual computed at the beginning of the lagged iteration, res 0 , gets larger as we get closer to T. This can be viewed as an indicator of incoming blow-up: indeed, by the initialization u (0) = u n , at the time step n + 1 we have res 0 = F(u n ) .Therefore, res 0 is the residual of F(u) = 0 as in (9) for u n at time n + 1.So, if the solution is smooth in the time variable (and ∆t is not excessively large), we do not expect res 0 to be large: in these hypotheses, the solution u n+1 will not be dramatically different from the solution u n .However, if we have a singularity at some time T, the solution vector will change faster and faster as we approach T, leading to progressively larger initial residuals res 0 .In this light, the analysis of res 0 gives information similar to the analysis of the derivative of u * with respect to t.

Effect of Refining the Time Grid
Taking into account the previous results, let us set ∆t = 10 −3 .We now start our analysis from t 0 = T − 0.1, so that we do not need an extremely large number of iterations in order to get close to the blow-up zone.
In Figure 6 we show how the computed solution looks like at the last computable time step before the blow-up; we call t f this final time.In Table 4 we report some data on the lagged iterative procedure at t = t f .The first thing we notice is the extremely large values of res 0 , which are even more remarkable considering that ∆t has been reduced.This can be explained as in the previous paragraph by considering that the solution changes more and more rapidly as we approach T. A gradual increase of res 0 with t is indeed what can be observed by considering the value of res 0 at different time levels.
The analysis of the plots reveals that the lagged diffusivity procedure is able to reproduce the behavior of all the solutions near their respective blow-up zones.On the other hand, the global errors gradually increase, but they generally remain in the order of 10 −3 ÷ 10 −2 .Indeed, only for u * 4 they exceed 10 −1 near T.
Notwithstanding this, maximum and minimum values of the solution vector u correspond to the maximum and minimum values of u * (computed in the inner grid points) to at least the second decimal digit in all the cases.It is to be noticed that the maximum of u is the point where the solution is nearest to the blow-up one; a good precision in evaluating the solution in this point is thus important for accurately reproducing the blow-up behavior.
Finally, we also notice that these maxima and minima of u are not so large after all: we are near to the blow-up time, but the maximum of the solution never exceeds values ∼ 30, nonetheless.This cannot be explained by numerical errors in the solution procedure, since we said that the maximum of the computed u is really similar to the value of u * in the grid point where the maximum occurs.Different causes concur to this: • the blow-up occurs on the boundary, where the solution is given by the Dirichlet boundary conditions.Thus, at inner points, where we compute the solution vector, the maximum of u is always smaller than the maximum of u at boundary points when t ∼ T; • as a consequence of the previous point, the smaller the number of grid points, the smaller max(u).Indeed, we get farther from the boundary and, thus, from the singularity.This is also shown in the next subsection, where we observe that max(u) increases when we increase the number of points of the space grid; • the singularity occurs abruptly when t = T is reached, so the solution is still quite small also at times quite close to T. Next, we show that by further reducing ∆t, we are able to get nearer to T, leading to an increase of max(u).
All these facts can be better viewed by considering a short example.Let us consider u * = u * 2 and N = 100.We remind that the blow-up here occurs at T = 1.00 in the point (1,0).The solution at t = t f = 0.993 at the inner point closest to (1, 0), which is (x n , y 1 ) = (100/101, 1/101) is: The algorithm, then, does not stall due to large values of u, but because its derivatives become large.Indeed, considering the largest first-order derivative in x, we have This makes so that the bound β on backward difference quotients increases as the blow-up approaches.Then, it is harder to satisfy the monotonicity condition of Theorem 1, which, in turn, affects the uniqueness of the solution and the convergence of the LDM.We can try to counterbalance the larger β by reducing ∆t (and, thus, τ, which appears in the condition of Theorem 1), as we verify in the next subsection.However, further decreasing ∆t means increasing the number of time levels and, thus, the computational time.Nonetheless, we can compensate this effect by starting our computations nearer to the blow-up time.

Effect of Refining the Space Discretization
Let us then analyze what happens when the space grid is refined; in this regard, the results obtained with N = 250 and ∆t = 10 −3 are reported in Table 5.
We notice that t f slightly decreases, while the maximum and the minimum of u in t f are about the same as previously computed for N = 100.Indeed, being the grid finer, the point nearest to the blow-up point is closer (in space) to the blow-up point than in the case N = 100.Thus, the solution gets larger earlier.(1,0).The solution at t = t f = 0.981 at the inner point closest to (1, 0), which is (x n , y 1 ) = (250/251, 1/251) is: ∼ 37.207, as computed.Analogous results can be easily verified with regard to the derivatives.We also notice that global errors decrease due to the finer space grid, while res 0 is of the same order of magnitude as above.
Lastly, we verify that further reducing ∆t allows us to better approach the blow-up time, as remarked at the end of the previous subsection.Thus, in Table 6 we report the results obtained by setting ∆t = 10 −5 and by using as t 0 the last computed time step t f of Table 5.In all cases we are able to get closer to the blow-up time, as expected.Moreover, global errors decrease.

Boundary Layer
Lastly, for completeness, let us consider a stationary case with non-constant ṽ which presents a boundary layer solution.
Let us consider the problem given by and with λ < 0 real parameter and boundary condition u = e λx + e λy on ∂Ω.The solution of this test problem is evidently u * = e λx + e λy , ( which presents a boundary layer as the absolute value of λ increases.Since we do not have any nonlinear reaction term, at each iteration of the LDM we just have to solve a linear system, which we do by the BiCGstab(l) algorithm.The results as |λ| increases are reported in Table 7 and the computed solutions are in Figure 7.We notice that the number of iterations of the linear solver decreases as |λ| increases.The global errors, instead, increase, but they never get larger than 10 −3 , confirming the good approximation of the solution.In particular, Figure 7c, correctly reproduces the profile of a boundary layer solution.

Conclusions
We analyzed the effect of the discretization on the LDM for the solution of nonlinear, non-steady convection-reaction-diffusion equations.Some insights on an efficient implementation of the algorithm are provided and are followed by numerical experiments, which highlight the differences arising when the space discretization is performed by central finite differences or by the upwind method.
The paper is completed by a numerical study on the application of the LDM to problems characterized by blow-up and boundary layer solutions.The numerical experiments show that the LDM is able to correctly reproduce the analytical solution also in proximity of the blow-up time or in presence of a boundary layer.

Figure 1 .
Figure 1.Summary of the nonlinear terms present at each time level and at each lagging iteration for different time discretizations.(a) θ−method; (b) IMEX scheme with explicit treatment of G(u); (c) IMEX scheme with non-constant velocity term ṽ(u) treated explicitly.Here A 1 (u) represents the part of A dependent on σ(u) and Ã(u) represents the part of A dependent on ṽ(u).

Figure 2 .
Figure 2. Starting vectors and tolerances of the iterative methods.

Figure 3 .
Figure 3. Systems arising from the discretization of the initial PDE and from the used iterative methods.

Figure 4 .
Figure 4. Blow-up solutions and representation of where and how the blow-up occurs.

Figure 5 .
Figure 5. Evolution of the computed solution as t approaches T in the cases u * = u * 2 , u * = u * 4 and u * = u * 5 .

Table 1 .
Comparison between discretization by central finite differences and by upwind scheme.

Table 2 .
Comparison between discretization by finite differences and by upwind scheme with high ṽ.

Table 4 .
Error and data of the lagged diffusivity procedure for the last computable solution for the problems in Figure4for ∆t = 10 −3 .

Table 5 .
Error and data of the lagged diffusivity procedure for the last computable solution for the problems in Figure4for ∆t = 10 −3 and N = 250.× 10 9 1.95 × 10 −5 1.04 × 10 −2 1.34 × 10 −3 51.To better see this, we can again consider the short example proposed earlier and analyze what happens for N = 250.Let us again consider u * = u * 2 with T = 1.00 in the point