# "Cleanest, most professional" approach to implement a finite differences scheme in dimension greater than two

Coding a finite difference algorithm in 1D does not require a complex mesh. In higher dimensions, you would need a mesh and its connectivity graph to compute the differential operators. Finite elements can be seen as a generalization of finite differences.

So my question is:

How is a finite difference code implemented? Are finite difference methods implemented with finite element algorithms on a Cartesian mesh using a product of P1 basis functions?

• I wouldn't consider finite elements to be a generalization of finite differences, the formulation of the weak problem is quite unique and not at all how finite differences operates. High order finite volume methods might be a more apt comparison to high order finite differences. There is the K-Exact method which generalizes high order finite volume methods to arbitrary order and unstructured meshes, though it still suffers from the fundamental problem of wide stencil methods that you have to figure out what to do at the boundaries. Commented May 10 at 16:08
• and it is possible to treat low order finite volume methods as equivalent to a piecewise constant discontinuous Galerkin method in certain situations. Commented May 10 at 16:11
• If it is a Cartesian mesh, the basis functions are tensor products of P1 basis functions (e.g. nodal basis), and if the continuity is enforced (e.g. if it is not a DG method), then yes, FEM and FDM will give you exactly the same matrix, the same right hand side and the same solution for all linear problems Commented May 12 at 2:29
• Implementing a FDM via FEM assembly seems a bit overkill as the primary advantage of FDM is that the operators and connectivity are known ahead of time so that the discretization can be chosen to have advantageous structure. Commented May 12 at 17:31

• @mle Classical FDM locally fits an interpolating polynomial and takes the derivatives of it. You can prove that this is the same result as from the Taylor expansion. As far as I am aware the reason it is called finite differences is because some derivatives can be expressed as differences between the point values in the tensor product grid setting, e.g. $d_x f(x) = \frac{f(x+h)-f(x)}{h} + O(h)$. But the above is just a special case of taking the derivative of a linear function interpolating $(x,f(x))$ and $(x+h, f(x+h))$. More generally you can take higher degree polynomials and more points. Commented May 30 at 10:19