As illustrated in Fig. 3.1 we divide the
interval **0<z<D** into **N** equal intervals to construct a mesh of
equally spaced points

where **h** is
the mesh width given by . Furthermore, we shall use the
notation . The number **N** should be chosen large
enough so that the modes are adequately sampled; usually **10** points
per wavelength are sufficient.

**Figure:** Finite-difference mesh.

We shall assume for the moment that the density is constant, yielding the modal problem

where the primes denote differentiation with respect to **z**.
Following a standard procedure for deriving finite
difference equations, we use the Taylor series expansion to obtain

Rearranging terms, we obtain a * forward-difference* approximation
for the first derivative,

An approximation to the first derivative is therefore,

An improved approximation is obtained by using the governing equation Eq. () to evaluate the first term in the forward difference approximation. That is, we substitute

This yields the approximation

Similarly, a * backward-difference* approximation is obtained starting
with the Taylor series

yielding, the approximation

and the approximation

Finally, adding Eqs. () and () we
obtain a * centered-difference* approximation to the second
derivative:

With these finite-difference approximations in hand, we can proceed to replace the derivatives in the continuous problem with discrete analogues. Let us recall the continuous problem:

Using the centered, forward and backward difference approximations for the ODE, the top and bottom boundary conditions we obtain:

We next write the first of these equations as

Then collecting the difference equations together we obtain an algebraic eigenvalue problem of the form:

Here, is the vector with components . These components are approximations of the
eigenfunctions of Eq. () evaluated at the mesh points. In
addition, **A** is a symmetric tridiagonal matrix defined
by

where the coefficients and are defined by

and,

We have consciously introduced a scaling factor of in every row here. For a constant-density problem with a single mesh width this would serve no purpose; however, later we shall consider multiple layers for which this scaling is the natural one.

Note that for a pressure-release surface the ratio which
appears in the boundary condition goes to infinity. In this case,
vanishes and we can simply delete the first line and column
from the matrix eigenvalue problem. Furthermore, if the functions
are independent of **k** (as happens for the
pressure release surface and rigid bottom conditions) then the above
problem is a standard algebraic eigenvalue problem and can be solved
using standard routines. In general, only the lower-order modes will
be sufficiently accurate: the higher order modes are undersampled by
the finite-difference mesh. Thus, routines which are designed to
extract a subset of the eigenvectors and eigenvalues are desired.

The problem is a non-standard eigenvalue problem because the eigenvalue enters in a functional form through bottom boundary conditions. (For perfectly rigid or free boundary conditions the problem reduces to a classical algebraic eigenvalue problem.) We shall discuss the solution of the eigenvalue problem in Sect. 3.3.

Tue Oct 28 13:27:38 PST 1997