Journal of Conference Abstracts

Optimally Accurate Numerical Methods for Computing Synthetic Seismograms

Robert Geller (bob@global.geoph.s.u-tokyo.ac.jp), Nozomu Takeuchi (takeuchi@geoph.s.u-tokyo.ac.jp) & Hiromitsu Mizutani (mizutani@global.geoph.s.u-tokyo.ac.jp)

Dept. of Earth and Planetary Physics, Faculty of Science, Tokyo University, Bunkyo-ku, Tokyo, 113-0033, Japan

Research on numerical methods in almost all fields of geophysics is faced with the problem of developing accurate and efficient numerical algorithms. All methods (e.g. finite difference or pseudo-spectral) use some approximation of the exact partial differential operators. It is easy to quantify the error of these numerical approximations (e.g. O(z²) or O(z⁴)). However, what we really want to know is the error of the numerical solutions (e.g. XX%) for a given choice of grid spacing for some particular wavelength. For example, if all of the numerical operators are O(z²), the error of the numerical solution will also be proportional to z², but, for a given wavelength, what will be the proportionality constant of z²? This is what governs for the user whether or not a given level of approximation is acceptable. Standard numerical techniques for estimating, in advance, the actual magnitude of the error of the numerical solution do not appear to exist, but we have devised a such a method. A relatively thorough search of the literature has failed to turn up any similar technique; we would greatly appreciate being informed of previous work if known to the reader. In any case, however, the techniques discussed here have proven highly useful in modeling seismic wave propagation. This is a linear problem, and our techniques fully exploit this fact. It is not clear to what extent our techniques will prove applicable to non-linear problems.

We begin by deriving a method for formally estimating the error of a numerical solution computed using some given approximate differential operators. To do this we formally expand the error of the synthetic seismograms in terms of the normal modes (whose explicit eigenfrequencies and eigenfunctions need not be known) (Geller and Takeuchi, 1995). The eigenfunction expansion is the key to our derivation, as it allows us to map from the space of the numerical operators to the space of the numerical solution. We find that in order to obtain stable estimates of the error of the numerical solution the various numerical operators must be "tuned" so that their errors cancel to the greatest possible extent. We call the resulting numerical operators "optimally accurate" for a given numerical scheme (e.g. O(z²) finite difference) and a given grid spacing. The above approach is generally applicable. It can be applied to problems in essentially arbitrarily heterogeneous media, and to computations in either the time-domain or the frequency domain. It can be applied to both purely local numerical methods (e.g., finite difference or finite element) and also to methods using global sets of trial functions (e.g., the pseudo-spectral method or methods using expansions in terms of spherical harmonics).

The initial application of the above results was to frequency domain calculations of synthetics for Earth models in spherical coordinates using the Direct Solution Method (DSM), a Galerkin weak-form approach. Methods for 1-D (spherically symmetric) models were presented for toroidal and spheroidal wavefields by Cummins et al. (1994) and Takeuchi et al. (1996) respectively, and methods for 3-D (laterally heterogeneous) models were presented by Cummins et al. (1997). More recently we have applied the same general approach to develop optimally accurate 2nd order time-domain finite difference operators for 1-D problems (Geller and Takeuchi, 1998, in press). We will also present more recent results for optimally accurate time-domain schemes for both finite difference and pseudo-spectral operators in 2-D or 3-D media. In all cases an order of magnitude (or greater) improvement in cost-effectiveness over conventional numerical schemes of the same type is achieved (defined by the ratio of CPU times required to achieve some given level of accuracy) and there are no obvious disadvantages.

Prospects for extending our approaches for developing optimally accurate numerical operators to non-linear problems will be discussed.

Cummins et al., GRL, 21, 533-536 (1994)

Cummins et al., GJI, 130, 1-16, (1997)

Geller & Takeuchi, GJI, 123, 449-470, (1995)

Geller & Takeuchi, GJI, in press, (1998)

Takeuchi et al., GRL 23, 1175-1178, (1996)