# On the Lagged Diffusivity Method for the Solution of Nonlinear Finite Difference Systems
## 1 Introduction  
Let us consider the nonlinear, non-steady reaction-convection-diffusion equation  
\[ \begin{array}{ll}\frac{\partial u}{\partial t}=\
abla\times(\sigma\
abla u)- \tilde{\
u}\times\
abla u-\alpha u-g+s&(x,y)\in\Omega,\,t>t_{0}&\text{Diff. Eq.}\\ u(x,y,t_{0})=u_{0}(x,y)&(x,y)\in\overline{\Omega},\,t=t_{0}\geq 0&\text{IC}\\ u(x,y,t)=u_{1}(x,y,t)&(x,y)\in\partial\Omega,\,t>t_{0}&\text{Dirichlet BCs}, \end{array}  \tag{1}\]  
where \(t_{0}\) is the initial time, \(u=u(x,y,t)\) is the density function, \(\sigma=\sigma(x,y,u)>0\) is the diffusion coefficient, \(\tilde{\
u}\) is the velocity vector (assumed constant), \(\alpha=\alpha(x,y)\geq 0\) is the absorption term, \(-g=-g(x,y,u)\) is the rate of change of the reaction and \(s=s(x,y)\) is the source term. This equation is important in various different fields and it can describe several processes and applications. 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 [PERSON]'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 [PERSON]'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 \(\theta\)-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 \(\bar{v}\) 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 \(\sigma\), \(\alpha\), \(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] (SS2.3) or [7]).