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.
