In this algorithmic and numerical analysis of Rayleigh-wave dispersion computations on a spherical, gravitating earth, we first simplify and optimize the existing algorithms based on the direct, Alterman-Jarosch-Pekeris (AJP) formulation; compare computation speed with the fastest existing methods; investigate the relation between integration step size when using direct-integration methods, and layer thickness when using homogeneous-layer methods; and give specific details concerning numerical limitations. Even in highly optimized form, the AJP, direct-integration formulation is slow relative to the best of the techniques based on the flat-, homogeneous-layer approximation. This motivates the development of an improved computational algorithm given as the final result of this report. Loss-of-precision problems appear to be intrinsic to the AJP formulation. At a fixed period, this results in the attainable accuracy of the phase velocity decreasing as radial mode number increases; for fixed accuracy in the phase velocity, as period decreases the maximum mode number that can be treated successfully decreases. Specific numerical relationships among period, mode number, and attainable accuracy are presented.
For gravitating structures, a solution to this loss-of-precision problem is given which is based on Gram-Schmidt orthogonalization. Our analytical and numerical studies show that: the dispersion function—roots of which determine dispersion properties—is invariant under orthogonalization, which ensures a smooth, regular variation of this function as roots are sought; the required modifications that orthogonalization introduces into the computation of actual components of motion are quite straightforward (full details are given); and, the addition of orthogonalization increases computation cost by only 5 to 10 per cent. Down to periods of 10 sec, for all radial mode numbers up to 90 to 100, the numerical tests indicate that orthogonalization brings the loss-of-precision problem under complete control. All analyses in our investigation show that this problem exists, with the same properties, when integrating the equations of motion in either upward or downward directions; there is no preferred direction of integration relative to this difficulty. Likewise, there is no preferred direction relative to accuracy or overflow problems. In fact, the overflow features are shown to be precisely the same, in form and magnitude, for the two directions. The advantage that we do find for downward integration is the degree to which it simplifies computational algorithms.
The problem of Rayleigh-wave propagation on a gravitating, spherical earth is characterized by a linear system of six differential equations in six unknowns; for a nongravitating earth this system reduces to four equations in four unknowns. We have succeeded in reducing the sixth-order system to fourth order, while retaining the effect of gravity. The speed of computation for a gravitating earth is thereby increased to about that of the nongravitating case, i.e., to that for nongravitating spheres and to the speed of the fastest of existing algorithms for flat, nongravitating structures treated with the homogeneous-layer approximation. Even greater speeds are possible with the new algorithm if delta matrices are used to control precision loss.