# Efficient Discretization of the Laplacian: Application to Moving Boundary Problems
## 2 Materials and Methods
### The Pointwise Optimization  
To achieve possibly the highest-order approximation in (2) and (3), these should hold with equality for some polynomials until a certain order.  
For example, if we have a uniform mesh with the grid size \(h\) and the coefficients are scaled as \(h^{-2}\), then to obtain an approximation order \(N-1\), the equality  
\[\Delta p(\mathbf{x}_{0})=a_{0}p(\mathbf{x}_{0})+a_{1}p(\mathbf{x}_{1})+\dots+ a_{8}p(\mathbf{x}_{8}). \tag{7}\]  
should hold for all polynomials \(p\) up to order \(N\). Henceforth, we focus to this approximation with nine grid points.  
At the same time, for general mesh geometries, the approximation order is decreased and, accordingly, (7) is not valid any more for all polynomials up to order \(N\). The principal problem here is to determine them, which should be taken into account in (3). Also, we demand that the right-hand side of (3) is zero for any constant function. This immediately gives the constraint  
\[a_{0}=-(a_{1}+a_{2}+\dots+a_{8})\]  
for the coefficients in (3). To satisfy (7), we will first consider the set of polynomials  
\[\begin{split} p_{1}(x,y)&=100x,\ p_{2}(x,y)=100y,\ p_{3}(x,y)=10 xy,\\ p_{4}(x,y)&=10x^{2},\ p_{5}(x,y)&=10( x^{2}-y^{2}),\end{split} \tag{8}\]  
which span the space of first- and second-order ones. Additionally, we take  
\[\begin{split} p_{6}(x,y)&=x^{3}-3x^{2}y,\quad p_{7} (x,y)&=y^{3}-3 xy^{2},\\ p_{8}(x,y)&=x^{3}y-xy^{3},\quad p_{9}(x,y)& =x^{4}+6x^{2}y^{2}-y^{4},\\ p_{10}(x,y)&=5x^{4}y-10x^{2}y^{3}+y^{5},\quad p_{11 }(x,y)&=x^{5}-10x^{3}y^{2}+5 xy^{4},\end{split}\]  
Figure 2: Location of non-zero elements in the discretization matrix \(D_{h,\text{FD}}\) corresponding to (4) (**left**) and (5) (**right**), respectively. The total matrix size \(36\times 36\) corresponds to the number of the interior grid points. The notation nz yields the number of non-zeros.  
which are the harmonic polynomials up to order five. To motivate this choice, we recall that in a number of real-life cases, we have to compute with harmonic functions, such that \(\Delta u=0\). In other words, we can only have a certain accuracy for approximating the solution of \(\Delta u=f\), if we can perform it in the case of \(f=0\). In this way, an approximation in (5) has the form  
\[a_{1}p_{j}(\mathbf{x}_{1})+a_{2}p_{j}(\mathbf{x}_{2})+\dots+a_{8}p_{j}( \mathbf{x}_{8})\approx\Delta p_{j}(\mathbf{x}_{0}),\quad j=1,2,\dots,11. \tag{9}\]  
Performing it for all polynomials \(p_{1},p_{2},\dots,p_{11}\) above, for the geometrical setup given with \(\mathbf{x}_{0},\mathbf{x}_{1},\dots,\mathbf{x}_{8}\), we obtain the following optimization task:  
* Find the vector of coefficients \(\mathbf{a}=[a_{1},a_{2},\dots,a_{8}]^{T}\) in (9) such that \(P\cdot\mathbf{a}\approx\Delta\mathbf{p}\), where
* \(P\in\mathbb{R}^{11\times 8}\) with \(P_{jk}=p_{j}(\mathbf{x}_{k})\);
* \(\Delta\mathbf{p}=[0,0,0,20,0,0,0,0,0,0]^{T}\in\mathbb{R}^{11}\), so that \(\Delta\mathbf{p}_{j}=p_{j}(\mathbf{x}_{0})\).  
Using the coefficient \(100\) in the first-order polynomials \(p_{1}\) and \(p_{2}\) and the coefficient \(10\) in the second-order ones \(p_{3},p_{4}\) and \(p_{5}\) in (8), we force more accuracy in the case of the lower-order polynomials. This procedure significantly improves the approximation in (9).  
The optimization task in (9) is solved in the least-square sense for each point \(\mathbf{x}_{0}\) using the np_linalg_lstsq subroutine in Python. This can be performed for all grid points independently. In this way, we can easily vectorize the procedure in the following way:  
* We create a matrix, where each row contains the geometry around a fixed grid point.
* The optimization procedure is applied (independently) to the rows of this matrix using the np_apply_along_axis subroutine in Python. Accordingly, this final step could also be executed with parallel processing.