In this paper, we adopt the CIP (Cubic Interpolated Profile) method to the simulation of elastic waves as a highly accurate and stable algorithm to solve characteristic equations. The basic computation schema was developed by Yabe and Aoki (1991) and had been used in the computational fluid dynamics and plasmas phenomena simulations. The key idea of the CIP is that not only the physical value itself but also its first spatial derivative obeys the same characteristic equations. Using this property the solution is interpolated by cubic polynomials and interpolation coefficients can be evaluated arithmetically.
We implemented this idea to elastic wave simulation by derivation of the characteristic equations from the basic equations of motion. The characteristic equations can be interpreted as the combined one-way wave equations for each mode of wave. Then, we define boundary conditions by using the merit of one-side propagation; free surface, solid-fluid boundary, irregular topography. From simulation study and the stability evaluation, we recognize the method of characteristics with the CIP is a very powerful simulation technique for the elastic wave propagation in geophysics. Numerical dispersion is negligible, requiring about half the number of grid cells per wavelength than other solvers. This allows accurate, high-frequency, full-wavefield simulation in models with highly variable, random elastic contrasts with fluid-solid mixed media and complex topographic media.