# On the Lagged Diffusivity Method for the Solution of Nonlinear Finite Difference Systems
## 2 Discretization
### Space Discretization
 \tag{4}\]  
Then, setting \(\hat{L}_{i,j}=L_{i,j}+\hat{L}_{i,j}\) (analogously for \(\hat{B}_{i,j}\), \(\hat{R}_{i,j}\), \(\hat{T}_{i,j}\) and \(\hat{D}_{i,j}\)), the right hand-side of Equations (2) and (3) can be written formally in the same way as  
\[\hat{B}_{i,j}u_{i,j-1}(t)+\hat{L}_{i,j}u_{i-1,j}(t)-\hat{D}_{i,j}u_{i,j}(t)+ \hat{R}_{i,j}u_{i+1,j}(t)+\hat{T}_{i,j}u_{i,j+1}(t)-g(x_{i},y_{j},u_{i,j}(t))+ s(x_{i},y_{j},t). \tag{5}\]  
In the following, we denote \(\sigma(x_{i},y_{j},u_{i,j})\) by \(\sigma(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  
\[\frac{\mathrm{d}\mathbf{u}(t)}{\mathrm{d}t}=-A(\mathbf{u}(t))\mathbf{u}(t)+\mathbf{b}(\mathbf{u}(t ))-\mathbf{G}(\mathbf{u}(t))+\mathbf{s}(t). \tag{6}\]  
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(\mathbf{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(\mathbf{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  
\[h<\min\left\{\frac{2\sigma_{\min}}{\left|\bar{\sigma}_{1}\right|},\frac{2\sigma_{ \min}}{\left|\bar{\sigma}_{2}\right|}\right\}. \tag{7}\]  
The vector \(\mathbf{b}(\mathbf{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})\in\partial\Omega_{h}\). Considering the Dirichlet BCs in Equation (1) and using, again, the row-lexicographic order, its components are  
\[b_{k} =L_{1,j}u_{1}(x_{0},y_{j}) \text{with }k=(j-1)N+1,\,j=1,...,M\] \[b_{k} =B_{i,1}u_{1}(x_{i},y_{0}) \text{with }k=i,\,i=1,...,N\] \[b_{k} =R_{N,j}u_{1}(x_{N+1},y_{j}) \text{with }k=jN,\,j=1,...,M\] \[b_{k} =T_{i,M}u_{1}(x_{i},y_{M+1}) \text{with }k=(M-1)N+i,\,i=1,...,N,\]  
where the coefficients \(L\), \(B\), \(R\), \(T\) depend on \(\mathbf{u}\) at the inner point which is neighbor to the considered point of \(\partial\Omega_{h}\). At corners, where the same inner point is neighbor to two boundary points, we then have two contributions: e.g., \(b_{1}=L_{1,1}u_{1}(x_{0},y_{1})+B_{1,1}u_{1}(x_{1},y_{0})\). Defining \(\mu=NM\), \(\mathbf{G}(\mathbf{u})\in\mathbb{R}^{\mu}\) is the nonlinear mapping of components \(G_{k}(\mathbf{u})=G_{k}(u_{k})=g(x_{i},y_{j},u_{k})\), with \(i=1,...,N,j=1,...,M\) and \(k=(j-1)N+i\). It is, then, easy to notice that \(\mathbf{G}(\mathbf{u})\) is a diagonal mapping [12] (p. 11). Finally, \(\mathbf{s}(t)\in\mathbb{R}^{\mu}\) contains the terms arising from the source term: \(s_{k}(t)=s(x_{i},y_{j},t)\) with \(i=1,...,N,j=1,...,M\) and \(k=(j-1)N+i\).