Wave atoms and time upscaling of wave equations

• Authors: Laurent Demanet;Lexing Ying
• Venue: Numerische Mathematik
• Year: 2009

We present a new geometric strategy for the numerical solution of hyperbolic wave equations in smoothly varying, two-dimensional time-independent periodic media. The method consists in representing the time-dependent Green’s function in wave atoms, a tight frame of multiscale, directional wave packets obeying a precise parabolic balance between oscillations and support size, namely wavelength ~(diameter).2 Wave atoms offer a uniquely structured representation of the Green’s function in the sense that the resulting matrix is universally sparse over the class of C ∞ coefficients, even for “large” times; the matrix has a natural low-rank block-structure after separation of the spatial indices. The parabolic scaling is essential for these properties to hold. As a result, it becomes realistic to accurately build the full matrix exponential in the wave atom frame, using repeated squaring up to some time typically of the form $${\Delta t \sim \sqrt{\Delta x}}$$, which is bigger than the standard CFL timestep. Once the “expensive” precomputation of the Green’s function has been carried out, it can be used to perform unusually large, upscaled, “cheap” time steps. The algorithm is relatively simple in that it does not require an underlying geometric optics solver. We prove accuracy and complexity results based on a priori estimates of sparsity and separation ranks. On a N-by-N grid, the “expensive” precomputation takes somewhere between O(N 3log N) and O(N 4log N) steps depending on the separability of the acoustic medium. The complexity of upscaled timestepping, however, beats the O(N 3log N) bottleneck of pseudospectral methods on an N-by-N grid, for a wide range of physically relevant situations. In particular, we show that a naive version of the wave atom algorithm provably runs in O(N 2+δ) operations for arbitrarily small δ—but for the final algorithm we had to slightly increase the exponent in order to reduce the large constant. As a result, we get estimates between O(N 2.5 log N) and O(N 3 log N) for upscaled timestepping. We also show several numerical examples. In practice, the current wave atom solver becomes competitive over a pseudospectral method in regimes where the wave equation should be solved hundreds of times with different initial conditions, as in reflection seismology. In academic examples of accurate propagation of bandlimited wavefronts, if the precomputation step is factored out, then the wave atom solver is indeed faster than a pseudospectral method by a factor of about 3–5 at N = 512, and a factor 10–20 at N = 1024, for the same accuracy. Very similar gains are obtained in comparison versus a finite difference method.