## 2 Materials and Methods
### Problem Statement and Mathematical Background
At the same time, in the course of the numerical solution of (2), this approach also leads to a linear system to solve:  
\[D_{h,\text{FE}}\mathbf{u}_{h,\text{FE}}=\mathbf{f}_{h,\text{FE}} \tag{7}\]  
Here, the unknown vector \(\mathbf{u}_{h,\text{FE}}\in\mathbb{R}^{N}\) consists of coefficients of a finite element basis, while \(\mathbf{f}_{h,\text{FE}}\) represents the scalar product of \(f\) with the basis elements. Usually, a lifting also has to be applied to incorporate the boundary conditions. In any case, again, the matrix \(D_{h}\in\mathbb{R}^{N\times N}\) can be regarded as an approximation of \(\Delta\) on the entire domain \(\Omega\). Accordingly, in a specific case, shown on the right of Figure 1, we can determine the coefficients as in (4) to obtain  
\[a_{0}=-\frac{1}{h_{x}^{2}}\cdot\frac{44}{9},\ \ a_{1}=a_{2}=a_{3}=a_{4}=\frac{1}{h _{x}^{2}}\cdot\frac{25}{18}\quad\text{and}\quad a_{5}=a_{6}=-\frac{1}{h_{x}^{2 }}\cdot\frac{1}{3}. \tag{8}\]  
In the course of computing these coefficients, we have used the ratio \(h_{s}=\frac{5}{3}h_{x}\). Note that this discretization leads to a second-order convergence with respect to the \(L_{2}(\Omega)\)-norm. To compare the above cases, we also give the position of the nonzero elements and the total number of the nonzeros in the matrices \(D_{h,\text{FD}}\) and \(D_{h,\text{FE}}\) corresponding to the discretizations in (4) and (8), respectively. To visualize the nonzero entries of these matrices, we have depicted them in Figure 2.  
Figure 1: Two regular constellations of grid points: a fixed one \(\mathbf{x}_{0}\) and its neighbors. They are presented on a uniform rectangular grid (**left**) and on a uniform simplicial (or triangular) grid (**right**). Neighboring points are connected. On the right-hand side, \(h_{s}=\frac{5}{3}h_{x}\). Whenever the finite element discretization is flexible regarding the geometry of the domain, compared to the simple finite difference approximation given by \(D_{h,\text{FD}}\), it has some drawbacks as well:  
* We need more nonzero entries in the matrices \(D_{h,\text{FE}}\), as one can observe in \[\Delta u(\mathbf{x}_{0})\approx\frac{u(\mathbf{x}_{1})-2u(\mathbf{x}_{0})+u( \mathbf{x}_{3})}{h_{x}^{2}}+\frac{u(\mathbf{x}_{2})-2u(\mathbf{x}_{0})+u( \mathbf{x}_{4})}{h_{y}^{2}}, \tag{3}\]0 It enhances both the assembling time of the matrices and the solution of the corresponding linear systems. * We need to generate and store the entire structure of the mesh, including a list of triangles with their nodes and faces. * In a nonuniform grid, we need to transform the gradient of basis functions on each of the triangles according to the transformation of the reference triangles computing Jacobian matrices and their inverses (see [13]). To avoid all of these, our objective is to develop a finite difference approximation, which has the beneficial properties in (i)-(iii). A direct way for constructing such a second-order approximation is, however, impossible even in the one-dimensional case using two neighboring grid points. This is stated, e.g., in an early contribution [17]. One should also note that, assuming a quasi-uniform one-dimensional mesh, the second-order accuracy can be preserved. Here, the difference in the neighboring grid distances should be one magnitude less compared to the grid size. In a general case, a possible idea to obtain a second-order (or even higher-order) finite difference scheme without using a wide stencil (data from more neighboring grid points) is offered by the so-called method of compact difference schemes introduced in [18]. In this framework, the coefficients are not given explicitly; they have to be computed using a linear system. In the corresponding equalities, on the left hand side the linear combination of function values are given as in (1) or (3). At the same time, on the right hand side, again, a linear combination of the second-order derivatives arises. This idea was further developed; see, e.g., [19].