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
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.