The above sections provide an essentially complete description of a numerical algorithm for the modal problem. The techniques described involve the simplest finite difference schemes and therefore lead to a relatively simple algorithm. However, as it stands, the method would be inefficient compared to certain methods based on higher-order difference schemes, e.g. Numerov's method. The simplicity of the basic method however can be retained while gaining great improvements in efficiency by using Richardson extrapolation. This is the technique used in KRAKEN and is described next.
It can be shown that the numerically derived eigenvalues vary as a function of the mesh width h, as follows:
where denotes the exact eigenvalue of the continuous problem. It is, of course, which is sought, however it is computationally expensive to evaluate with small values of h. Instead, we solve the discretized problem for a series of meshes and fit a polynomial in through the mesh points. The value of the polynomial evaluated at h=0 provides the Richardson extrapolation of the eigenvalue.
The extrapolation can be done with little effort. We denote the Richardson extrapolation using meshes by . The Richardson extrapolation for the first mesh is trivial since the polynomial fit degenerates to a constant. Thus, is identically equal to the value obtained at the first mesh. Subsequent extrapolations are obtained recursively as follows:
where h denotes the mesh width for which the extrapolation is desired.
As discussed in Chap. 2, a water/sediment interface can be handled by setting up an independent mesh in each medium and applying matching conditions at the interface. In such cases there are two mesh spacings and . As long as the ratio is kept constant when the mesh is refined one can pick either or in Eq. () and obtain the same result.
The extrapolation is implemented in an adaptive fashion. An error tolerance in the eigenvalues is computed based on the maximum range for which the field will be calculated. Then, the problem is solved for successively refined meshes and at each step the extrapolation to zero mesh width (h=0) is performed using the above recursion. When the difference between two successive extrapolations is less than , the extrapolated eigenvalue is accepted and the process is terminated. For simplicity, this convergence test is performed on a ``key-value'' which is a particular eigenvalue chosen roughly in the middle of the spectrum.
The solution of the discretized problem for each mesh can be performed as described above, however, for subsequent meshes the circumstances are changed in that we have a good estimate of the eigenvalue and eigenfunction. This information is exploited by using Richardson extrapolation in a second fashion to provide an estimate of the eigenvalue for each new mesh. In order to be confident that we extrapolations are sufficiently accurate for the root-finder this process is adopted only for the third and subsequent meshes.
A priori it is not at all obvious how well the extrapolation process should work. It depends on how rapidly the coefficients in the Taylor series for decay. Consequently, it is necessary to examine the merits of the method for individual classes of problems. Our initial experience with smooth single layer ocean acoustic problems was extremely favorable. For more complicated multilayer problems one must take care that the mesh does not straddle an interface. For this reason the difference equations have been set up to allow an independent mesh within each medium.
An additional issue arises when the sound speed is provided in a tabular form as is typical for ocean acoustic problems. If piecewise-linear interpolation of the sound speed profile is used, then the mesh must be set up to land precisely on those depths where the sound speed is tabulated. Since this must be done at each mesh, the mesh refinement is refined by simple halving. Alternatively, a new ``medium'' can be introduced for each piecewise linear layer. Finally, one can bypass both these alternatives by simply using a very smooth fit to the sound speed profile. In the present version, one has the option of using a spline fit, which however can introduce its own artifacts.
The reader who is having difficulty remembering the constraints and the various alternatives may be comforted to know that in practice the extrapolation works quite well even when the theoretical requirements are not satisfied. At one extreme, if you have a single medium with numerous discontinuities within the medium, you would probably find that the extrapolation provided no improvement over the answers obtained at the finest mesh. On the other hand, for extremely smooth profiles the extrapolation yields typically 3 additional significant digits per extrapolation.
We should also mention for the reader that contemplates modifying the code, that we have found that seemingly innocuous changes can introduce terms of odd power in h to the Maclaurin series for . The result then is to eliminate the usefulness of polynomial extrapolation in . Such terms can easily be introduced at interfaces if care is not taken that the interface condition is and that higher-order terms of odd power are not present.