Multistep and Multistage Convolution Quadrature for the Wave Equation: Algorithms and Experiments

We describe how a time-discretized wave equation in a homogeneous medium can be solved by boundary integral methods. The time discretization can be a multistep, Runge–Kutta, or a more general multistep-multistage method. The resulting convolutional system of boundary integral equations belongs to the family of convolution quadratures of Lubich. The aim of this work is twofold. It describes an efficient, robust, and easily parallelizable method for solving the semidiscretized system. The resulting algorithm has the main advantages of time-stepping methods and of Fourier synthesis: at each time-step, a system of linear equations with the same system matrix needs to be solved, yet computations can easily be done in parallel; the computational cost is almost linear in the number of time-steps; and only the Laplace transform of the time-domain fundamental solution is needed. The new aspect of the algorithm is that all this is possible without ever explicitly constructing the weights of the convolution quadrature. This approach also readily allows the use of modern data-sparse techniques to efficiently perform computation in space. We investigate theoretically and numerically to which extent hierarchical matrix ($\mathcal H$-matrix) techniques can be used to speed up the space computation. The second aim of this article is to perform a series of large-scale 3D experiments with a range of multistep and multistage time discretization methods: the backward difference formula of order 2 (BDF2), the Trapezoid rule, and the 3-stage Radau IIA methods are investigated in detail. One of the conclusions of the experiments is that the Radau IIA method often performs overwhelmingly better than the linear multistep methods, especially for problems with many reflections, yet, in connection with hyperbolic problems, BDFs have so far been predominant in the literature on convolution quadrature.