  An FFT-based algorithm is developed to simulate the propagation of elastic
waves in heterogeneous $d$-dimensional rectangular shape domains. The method
allows one to prescribe the displacement as a function of time in a subregion
of the domain, emulating the application of Dirichlet boundary conditions on an
outer face. Time discretization is performed using an unconditionally stable
beta-Newmark approach. The implicit problem for obtaining the displacement at
each time step is solved by transforming the equilibrium equations into Fourier
space and solving the corresponding linear system with a preconditioned Krylov
solver. The resulting method is validated against analytical solutions and
compared with implicit and explicit finite element simulations and with an
explicit FFT approach. The accuracy of the method is similar to or better than
that of finite elements, and the numerical performance is clearly superior,
allowing the use of much larger models. To illustrate the capabilities of the
method, some numerical examples are presented, including the propagation of
planar, circular, and spherical waves and the simulation of the propagation of
a pulse in a polycrystalline medium.
