Optimally accurate operators for partial differential equations

In this contribution, we generalize the concept of \textit{optimally accurate operators} proposed and used in a series of studies on the simulation of seismic wave propagation, particularly based on Geller \& Takeuchi (1995). Although these operators have been mathematically and numerically proven to be more accurate than conventional methods, the theory was specifically developed for the equations of motion in linear elastic continuous media. Furthermore, the original theory requires compensation for errors from each term due to truncation at low orders during the error estimation, which has limited its application to other types of physics described by partial differential equations. Here, we present a new method that can automatically derive numerical operators for arbitrary partial differential equations. These operators, which involve a small number of nodes in time and space (compact operators), are more accurate than conventional ones and do not require meshing. Our method evaluates the weak formulation of the equations of motion, developed with the aid of Taylor expansions. We establish the link between our new method and the classic optimally accurate operators, showing that they produce identical coefficients in homogeneous media. Finally, we perform a benchmark test for the 1D Poisson problem across various heterogeneous media. The benchmarks demonstrate the superiority of our method compared to conventional operators, even when using a set of linear B-spline test functions (three-point hat functions). However, the convergence rate can depend on the wavelength of the material property: when the material property has the same wavelength as that of the field, the convergence rate is O(4), whereas it can be less efficient O(2) for other models.