For a purely acoustic problem a bisection approach is generally the most attractive solution. For this fairly broad class of problems the method will calculate a particular subset of the modes in a predictable amount of time without fail. However, for more complicated problems, for instance problems involving elasticity or complex wave speeds this technique is not applicable.
In brief, the bisection algorithm relies on the interlace property of the Sturm-Liouville problem that the number of zeroes in a trial eigenfunction increases monotonically as k decreases. (Recall that the mth mode has m zeroes.) A discrete analogue of a trial eigenfunction is the vector of principal minors of A. That is,
For any fixed the number of sign changes in the sequence indicates the number of eigenvalues greater than that particular . In addition, the function is the characteristic equation, the roots of which are the eigenvalues. The first property is used to construct isolating intervals for the eigenvalues. Thus, for each eigenvalue one computes interval endpoints which contain exactly one eigenvalue. The second property is used by a root finder to refine the eigenvalues within their isolating intervals.
In more detail this process proceeds as follows. First we compute an upper bound on the eigenvalues using Gerschgorin's theorem. (Incidentally, Gerschgorin's theorem applied to the discrete problem yields the bound where is the lowest sound speed in the problem.) This provides an upper bound for the mode search. There are an infinite number of modes so that a lower bound must be selected in some fashion. This bound is user specified, but if it exceeds the halfspace velocity in the problem, the bound is reduced to eliminate leaky modes from the problem.
Next, we take the midpoint of the interval and compute the number of modes to the right of the midpoint. Based on the number of zero-crossings in the trial eigenfunction one may decide whether the first eigenvalue lies to the left or to the right of the midpoint. The midpoint then becomes either a new lower bound or a new upper bound for the eigenvalue. This process of interval halving is repeated until the interval contains precisely one eigenvalue. With the isolating interval computed for the first mode, one then performs the same process for the second mode and so on. For subsequent modes an upper bound is available from the lower bound of the previous mode. In addition, information generated during the bisection for the first mode provides useful bounds for higher order modes. As a result the generation of all M of the isolating intervals typically requires little more than M bisection steps.
In the second stage these bracketing intervals are passed to sophisticated root finder (Brent's method ) which combines bisection, linear interpolation and inverse quadratic interpolation to provide an estimate of the eigenvalue. This latter root finder is guaranteed to converge given an isolating interval for the root. The whole process is very efficient (compared to brute force linear searching) and robust (compared to techniques which rely on asymptotic estimates to provide initial guesses for a root finder).
In practice, the sequence is replaced by where,
The sequence then satisfies the following recursion:
The number of sign changes in is the same as that for the original sequence, since has no sign changes. However, is well scaled and requires half as many floating-point operations to compute. Incidentally, this latter sequence is precisely equivalent to what one obtains in a simple one-sided shooting scheme. Instabilities often occur in this sequence which are removed by simply scaling it down as it gets too large. The scale factors are retained and passed separately to the root finder. Surprisingly, this instability does not in the least affect the accuracy of the eigenvalues. This is proven in Wilkinson's text .
We have glossed over a difficult aspect of this method. Namely, the information provided by the Sturm sequence (regarding the number of modes to the right of some trial eigenvalue) is only justified for purely acoustic problems with free or rigid boundaries. However, for homogeneous half-spaces or elastic half-spaces there exists a simple correction to the count[16,71]. For elastic or acousto-elastic problems with more than one elastic layer or elastic gradients there is no useful extension and we are forced to use the deflation method described next.