## 3 Theoretical Baseline
### Finite Difference Method
We thus need to evaluate the weak form in a way such that the resulting matrix \(\phi_{i}\cdot\phi_{j}\) becomes sparse. It turns out that this can be achieved by choosing a collocation method for quadrature. In that case, quadrature points are chosen to coincide with the node coordinates, and as a consequence, the mass matrix \(\phi_{i}\cdot\phi_{j}\) becomes the identity matrix. From the family of Gaussian quadrature schemes, one can achieve this using a Gauss-Lobatto quadrature of order equal to the polynomial order of the finite element. We note at this point that choosing this particular combination of quadrature method and number of nodes will lead to inexact integration and, thus, a numerical error is introduced. In this case, to produce a collocated scheme, one must pick the second-order Gauss-Lobatto variant using two quadrature nodes per coordinate direction. As this type of integration is known to be accurate up to degree \(2n-3\), this scheme will only integrate linear polynomials exactly. However, the above-listed basis functions are bilinear and have a combined polynomial order of 2. Thus, integration will not be accurate in this case. However, the modern FEM in general does not prescribe any particular method of integrating the weak formulation, per se [31]. Thus, although the results of a properly implemented FEM in the sense of exact integration will be slightly more accurate, one can still show equivalence regarding a particular instance of the FEM. We now evaluate the element-wise stiffness matrix \(-\int_{\Omega^{(i)}}\
abla_{k}\phi_{i}\
abla_{k}\phi_{j}\,dx\) within the reference domain \([0;1]\times[0;1]\) for the given first order Lagrange element using Gauss-Lobatto quadrature, more specifically the variant using two quadrature points per coordinate direction. This results in  
\[\underline{\underline{K}}^{(e)}= \begin{bmatrix}-1&1/2&1/2&0\\ 1/2&-1&0&1/2\\ 1/2&0&-1&1/2\\ 0&1/2&1/2&-1\end{bmatrix} . \tag{14}\]  
As such, \(\underline{\underline{K}}^{(e)}\) does not yet equal Equation (8). The final step consists of assembling the linear system in the physical domain using the reference stiffness matrix. In a Cartesian mesh in two dimensions, an interior node is owned by four quadrilateral elements. If one carries out this assembly process, an equivalent formulation can be obtained:  
\[\underline{\underline{K}}= \begin{bmatrix}\ddots&\ddots&&&\vdots&\vdots&\vdots \\ \ddots&\ddots&\ddots&&\vdots&1&\vdots\\ &\ddots&\ddots&\ddots&\vdots&\vdots&\vdots\\ &&\ddots&\ddots&\ddots&1&\vdots\\ \ldots&\ldots&\ldots&\ddots&\ddots&\ddots&\vdots&\ldots&\ldots&\ldots\\ \ldots&1&\ldots&1&\ddots&-4&\ddots&1&\cdots&1&\ldots\\ \ldots&\ldots&\ldots&\ldots&\vdots&\ddots&\ddots&\ddots&\ldots&\ldots\\ &&&\vdots&1&\ddots&\ddots&\ddots&&\\ &&&\vdots&\vdots&\vdots&\ddots&\ddots&\ddots&\\ &&&\vdots&1&\vdots&&\ddots&\ddots&\ddots\\ &&\vdots&\vdots&\vdots&&&\ddots&\ddots&\ddots\\ &&\vdots&\vdots&\vdots&&&\ddots&\ddots&\ddots\\ \end{bmatrix} . \tag{15}\]  
The exact position of the one entry in the typically large and sparse matrix depends on the mesh topology as well as the global numbering of the degrees of freedom. For the discretization of other operators, a similar argument holds, as the shown procedure is irrespective of the choice of weak form or basis function. For example, one could discretize the gradient of a function \(\
abla u\) using an upwind finite difference formulation in fluid mechanics for resolving convective terms. An equivalent finite element method can be assembled by producing a weak form, as given in the above example, choosing the same collocation method and carrying out the integration numerically. However, one important difference is that one cannot choose the test space to be equivalent to the trial space. This would yield a symmetric system that does not correspond to an upwind finite difference formulation and is also not stable in the case of solving a pure advection equation.