Optimally accurate operators for arbitrary PDEs using interpolated Taylor expansions

We present a modified version of optimally accurate operators for partial differential equations of arbitrary order and arbitrary dimension. Optimally accurate operators were originally proposed for seismic wave propagation in homogeneous media by Geller and Takeuchi (1995), who derived compact operator coefficients for specific wave equations in specific dimensions. Fuji and Duretz (2025) showed that these coefficients can be obtained by formulating a weak form of the PDE using basis functions defined as Taylor expansions about the centred grid point, yielding a reduction of the error by a factor of 100 for the 1D Poisson equation with three collocated grid points. However, the convergence rate varies from O4 to O2 depending on the degree of heterogeneity. Here, we generalise the theory by using Taylor expansions at all grid points involved as basis functions. The construction of symbolic expressions for the resulting coefficients requires nested loops over all grid points in space and time, which becomes intractable without GPU acceleration. In this contribution, we present the theoretical framework, benchmark results, and a Julia notebook implementing the proposed method.