Optimally accurate second-order time-domain ¢nite di¡erence scheme for the elastic equation of motion: one-dimensional case

SUMMARY We previously derived a general criterion for optimally accurate numerical operators for the calculation of synthetic seismograms in the frequency domain (Geller & Takeuchi 1995). We then derived modi¢ed operators for the Direct Solution Method (DSM) (Geller & Ohminato 1994) which satisfy this general criterion, thereby yielding signi¢cantly more accurate synthetics (for any given numerical grid spacing) without increasing the computational requirements (Cummins et al . 1994; Takeuchi, Geller & Cummins 1996; Cummins, Takeuchi & Geller 1997). In this paper, we derive optimally accurate time-domain ¢nite di¡erence (FD) operators which are second order in space and time using a similar approach. As our FD operators are local, our algorithm is well suited to massively parallel computers. Our approach can be extended to other methods (e.g. pseudo-spectral) for solving the elastic equation of motion. It might also be possible to extend this approach to equations other than the elastic equation of motion, including non-linear equations.


INTRODUCTION
Waveform inversion (e.g. Tarantola 1984;Geller & Hara 1993) is a promising approach for determining 3-D Earth structure, as it directly uses all of the information in the recorded seismograms. Practical applications of waveform inversion require e¤cient methods for computing synthetic seismograms. It is also important to be able to quantify the accuracy of the synthetics. Reliable advance estimates of the accuracy that will be attained for a given numerical scheme and grid size are particularly desirable.
Standard techniques of numerical analysis provide estimates of the error of discretized numerical operators rather than the error of the numerical solutions (i.e. synthetic seismograms) computed using these operators. For example, the error of the three-point central ¢nite di¡erence (FD) operator for the second derivative is well known to be the second term of the following Taylor series: 1 *z 2 [u(z{*z){2u(z)zu(zz*z)]~d 2 u dz 2 z *z 2 12 The error of synthetic seismograms computed by schemes based on eq. (1) will thus also be proportional to *z 2 . However, this information is of little practical value unless the coe¤cient of *z 2 is also known. Geller & Takeuchi (1995), cited hereafter as GT95, developed a method that uses an eigenfunction expansion to make formal estimates of the error of synthetic seismograms computed using a given numerical scheme. Their results allow numerical operators to be`tuned' to produce optimally accurate numerical schemes of a given type and order. GT95 obtained the following estimate of the relative error of synthetics computed using an optimally accurate 1-D second-order FD scheme in the frequency domain: relative error& k 2 z *z 2 12 , much CPU time. Many other numerical schemes for computing synthetics, e.g. fourth-order FD in space and second-order FD in time (Levander 1988) or pseudo-spectral in space and second-order FD in time (Fornberg 1987), have also been presented. Mizutani, Geller & Takeuchi (1997) used basically the same method that we present in this paper to derive an optimally accurate scheme that is pseudo-spectral in space and second-order FD in time.
The question which method is the best for a given problem involves not only CPU time, but also memory requirements and parallelization, and is beyond the scope of this paper. Also, optimally accurate algorithms should be derived for each type of scheme before comparisons are made.
In this paper we derive modi¢ed time-domain FD operators using the results of GT95. The modi¢ed FD operators lead to an implicit scheme; that is, the solution of a set of simultaneous linear equations is required at each time step to obtain the displacements at the next time step. To obviate this time-consuming computation, we approximate the implicit scheme by a two-part explicit (`predictor^corrector') scheme. We present results for the 1-D case in this paper. Results for 2-D and 3-D cases in Cartesian and curvilinear coordinates will be presented in future work. The derivations in this paper are not self-contained; results obtained by GT95 are in general used here without proof.
In our frequency domain calculations of synthetics we routinely include anelastic attenuation in the elastic moduli using a superposition of standard linear solids (Liu, Anderson & Kanamori 1976). For simplicity we consider a purely elastic medium in this paper. It appears possible to use an approach similar to that of Emmerich & Korn (1987) to extend our methods to the anelastic case, but we have not yet veri¢ed this.
In many FD implementations (e.g. Vireux 1986) velocity and stress are the dependent variables, and a coupled system of ¢rstorder di¡erential equations in space and time is obtained. This system is then solved using the`staggered grid' approach. In contrast, displacement is the only dependent variable in both the conventional and modi¢ed operators in this paper. As discussed below in Section 6, it does not appear possible to derive modi¢ed operators for staggered grid schemes.

FINITE DIFFERENCE OPERATORS
The strong form of the time-domain equation of motion for the 1-D case is as follows: where u is the displacement, o is the density, k the elastic modulus and f the external force. Following the normal FD approach, we discretize the unknown displacement u, which is a function of space and time, as follows: The FD equation of motion is as follows: where A is an FD operator for o(L 2 u/Lt 2 ), K is an FD operator for L/Lz[k(Lu/Lz)] and f is the discretized external force term. Summation over repeated indices is implied throughout this paper. A and K in eq. (5) are time-domain FD operators. We now compute their Fourier transforms, which are denoted by B and L respectively. In the frequency domain, the displacement is discretized as follows: The FD equation of motion in the frequency domain (the transform of eq. 5) is as follows: As eq. (7) and eq. (2.1) of GT95 have essentially the same form, we can use the theory in Section 2 of GT95 to obtain formal error estimates for d. The only di¡erence is that the errors of the operators are now u-dependent as well as k-dependent, because of numerical dispersion in time as well as space.

Homogeneous 1-D problem
We begin by considering a homogeneous 1-D problem with a constant temporal grid interval *t and a constant spatial grid interval *z. In this paper the superscript 0 denotes the conventional operators rather than the exact operators. The conventional FD operators A 0 and K 0 are as follows: with the format Throughout this paper blank spaces in the FD stencils denote zeros.
A 0 and K 0 have the numerical dispersion of normal three-point second-derivative operators. The operator error for the conventional operators is thus as follows: Note that throughout this paper the error expressions are given to O(*t 2 ) and O(*z 2 ), with higher-order terms omitted. We de¢ne B 0 , L 0 , B exact , L exact , äB 0 , and äL 0 to be the Fourier transforms of A 0 , K 0 , A exact , K exact , äA 0 , and äK 0 respectively. These operators are related as follows: Transforming eq. (10) into the frequency domain, we obtain The quantity on the rhs of eq. (15) is the`basic error' of the operators (see GT95). When the operand is an eigenfunction and the frequency approaches the corresponding eigenfrequency, the basic error given by eq. (15) will not in general equal zero. Thus the conventional operators do not in general satisfy eq. (2.20) of GT95, and are therefore not optimally accurate. To satisfy eq. (2.20) of GT95, the basic error should be zero when u is an eigenfunction and u is equal to the corresponding eigenfrequency. We therefore must derive modi¢ed operators L mn and B mn for which, rather than eq. (15), the basic error is instead given by As the quantity inside the square brackets in eq. (16) is the lhs of the exact equation of motion, it is zero for every eigenfunction when u is equal to the corresponding eigenfrequency, thus the basic error is zero. By taking the inverse Fourier transform, the time-domain A~t representation of eq. (16) is obtained: ( 1 7 ) where A and K are the modi¢ed operators and äA and äK are their errors. It is easier to construct the numerical operators using the rhs of the ¢rst line of eq. (17), but the meaning of eq. (17) is more clearly shown by the second line. The latter shows clearly that the basic error of the modi¢ed operators is given by derivatives of the homogeneous equation of motion (the term in square brackets). Thus when u(z, t)~u m (z) exp (iu m t), where u m (z) is the eigenfunction of a mode with eigenfrequency u m , the bracketed term (and hence its derivatives) will be zero.
Omitting the details of the derivation, the modi¢ed operators that yield errors of the form of eq. (17) are as follows: Note that if we sum horizontally for A and sum vertically for K we obtain the conventional operators in eq. (8). An intuitive explanation of eq. (18) is that we smear out the discretized second time-derivative operator in space, and smear out the discretized second spatial-derivative operator in time, so that the numerical dispersion (error of the phase velocity) of the discretized equation of motion is zero to second order in *t 2 and *z 2 . We omit discussion of the boundary error in this paper because, as shown in Sections 2 and 3 of GT95, the boundary error does not have an important e¡ect on the error of the numerical solution.

Heterogeneous 1-D problem
The conventional and modi¢ed time-domain operators for an inhomogeneous medium can be derived using the procedures given in Section 3, eqs (3.48)^(3.52), of GT95. A ¢rst-order discontinuity in elastic properties is handled by overlapping (see Fig. 3 of GT95). Details are not given here.
The explicit forms of the conventional operators A 0 and K 0 for an inhomogeneous medium are as follows: where o m and k m are respectively the density and rigidity at the mth node. The spatial variation of the elastic properties is assumed to be reasonably smooth. The explicit forms of the modi¢ed operators A and K are as follows: The de¢nition of A in eq. (20) corresponds to T' right as de¢ned in eq. (3.49) of Geller & Takeuchi (1995). We could also have de¢ned A using T' left , or any other linear combination cT' right z(1{c)T' left .
Omitting details (see GT95) the basic error for the modi¢ed operators in eq. (20) is ( 2 1 )

Operators for boundaries
For completeness we give the operators for a left-hand boundary with a free-surface natural boundary condition. If this is an internal boundary rather than a free boundary, the operators are overlapped with those in the adjacent segment (see GT95 for details).
Operators for a right-hand boundary can be readily obtained from those given below. The boundary terms for the conventional operators for a heterogeneous medium are The modi¢ed operators for a left-hand boundary for a heterogeneous medium are

PREDICTOR^CORRECTOR SCHEME USING THE MODIFIED OPERATORS
The modi¢ed operator (A mMnN {K mMnN ) given by eqs (20) and (23) has multiple non-zero elements for time tz*t. If we use these modi¢ed operators in a time-marching scheme to solve eq. (5), the FD equation of motion, we obtain an implicit scheme, rather than the explicit scheme for the conventional operators in eq. (8). To obviate the need to solve a system of simultaneous linear equations at each time step, we use the modi¢ed operators in an approximate (predictor^corrector) scheme based on the ¢rst-order Born approximation. A detailed discussion of the implementation is given in the Appendix. First we predict the wave¢eld at the next time step using the conventional operators A 0 and K 0 de¢ned in eq. (8): where c 0 n(Nz1) , the predicted wave¢eld at time tz*t, is obtained explicitly from eq. (24). Next we compute äc, the correction to the displacement at time tz*t, using the ¢rst-order Born approximation. In the previous section we used äA and äK to denote the di¡erence between the numerical and exact operators. However, in this and later sections we denote the di¡erence between the conventional operators A 0 , K 0 and the modi¢ed operators A, K by äA, äK respectively. To obtain the correction we thus solve As the lhs of eq. (25) uses the conventional operators, we can solve explicitly for the value of äc at time tz*t. Note that dc nN~0 and dc n(N{1)~0 in eq. (25). We compute the corrected displacement c n(Nz1 ) after each time step using c 0 computed by eq. (24) and äc computed by eq. (25): Finally, before advancing to the next time step we rede¢ne c 0 : Note that we use the displacements given by eq. (27) as the values for c 0 nN and c 0 n(N{1) when we evaluate eq. (24) for the next time step. äA and äK are obtained by taking the di¡erence of eqs (20) and (19): äA and äK for the boundary elements are obtained by di¡erencing eqs (23) and (22).

Stability
It is well known that the time step, *t, must be less than or equal to the Courant limit, *t courant , in order to obtain stable solutions using the conventional operators. In this section we show that the Courant stability limits for the conventional and modi¢ed operators are equal for a homogeneous medium.
The stability condition for both modi¢ed and conventional operators is the condition that there cannot be any exponentially growing modes in the FD equation of motion (eq. 5). To derive the stability condition we solve a generalized eigenvalue problem (e.g. Section 7.7 of Golub & Van Loan 1989). We begin by formulating the stability condition for the conventional operators. We consider harmonic solutions of eq. (5) with spatial dependence c and temporal dependence exp (iut). We substitute this solution into eq. (5) to ( 2 9 ) obtain the following:6.5pt Hc~2 *t 2 (1{ cos u*t)Tc , where T and H are the spatially dependent parts of A and {K. T and H for a homogeneous medium are as follows: ( 3 1 ) We de¢ne eigenvalues j n and eigenvectors c n (Note that in this section n refers to the eigenvector and eigenvalue rather than to a particular node.) for T and H as follows: Hc n~jn Tc n , where 0¦j 0`j1`Á Á Á`j max . Comparing eqs (30) and (32) we get 2 *t 2 (1{ cos u*t)~j n or cos u*t~1{ *t 2 2 j n .
Since j n is real and j n §0 for all n, the requirement for real roots is cos u*t §{1 .
The most severe constraint is imposed by j n~jmax . Combining eqs (34) and (35), we thus require The eigenvectors of eq. (32) for the matrices for a homogeneous medium (eq. 31) are H where N is the number of grid intervals. The corresponding eigenvalues are given by where b~ k/o p is the wave velocity. We thus see from eqs (39) and (37) that From eqs (37) and (40) the stability condition is where Thus the Courant stability condition holds rigorously for the conventional operators in a 1-D homogeneous medium with free boundaries. We now consider the stability condition for the modi¢ed operators for a homogeneous medium. The spatial operator T' corresponding to the modi¢ed operator A is given by T'~o 5/12 1/12 1/12 10/12 1/12 F F F F F F F F F 1/12 10/12 1/12 1/12 5/12 The generalized eigenvalue problem for H and T' is given by where c n is given by eq. (38). We use the generalized eigenvalue problem for T' and T as an intermediate result: T'c n~cn Tc n , where c n~5 6 z 1 6 cos nn N and c n is given by eq. (38). Using eqs (39) and (44) We now consider harmonic solutions of eq. (5) for a homogeneous medium. We obtain 5 6 z 1 6 cos u*t Hc n~2 *t 2 (1{ cos u*t)T'c n .
As the condition for stability is cos u*t j j ¦{1, we obtain from eq. (52) Thus both the modi¢ed and conventional operators must satisfy the same stability condition.
j max can be determined for heterogeneous media by numerically ¢nding the maximum eigenvalue of eq. (32) using the Sturm sequence property (Peters & Wilkinson 1969). We can then obtain the stability condition for the inhomogeneous case following the same procedures. In general, the stability condition for an inhomogeneous medium is for both conventional and modi¢ed operators, where is a ¢nite but relatively small number whose amplitude and sign depend on the exact nature of the problem. We omit further details.

Accuracy
The analytic eigenvalues for a homogeneous 1-D medium are the squares of the eigenfrequencies, where L~N*z is the length of the medium. If we expand cos (nn/N) in a Taylor series, we obtain the following from eqs (39) and (47) respectively: Note that nn/N~k*z, where k is the wavenumber. Thus we see that the eigenfrequencies for the modi¢ed operators have an error of O(k 4 *z 4 ), while the eigenfrequencies for the conventional operators have an error of O(k 2 *z 2 ). This is an expected result, because, as discussed by GT95, setting the basic error of the modi¢ed operators to zero is equivalent to requiring the error of the eigenfrequencies of the modi¢ed operators to be zero to O(*z 2 ), excepting possible minor O(*z 2 ) boundary errors.

STAGGERED GR ID APPROACHES
For the 1-D problem in a homogeneous medium the ¢rst-order staggered grid approach of Vireux (1986) consists of solving the two following coupled ¢rst-order equations: where p is stress and o is velocity. Eq. (58) is approximated by the ¢nite di¡erence scheme We now consider the conventional operators for an inhomogeneous medium. It can be shown that the conventional operators approximately satisfy the general criterion for optimally accurate operators when for all grid points m, where k m is the wavenumber at z~z m and *z m is the grid interval between the mth and (m+1)th grid points. On the other hand, the Courant condition for an inhomogeneous medium is for all m, where b m is the wave velocity at z~z m . To obtain accurate synthetics, both eqs (69) and (70) must be satis¢ed. This can be achieved only if variable grid spacing is chosen so that the number of nodes per wavelength is approximately constant everywhere, and *t is chosen to be slightly less than the Courant limit. We now present numerical examples. First, we study the accuracy as a function of the time step normalized by the Courant limit. We consider a homogeneous medium (Fig. 1a) and a two-layered medium (Fig. 1b). We calculate synthetics with a duration of 750 s.  Fig. 2. (c) Density and velocity structure used in the numerical example for Fig. 3. Figure 2. Error of the waveforms obtained by the conventional (black diamonds) and modi¢ed (grey squares) operators. The time step is normalized by the Courant stability limit. (a) Errors for a homogeneous medium with a homogeneous grid. (b) Errors for a two-layered medium with a homogeneous grid. (c) errors for a two-layered medium with an inhomogeneous grid chosen so that the number of nodes per wavelength is constant throughout the medium. The dependence of the errors for the conventional operators for (a) and (c) matches the *t 2 courant {*t 2 À Á dependence predicted by eq. (68).
The source is a single force with a Ricker wavelet time history whose central frequency is 10 s. The source is at z~500 km and the receiver is at z~300 km. For the homogeneous medium we use a spatial grid with *z~1 km. For the two-layered medium we use two spatial grids: (1) a homogeneous grid with *z~1 km, and (2) an inhomogeneous grid with *z~0X5 km in the upper layer and *z~1 km in the lower layer. The latter grid has a constant number of nodes per wavelength. Fig. 2 shows the error of the synthetics obtained using the conventional (black diamonds) and modi¢ed (grey squares) operators. Fig. 2(a) shows the error (rms residual) for the homogeneous medium (Fig. 1a) with a homogeneous grid. Fig. 2(b) shows the error for the two-layered medium (Fig. 1b) with a homogeneous grid. Fig. 2(c) shows the error for the two-layered medium (Fig. 1b) with an inhomogeneous grid (chosen so that b*z~constant). The error for all of the cases considered in this paper is computed using the numerical solution for an extremely ¢ne grid as the reference solution.
As theoretically predicted, the error for the conventional operators is much worse in general than that for the modi¢ed operators, but approaches that for the modi¢ed operators if and only if *t is close to the Courant limit everywhere in the medium. Note that the errors for the conventional operators in both Figs 2(a) and (c) show the *t 2 courant {*t 2 À Á dependence predicted by eq. (68). On the other hand, the error of the solutions obtained using the modi¢ed operators is essentially uniform regardless of the choice of *t and spatial grid.
The reader might conclude that the conventional operators with a spatially varying grid interval chosen to give an essentially constant number of nodes per wavelength and a time step close to the Courant limit are su¤cient for accurate computation of synthetics, but in general this is not the case. For 2-D or 3-D problems waves propagate in various directions, and both P and S waves exist. No spatial grid for the conventional operators will yield optimally accurate solutions for such problems. On the other hand, the modi¢ed operators will yield optimally accurate solutions in general. Fig. 2(b) (uniform gridding for a heterogeneous medium) is thus probably most representative of the advantage of the modi¢ed operators for real problems. Explicit results for 2-D and 3-D cases will be presented in future papers.  We now compare the error and the CPU time for solutions obtained using the modi¢ed and conventional operators for the heterogeneous medium shown in Fig. 1(c). The results are shown in Table 1. The length of the medium is 1000 km as shown in Fig. 1(c), and the length of the time series is 500 s for each case. We use homogeneous temporal and spatial gridding. Thus, for example for the case of 500 nodes and 5000 time steps, *z~2 km and *t~0X1 s. For each case the ratio of the spatial and temporal grid intervals is *ta*z~0X05 s km {1 , whereas the Courant condition is roughly *ta*z¦0X01 s km {1 . Thus *ta*t courant &0X5.
As shown in the Appendix, the modi¢ed operators require twice the number of £oating point multiplications and three times the number of £oating point additions. As multiplication operations require more CPU time than addition operations, the theoretically predicted CPU time for the new method is about twice that of the conventional method. The actual CPU times (Table 1) follow this prediction. On the other hand the new method yields waveforms that are about 60^100 times more accurate than the conventional method. Fig. 3 shows a comparison of the accuracy of the waveforms for the case of 500 nodes and 5000 time steps. As expected based on the results in Section 2 of GT95, the modi¢ed operators yield accurate phase velocities, whereas the conventional operators cause large frequency-dependent phase errors. We can see that the modi¢ed operators are especially e¡ective for later phases for this reason. The improvement factor is 69 for this time-series, but would be larger for longer time-series. Fig. 3 raises an important general point. A visual comparison of the synthetic for the conventional operator (a) and the reference solution (c) might lead to the erroneous conclusion that there was`good agreement', because all of the expected arrivals are present at about the right arrival times. But, as shown by the residual (the second trace of a), the rms residual is actually 22 per cent. The large rms error is due to the error in phase in the synthetics for the conventional operators. The human eye is unfortunately not well suited to evaluating such phase errors. Thus a quantitative evaluation of the rms error is an essential step in the evaluation of any method for computing synthetic seismograms.
In summary, the modi¢ed operators require about twice the CPU time, but reduce the error by a factor of about 60^100 for reasonable choices of *t and *z. This means that by using the modi¢ed operators either (1) the required CPU time can be reduced by a factor of about 40 to obtain waveforms with the same accuracy, or (2) the error can be reduced by a factor of about 40 if the CPU time is kept constant.
Finally, we present an eigenvalue calculation. We numerically solve the eigenvalue problems de¢ned by eqs (32) and (45) for the inhomogeneous structure in Fig. 1(c). The computed values of as de¢ned in eq. (54) are shown in Table 2. We thus con¢rm that the stability condition given by eq. (54) is correct for both the modi¢ed and conventional operators. As a further check, we conducted calculations for values of *t slightly greater and slightly less than *z b min (1z). We con¢rmed that stable solutions were obtained for a heterogeneous medium for *z b min ¦*t¦ *z b min (1z), but that exponential instability occurred when *t b *z b min (1z).
Numerical experiments (not presented here) for optimally accurate schemes for other problems (e.g. 2-D P^SV) have shown that j j is small, but that its sign is sometimes negative.

DISCUSSION
The error of FD or other numerical schemes has generally been evaluated by considering the numerical dispersion (error of phase velocities) as a function of the number of grid points per wavelength (e.g. Alford, Kelly & Boore 1974). It is frequently assumed that the errors due to spatial and temporal discretization can be considered separately (see e.g. Fig. 4 of Fornberg 1987), but this is not the case. The key point of this paper is that the net error of the synthetics due to the combined e¡ects of temporal and spatial discretization must be considered as a single quantity. This error can be minimized by tuning the operators so that the errors due to spatial and temporal discretization come as close as possible to cancelling each other.
Second-order (in space and time) FD schemes have been deemed inferior both to FD schemes which are fourth order in space and second order in time and to schemes which are pseudo-spectral in space and second order in time (e.g. Fornberg 1987). This evaluation seems correct for conventional second-order FD schemes. However, modi¢ed second-order FD schemes of the type presented in this paper may be preferable for many applications.

ACKNOWLEDGMENTS
We thank Heiner Igel, Tatsuhiko Hara and two anonymous referees for helpful comments that greatly contributed to improving the manuscript. This research was partly supported by a grant from the Japanese Ministry of Education, Science and Culture (No. 08640525), and by the ISM Cooperative Research Program (ISM97 . CRP-A-59). NT was supported by a JSPS Fellowship for Young Scientists.