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[13]. For more complicated
multilayer problems one must take care that the mesh does not
straddle an interface[64]. 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.
