## 2 Materials and Methods
### Problem Statement and Mathematical Background  
We investigate the numerical approximation  
\[\sum_{j=0}^{l}a_{j}\mu(\mathbf{x}_{j})\approx\Delta u(\mathbf{x}_{0}) \tag{1}\]  
of the Laplacian operator at a given point \(\mathbf{x}_{0}\). Here, \(\{\mathbf{x}_{j}\}_{j=1}^{J}\) denotes the adjacent grid points to \(\mathbf{x}_{0}\) with given function values \(\{u(\mathbf{x}_{j})\}_{j=1}^{J}\) and \(u(\mathbf{x}_{0})\), respectively. In practice, the coefficients \(\{a_{j}\}_{j=0}^{J}\) are to be optimized, which depend on the position of the grid points. Indeed, we apply this approximation for a spatial discretization of a domain \(\Omega\), which can be regarded as a grid. The aim of such an approximation is the numerical solution of some PDE. To demonstrate the principles, we consider the simple [PERSON]'s problem  
\[ \begin{cases}\Delta u(\mathbf{x})=f(\mathbf{x})&\mathbf{x}\in\Omega\\ u(\mathbf{x})=g(\mathbf{x})&\mathbf{x}\in\partial\Omega,\end{cases}  \tag{2}\]  
where \(\Omega\subset\mathbb{R}^{2}\) denotes the computational domain, \(u:\Omega\to\mathbb{R}\) is the unknown function and \(g:\partial\Omega\to\mathbb{R}\) is given. The well-posedness of (2) is ensured in \(H^{1}(\Omega)\) provided that \(\Omega\) is a bounded Lipschitz domain, \(f\in H^{-1}(\Omega)\) and \(g\in H^{1}(\partial\Omega)\) (see [16]). Concerning the numerical solution of (2), in the most easy case, using a uniform rectangular grid, we have locally the geometric setup shown on the left in Figure 1. Here, for the approximation of \(\Delta u(\mathbf{x}_{0})\), we use the function values at \(\mathbf{x}_{1},\mathbf{x}_{2},\mathbf{x}_{3}\) and \(\mathbf{x}_{4}\), which are considered the neighbors of \(\mathbf{x}_{0}\). Then, in a given grid point \(\mathbf{x}_{0}\), we have the classical five-point approximation  
\[\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}\]  
so that in this case, in (1), we will have  
\[a_{0}=-\frac{2}{h_{x}^{2}}-\frac{2}{h_{y}^{2}},\ \ a_{1}=a_{3}=\frac{1}{h_{x}^{ 2}}\quad\text{and}\quad a_{2}=a_{4}=\frac{1}{h_{y}^{2}}. \tag{4}\]  
This has the approximation order 2, i.e., for the difference of the two sides in (3), the choice in (4) delivers the approximation order 2 with respect to both variables:  
\[\Delta u(\mathbf{x}_{0})-\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}}=\mathcal{O}(h_{x}^{2})+\mathcal{O}(h_{y}^{2}), \tag{5}\]  
provided that the function \(u\) is four-times continuously differentiable. Whenever it seems to be a strict assumption, for a number of PDEs (such as diffusion problems), including the Laplacian operator, we have smooth analytic solutions. In practice, the pointwise approximations in (3) are assembled into a matrix \(D_{h,\text{FD}}\). Incorporating also the boundary condition and using the vector \(\mathbf{f}_{h}\) for the pointwise values of \(f\) in the grid points, we finally have a linear algebraic problem to solve:  
\[D_{h,\text{FD}}\mathbf{u}_{h,\text{FD}}=\mathbf{f}_{h} \tag{6}\]  
The solution \(\mathbf{u}_{h,\text{FD}}\) is then called a finite difference (FD) approximation of \(u\) in (2). In another well-known case, one can discretize any domain into (possibly curved) triangles and, based on this tessellation, a first-order finite element (FE) approximation can be applied. This situation in a special geometric case is presented on the right in Figure 1. Here, the neighbors of \(\mathbf{x}_{0}\) are the points \(\mathbf{x}_{1},\mathbf{x}_{2},\ldots,\mathbf{x}_{6}\). In this case, instead of the direct form of the Laplacian operator, its weak form is discretized.