# Efficient Discretization of the Laplacian: Application to Moving Boundary Problems
## 2 Materials and Methods
### Application to the Numerical Solution of the Stefan Problem
#### 2.5.1. STEP (i) (b)  
We now discuss step (i) (b). According to the implicit Euler time discretization, we have the general scheme  
\[\mathbf{u}^{n+1}=\mathbf{u}^{n}+\delta\cdot(\Delta_{h}\mathbf{u}^{n+1}+ \mathbf{f}^{n+1}), \tag{16}\]  
where the matrix \(\Delta_{h}\) contains data from the virtual outward points. To incorporate them into the scheme, we will use the given Neumann boundary data. To simplify the notations, the algorithm will be presented on the left side of the computational domain. Here, a generic interior node has the index \(K\) and its left virtual neighbors have the indices \(K-Nx-3,K-1\) and \(K+Nx+1\). Yielding the components of matrices and vectors in the subscript, and using the notation  
\[J =\{K-Nx-3,K-Nx-2,K-Nx-1,K-1,K,\] \[K+1,K+Nx+1,K+Nx+2,K+Nx+3\},\]  
we can express the component \(K\) of the product \(\Delta_{h}\mathbf{u}^{n+1}\) as  
\[(\Delta_{h}\mathbf{u}^{n+1})_{K}=\sum_{j\in J}(\Delta_{h})_{K,j}\mathbf{u}_{j} ^{n+1}. \tag{17}\]  
Here, we have to eliminate the sum  
\[(\Delta_{h})_{K,K-Nx-3}\mathbf{u}_{K-Nx-3}^{n+1}+(\Delta_{h})_{K,K-1}\mathbf{ u}_{K-1}^{n+1}+(\Delta_{h})_{K,K+Nx+1}\mathbf{u}_{K+Nx+1}^{n+1}. \tag{18}\]  
Using the second-order approximations  
\[\mathbf{u}_{K-Nx-3}^{n+1}-\mathbf{u}_{K-Nx-2}^{n+1}\approx\partial _{v}u((n+1)\delta,0,y_{K-1})\] \[\mathbf{u}_{K-1}^{n+1}-\mathbf{u}_{K}^{n+1}\approx\partial_{v}u( (n+1)\delta,0,y_{K})\] \[\mathbf{u}_{K+Nx+1}^{n+1}-\mathbf{u}_{K+Nx+2}^{n+1}\approx \partial_{v}u((n+1)\delta,0,y_{K+1}),\]  
we can rewrite (17) as  
\[(\Delta_{h}\mathbf{u}^{n+1})_{K}=(\Delta_{h})_{K,K-Nx-3}(\mathbf{ u}_{K-Nx-3}^{n+1}-\mathbf{u}_{K-Nx-2}^{n+1})+(\Delta_{h})_{K,K-1}(\mathbf{u}_{K-1 }^{n+1}-\mathbf{u}_{K}^{n+1})\] \[+(\Delta_{h})_{K,K+Nx+1}(\mathbf{u}_{K+Nx+1}^{n+1}-\mathbf{u}_{K +Nx+2}^{n+1})\] \[J =\{K-Nx-3,K-Nx-2,K-Nx-1,K-1,K,\]0 \[J =\{K-Nx-3,K-Nx-2,K-Nx-1,K-1,K,\]1 \[J =\{K-Nx-3,K-Nx-2,K-Nx-1,K-1,K,\]2 \[J =\{K-Nx-3,K-Nx-2,K-Nx-1,K-1,K,\]3 \[J =\{K-Nx-3,K-Nx-2,K-Nx-1,K-1,K,\]4 \[J =\{K-Nx-3,K-Nx-2,K-Nx-1,K-1,K,\]5 \[J =\{K-Nx-3,K-Nx-2,K-Nx-1,K-1,K,\]6 \[J =\{K-Nx-3,K-Nx-2,K-Nx-1,K-1,K,\]7such that we can eliminate the sum in (18). Accordingly, we perform the following:  
* on the right-hand side of (16), we have to add \(\delta\) times the first three terms in (19) modifying \(\mathbf{f}^{n+1}\);
* we have to modify the matrix elements \((\Delta_{h})_{K,K-Nx-2},(\Delta_{h})_{K,K}\) and \((\Delta_{h})_{K,K-Nx+2}\) by adding the entries \((\Delta_{h})_{K,K-Nx-3},(\Delta_{h})_{K,K-1}\) and \((\Delta_{h})_{K,K-Nx+1}\), respectively, corresponding to the virtual grid points as shown in the next three terms in (19). Additionally, we have to incorporate the Dirichlet boundary conditions. Regarding the common boundary, the situation is easy: we do not have to change anything except for leaving the virtual points located here. On the top of the top domain, we have to collect the Dirichlet data \(\mathbf{u}_{K+Nx+1}^{n+1},\mathbf{u}_{K+Nx+2}^{n+1},\mathbf{u}_{K+Nx+3}^{n+1}\) associated with the interior point with index \(K\). Then, on the right-hand side of (16), we have to keep  
\[J =\{K-Nx-3,K-Nx-2,K-Nx-1,K-1,K,\]8  
modifying with this the source term \(\mathbf{f}^{n+1}\) again, while the rest has to be moved to the left-hand side. The new source term is denoted with \(\tilde{\mathbf{f}}^{n+1}\). Performing a similar process on the bottom face of the bottom domain, we finally obtain the linear system  
\[J =\{K-Nx-3,K-Nx-2,K-Nx-1,K-1,K,\]9  
where  
* \(\tilde{\Delta}_{h}\) contains only entries corresponding to interior elements, which are modified according to the Neumann boundary data in (17). * Similarly, the new source term, according to the incorporation of [PERSON] and [PERSON] data, is denoted with \(\tilde{\mathbf{f}}^{n+1}\).