Path-Following Methods for Calculating Linear Surface Wave Dispersion Relations on Vertical Shear Flows

The path-following scheme in Loisel and Maxwell (SIAM J Matrix Anal Appl 39(4):1726–1749, 2018) is adapted to efficiently calculate the dispersion relation curve for linear surface waves on an arbitrary vertical shear current. This is equivalent to solving the Rayleigh stability equation with linearized free-surface boundary condition for each sought point on the curve. Taking advantage of the analyticity of the dispersion relation, a path-following or continuation approach is adopted. The problem is discretized using a collocation scheme, parametrized along either a radial or angular path in the wave vector plane, and differentiated to yield a system of ODEs. After an initial eigenproblem solve using QZ decomposition, numerical integration proceeds along the curve using linear solves as the Runge–Kutta F(·)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F(\cdot )$$\end{document} function; thus, many QZ decompositions on a size 2N companion matrix are exchanged for one QZ decomposition and a small number of linear solves on a size N matrix. A piecewise interpolant provides dense output. The integration represents a nominal setup cost whereafter very many points can be computed at negligible cost whilst preserving high accuracy. Furthermore, a two-dimensional interpolant suitable for scattered data query points in the wave vector plane is described. Finally, a comparison is made with existing numerical methods for this problem, revealing that the path-following scheme is the most competitive algorithm for this problem whenever calculating more than circa 1,000 data points or relative normwise accuracy better than 10-4\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-4}$$\end{document} is sought.


Introduction
We develop a path-following method to numerically calculate the dispersion relation curve for linear surface waves atop a horizontal current of arbitrary depth dependence, an adaptation of the scheme developed by Loisel and Maxwell [29].

Background and Motivation
We consider a classical problem from the study of wave-current interactions [38][39][40], [36,sec. IV], [12,34,35,43], that of linear surface waves upon a vertical shear current. In other words, waves are considered as perturbations to first-order in wave steepness = ka upon some zeroth-order depth-dependent horizontal background flow that, in general, has non-constant vorticity. The problem geometry is shown in Fig. 1. Waves in this regime are dispersive with their behaviour being entirely characterized by the dispersion relation, i.e. the dependence of phase velocity c(k) on the wave vector k whose modulus k is the wavenumber.
When a vertically varying shear current is present beneath a water surface, the dispersion of water waves riding atop it is affected in a complicated way. The well-known expression for the phase velocity of a linear surface wave in inviscid, initially quiescent water of depth h, c 0 (k) = c 0 (k) = √ (g/k) tanh(kh), where g is gravitational acceleration, does not hold when considering general horizontal depth-dependent currents of the form U(z) where z is the depth. In the general case, c(k) is an eigenvalue of the boundary value problem formed by the Rayleigh stability equation (inviscid Orr-Sommerfeld equation) with a bottom Dirichlet boundary condition and a freesurface boundary condition. Only in the cases of quiescent water, uniform current, and strictly linearly varying current is a closed-form expression for c(k) known for all permissible values of k (and vice versa). Figure 2a In a theory setting, the typical situation is that c(k) is desired for a given k, such as when initial value problems are solved via a Fourier transformation in horizontal spatial directions. In a laboratory setting, however, wave-makers operate at a given frequency ω = kc, with k as the dependent parameter. Likewise, when a monochromatic wave propagates in an environment of varying current and bathymetry, its frequency is the conserved quantity whereas the wave vector k changes according to local conditions. The former setting of determining c(k) from k, we refer to as the 'forward problem'; the latter, finding k(c), we refer to as the 'backward problem' (distinct from the 'inverse problem' which aims to infer the current profile when c(k) is known). Our primary interest in this paper is in the forward problem but, for purposes of completeness, we also describe how to solve for k(c) using a collocation scheme.

Existing Methods to Calculate the Dispersion Relation
Unsurprisingly, a number of approaches exist for determining c(k). The oldest approach was to approximate a depth-varying unidirectional velocity profile U (z), where z is the vertical coordinate, by a piecewise linear profile [13,38,39,50]. U (z) could realistically be divided only into a small number of pieces before the formalism would become impractical for manual calculation. The availability of modern computing permits extending the piecewise linear approximation to arbitrarily many pieces [45,57]. Despite the procedure being physically transparent, it converges only as N −2 [57], where N is the number of pieces.
Another popular approach is integral approximation methods. An expression was derived by Stewart and Joy [47] (infinite depth) and then generalized by Skop [44] (finite depth) and further extended by Kirby and Chen [26] (finite depth to second order); these approaches are usually valid to within a few percent for typical oceanographic current profiles. An alternative able to tackle cases with extreme shear was presented recently by Ellingsen and Li [21].
Numerical techniques are applicable to the problem, for example, a shooting method was used by Dong and Kirby [14]. The problem is also amenable to solution with a naive shooting method using a Runge-Kutta numerical integrator and Newton's method. A recent addition is the so-called Direct Integration Method (DIM) due to Li and Ellingsen [27] which is versatile and arguably the best available method when not so many points are sought on the dispersion relation curve and only moderate accuracy is required. There are various ways to apply spectral methods; example usage of a collocation method is described later in Sect. 3.
It is also possible to calculate a high-order interpolant of the curve using the collocation method for pointwise sampling; this is used as a comparison in our performance tests in Sect. 7 using Chebfun [17] to construct the interpolant.
The path-following method presented in Sect. 4 calculates not only a piecewise Hermite interpolant but also exploits the local analyticity by integrating along the curve. This confers a significant efficiency advantage in most cases.

Methodology and Approach
Engineering and the natural sciences are replete with eigenproblems for ordinary differential operators which depend on a finite set of parameters. We are interested in problems which are parametrized by a single real variable.
The canonical solution approach involves conversion to an algebraic problem via spatial discretization, which often leads to polynomial or even nonlinear eigenproblems of potentially large dimension. These can be solved using classical techniques for each sought parameter value. This strategy may become prohibitively expensive when the computation must be repeated many times. It can also be difficult to take advantage of the nearness of solutions for small parameter variations, forcing full calculations for each point in the parameter space. An alternative approach is to solve the discretized eigenproblem once for a fixed parameter value then use the local piecewise analyticity of the eigenvalue and eigenvector [4,24,37] to calculate along the solution curve using a path-following or numerical continuation algorithm.
In a more general setting, this comprises a numerical continuation method whereby the parameter-dependent solution is calculated as an implicitly defined curve [3]. Homotopy methods have a similar philosophy but introduce an artificial parameter to parametrize a convex homotopy to map from the solution of an 'easy' problem to the solution of the actual problem [28]. These methods tend to use predictor-corrector schemes such as pseudo-arclength continuation or similar approaches. We make reference to the homotopy method in [30] and the invariant subspace methods in [8], [10] as relevant examples. For a recent approach that shares a strong philosophical similarity with the material herein for working with time-varying matrix eigenproblems, albeit using different techniques (look-ahead finite difference formulas), see [54] and [53]. This paper is concerned with repurposing the path-following technique used in [29] to solve for the dispersion relation associated with the geometry described in Sect. 1.1. This problem is particularly suited as a motivating example of the technique: it is conceptually simple, it has an eigenvalue-dependent boundary condition, it is well known from both the waves literature and hydrodynamic stability, and there is a practical requirement for efficient numerical solution.
We summarize our approach as follows. The original eigenproblem is spatially discretized using a collocation method, implicitly incorporating the boundary conditions, to obtain a parameter-varying system of equations that are then differentiated to yield an under-determined system of ODEs. An additional constraint is then included. After performing an initial eigenproblem solve, numerical integration can proceed along the solution curve using linear solves as the Runge-Kutta F(·) function. Specifically, the Dormand-Prince RK5(4)7M method [15, p. 23], which incorporates adaptive step-size control, is used. A piecewise polynomial interpolant provides dense output.
This approach is in contrast to standard predictor-corrector style numerical continuation, see, for example, the description in [3]. Rather than applying a correction at each step, reliance is placed entirely on the adaptive step-size control in the Dormand-Prince 5th-order scheme to ensure that the error is maintained within the desired tolerance. As demonstrated later in Sects. 5 and 7, this approach produces an algorithm that efficiently computes the solution with high accuracy.

Outline of the Paper
We begin by formally describing the geometry of the physical problem and some problem-specific background in Sect. 2. The collocation scheme used is briefly described in Sect. 3. The path-following method is described in Sect. 4 along radial and angular one-dimensional paths in the k wave vector plane, and also for general scattered data points in the two-dimensional k plane. In Sect. 5, we provide numerical results to determine the expected accuracy of the collocation and path-following methods. In Sect. 6, we describe an adaptive partition of unity method to address the problem of the eigenvector becoming numerically singular at the surface with increasing k. In Sect. 7, we evaluate the relative performance characteristics of the various methods and in Sect. 8, we describe how to choose optimal parameters. Finally, in Sect. 9, we provide some conclusions.

Problem Description
Wave-current interaction problems are often studied by adopting a modal linear perturbative approach: waves are considered as first-order perturbations of a stationary, incompressible, and inviscid bulk fluid flow. In this context, waves are dispersive with the phase velocity of a wave dependent on the wave vector in a nonlinear manner. The relationship between the wave vector, k, and the phase velocity, c, is termed the dispersion relation and is determined by factors such as water depth and background current.
In our context of first-order free-surface waves atop a vertical shear flow, the problem reduces to finding solutions of the eigenvalue problem formed from the Rayleigh stability equation and appropriate boundary conditions. The Rayleigh equation is a second-order ODE that is equivalent to the Orr-Sommerfeld equation when viscous effects are neglected. A solution to the eigenproblem will yield a {k, c, w(z)} triplet where w(z) is the associated eigenfunction. The eigenfunction can be used to reconstruct the velocity and pressure field for the corresponding wave vector [27, sec. 4.4]. Notably, there is substantial overlap with the literature on linear stability theory, e.g. [41, ch. 2] or [16], albeit with different boundary conditions.
Closed-form expressions for the dispersion relation for free-surface waves exist only in specific scenarios such as quiescent water [33], [ [19,20].
To make matters more concrete, we adopt the general approach used in [27,43] and refer the reader to the derivations therein for full detail. For expediency, the approach is only summarized here. The physical model is depicted in Fig. 1a. Dimensional quantities are denoted with an acute accent, e.g.ḱ.
The background flow is specified by a velocity profileÚ(ź) = (Ú x (ź),Ú y (ź)): a two-dimensional vector field describing the bulk fluid velocity in the horizontal plane for a givenź ∈ [−h, 0] whereh is the constant fluid depth and the unperturbed surface is atź = 0. Let We useh as a characteristic length scale andÚ m as a characteristic velocity, to arrive at the following nondimensionalization: so that a notional shear strength can be expressed with Froude number, F = U m / gh. The velocity profile must be suitably regular, so we impose that U x , U y ∈ C 2 ([−1, 0], R). We also require c parallel to k and we henceforth only refer to scalar c = |c|. It is assumed that the current can influence the waves but not conversely and, for clarity of exposition, we neglect surface tension. We adopt the Ansatz that perturbations are plane waves proportional to exp[i(k · x − ωt)], where the (scalar) wave angular frequency is ω = kc with k = |k|, and use a Fourier transform in the horizontal plane (coordinate space quantities indicated with a tilde), The velocity perturbations along the x, y, and z axes are, respectively, u = u(k, z), v = v(k, z), and w = w(k, z), whilst the pressure is p = p(k, z). The governing equations are the linearized Euler equations and incompressibility condition, e.g. recall [43]: After rearranging, the Rayleigh equation (4a) is obtained which we write along with the relevant free-surface condition Eq. (4b) (combined kinematic and dynamic boundary condition) and Dirichlet boundary condition for the fluid bottom Eq. (4c), for μ = k 2 . The velocity and pressure field for a specific wave vector can be recovered by substituting the appropriate eigenpair into [27, sec. 4.4], Often, the physical problem is posed so that a scalar U (z) is aligned along the x-axis as U(z) = (U (z), 0), with scalar k and c, see Fig. 1b. However, there is no difficulty in solving with a general U(z) by simply projecting U(z) along k as described in Sect. 2.4. Full three-dimensional considerations only come to the fore when calculating the velocity and pressure field; the eigenvalue problem described by Eq. (4) is always one-dimensional in the sense that the eigenfunctions are a function of one variable, z. This is all straightforward but it is worth stating explicitly to avoid the inherent ambiguities and possible confusion caused by '1d' / '2d' or '2d' / '3d' descriptions, e.g. [21, sec. 1.1]. In particular, we refer to the problem with scalar U (z) as the reduced problem and with vector U(z) as the general problem. In Sect. 2.4, we note that the reduced problem is equivalent to solving the general problem in a radial 'slice' at some fixed angle θ 0 .
So far, we have deliberately avoided specifying which variable is the sought eigenvalue in Eq. (4): it can be chosen as either μ or c, with its counterpart parametrizing the problem and always chosen to be real valued, in a similar manner to [11, sec. 7.4]. Since we are always choosing the parameter to be real valued, we are concerned with a subset of the parametrized spectrum in each case and can plot this arrangement as a function of the parameter. • The spectrum for μ(c), for c in some suitable interval, is comprised of a countably infinite set of eigenvalues. The dominant eigenvalue, μ 1 = k 2 , in this case is the only positive eigenvalue, and corresponds to a propagating wave (for ±k). The negative eigenvalues correspond to the countably infinite set of discrete k arranged along the imaginary axis and are not mentioned further. See Fig. 3a. • The spectrum for c(k), for k > 0, has both discrete and essential parts (c such that U (z) − c = 0, regular singular points of the ODE). In particular, there are two discrete eigenvalues usually labelled c + (k) and c − (k). Further discussion is available in [16, p. 137]; in particular, our interest shall be in computing the dominant 'regular inviscid solution', the smooth c + (k) curve. For some k, c + (k) may be located within the same interval as the essential spectrum and, therefore, can be difficult to identify within numerical solution sets, see Fig. 3b.
The c ± exhibit odd symmetry with respect to k: −c + (k) = c − (−k); the associated eigenfunctions are also odd symmetric, −w + (k, z) = w − (−k, z). Furthermore, it should be noted that the one-sided limits lim k→0 + c(k) and lim k→0 − c(k) exist but depend on angle θ . So c ± (k) is, in general, undefined at the origin.

Problem Types: Forward, Backward, and Inverse
We distinguish three types of problem.

For velocity profile, U (z), and collection of wavenumbers, {k
We denote this as the forward problem. 2. For velocity profile, U (z), and collection of velocities, {c ( j) } J j=1 , calculate associated wavenumbers {k(c ( j) )} J j=1 . We denote this as the backward problem. 3. For collection of wavenumber and phase velocity pairs, {(k ( j) , c ( j) )} J j=1 , determine velocity profile, U (z). We denote this as the inverse problem, which is of an entirely different nature and herein not considered further.
Both the forward and backward problems usually amount to calculating sufficiently many {k ( j) , c ( j) } pairs as to adequately specify the full dispersion relation for a given velocity profile. For practical purposes, these problems are almost always posed as Eq. (3) with exponent of the form i(kx − ω(k)t) (see [36, p. 77, eqn. 4.6], [26, eqn. 1], [43, eqn. 2.4], etc). This, by default, presupposes solving the forward problem. Furthermore, handling of critical layers as in Sect. 3.4 is, from a numerical standpoint, easier for the forward problem. Therefore, most of this paper concerns solution of the forward problem. Unless otherwise specified, we shall compute c + since c − follows automatically from symmetry.
There are a few exceptions when solving the backward problem is more suitable such as for wave problems with periodic or stationary time dependence, e.g. in ship waves. Hence, for purposes of completeness, we also describe solving the backward problem using a rudimentary collocation scheme and the basic properties of the spectrum.

Summary of Numerical Schemes
For ease of reference, we denote the various numerical schemes used or described in this paper: • CL-c Collocation scheme for the forward problem, see Sect. 3.3.
• CL-k Collocation scheme for the backward problem, see Sect. 3.2.
• PF-R-r-c A path-following scheme with dense output to solve the forward reduced problem along a fixed angle θ 0 in the k-plane, see Sect. 4.2.1. • PF-R-a-c A path-following scheme with dense output to solve the forward reduced problem for c k 0 (θ ) along a fixed circle of radius k = k 0 with varying θ in the kplane, which we term the angular solution, see Sect. 4.2.2. The CL-c and CL-k schemes incur an eigenvalue calculation for each query point, so the computational cost will increase linearly with the number of query points. The arrangement of points in the k-plane for the CL-schemes can be random (scattered data) without affecting computational cost.
The PF-R-r-c path-following algorithm is two staged: it first performs numerical integration to calculate control points along a radial 'slice' at fixed angle θ 0 , which incurs a nominal initial computational cost; query points on that slice are then calculated using a Hermite interpolant. Although the computational cost of interpolation is linear in the number of query points, it is so lightweight as to be of almost negligible cost in most situations: so after the initial computation, very many query points can be calculated efficiently. The angular PF-R-a-c scheme is similar but instead calculates along a circular path at a fixed radius k 0 .
The PF-G-c scheme is more involved because we accept query points in the kplane with no assumption on arrangement, i.e. scattered data. A naive approach would incur a complete first stage calculation of PF-R-r-c for every query point, which is unacceptable. The PF-G-c scheme instead precalculates a two-dimensional polar grid of suitable control points and then can interpolate for query points at negligible cost.
Note that all methods presented can also make available the eigenfunction w(z) so that the velocity and pressure field can be reconstructed using Eq. (5).

General to Reduced Problem
The general problem Eq. (4) can be simplified to a reduced problem by projecting U along k, see, for example, [36, p. 77], [43, p. 566]. Define the scalar velocity profile for the reduced problem as where θ is taken to be the standard angular polar coordinate for k.

Velocity Profiles and Parameters Used
For later numerical tests, we define a test velocity profile U T for the reduced problem, as shown in Fig. 1b, We choose the physical depthh = 20 and shear Froude number as F 2 = 0.05. We choose nondimensional k ∈ I k := [ 0.025, 250]. This corresponds broadly to gravity waves in the air-water interface regime [33, p. 4] with shortest period ≈ 0.2s. The function is chosen as a suitable test candidate because it has several stationary points and cannot be approximated exactly over a finite-dimensional polynomial basis. For the figures produced from PF-G-c shown later in Fig. 7, we use U T (z) along the x axis and an approximation of a flow in the Columbia River on the y axis, which we denote U CR . U CR is defined by a sixth-degree polynomial as in [32]. For the plots in Sect. 4.3, U CR is scaled to have F 2 ≈ 0.01; it is also used in Sect. 3.4 with F 2 = 0.05.
For use later, define

Discretization of the Equations
We use sans serif notation to indicate matrices (uppercase) and vectors (lowercase), e.g. U or w, to distinguish from their continuous counterparts. Let ζ j = cos(( j − 1)π/N z ), j = 1 . . . N z + 1 be the Chebyshev-Gauss-Lobatto collocation points on [−1, 1] (second-kind points). We use the change of variable z = (h/2)(ζ − 1) to map ζ j to z j ∈ [−1, 0] in nondimensional coordinates and let z := [ z 1 , . . . , z N z +1 ] T be the associated column vector. Let D be a corresponding square differentiation matrix (in practice, we calculate the first-and second-order differentiation matrices, D and D 2 , using poldif.m from the Weideman-Reddy library [55] and then apply the 'negative sum trick' as detailed in [6]). We define vector discretizations of the velocity profile, Quantities u and u are similarly defined. Let w := w(z). Define the diagonal matrices U := diag(u), U := diag(u ), and U := diag(u ).
The problem is a two-point boundary value problem so is amenable to the standard 'row-replacement' strategy, see, for example, [52, ch. 7]. Specifically, we aim to construct eigenvalue equations which discretize the governing equation Eq. (4a) using the 'interior' rows 2 through N z of the differentiation and coefficient matrices. The free-surface boundary condition Eq. (4b) is incorporated in the first row of the matrices. The bottom Dirichlet Eq. (4c) boundary condition is accounted for by eliminating the last row and column in the matrices, as in [52, ch. 7]. For notational convenience, we define 'interior' differentiation and velocity profile matrices as (in other words, eliminating the first and last rows, and last column). We also define a free-surface differentiation vector as the first row of D, d f := D 1m , m = 1 . . . N z , again without the last column.

Backward Reduced Problem (CL-k)
Treating c as a parameter and k as the eigenvalue then using the projection from Sect. 2.4, we obtain a regular Sturm-Liouville problem on z ∈ [−1, 0] with eigenvalue μ = k 2 , The discretization of Eq. (7a) proceeds in an obvious manner, Discretizing Eq. (7b) into a row vector gives, Write to obtain the generalized eigenvalue problem, Note that the only effect of B is to ensure that the row of A with the free-surface boundary condition is set equal to zero and is not dependent on the eigenvalue. Equation (11) can be solved in several ways, e.g. using MATLAB's implementation of QZ as eig (A,B). For a given c, the solution of Eq. (11) yields the greatest N z eigenvalues from the countably infinite set of discrete μ j eigenvalues, ordered μ 1 > μ 2 > . . ., associated with the infinite-dimensional operator. The only positive eigenvalue is μ 1 . Since k = √ μ, μ 1 gives the only real k. Our interest is in propagating waves; therefore, k ± = ± √ μ 1 is the sought solution. This is shown in Fig. 3a.

Forward Reduced Problem (CL-c)
Now, treating k as a parameter and c as the eigenvalue, we rewrite the reduced problem Eq. (7) to emphasize the quadratic eigenvalue dependence in the free-surface boundary condition, We initially discretize Eq. (12a) as We proceed by expressing the free-surface condition as coefficients of the powers of c, In the same manner, we now separate Eq. (13) into powers of c, Define To obtain the sought solution, we solve the quadratic eigenproblem, There are several techniques to solve the quadratic eigenvalue problem. A linearization and then using a QZ decomposition works well in this setting. MATLAB's polyeig(A2,A1,A0) implements such a linearization although some care must be taken. In particular, the A 2 matrix is badly rank-deficient. As a consequence, the QZ algorithm will return infinite and large-but-finite eigenvalues which are merely artefacts of the numerical method and must be removed.
From the solution of Eq. (17), it is the c + eigenvalue that is sought. The result set will also include an eigenvalue for c − and other finite eigenvalues that are regular singular points of the underlying ODE. These regular singular points, or critical layer points, are sampled from the essential spectrum: c such that U (z) − c = 0 for some z ∈ [−1, 0], see [16] and Fig. 3b.
Let α j for j = 1, . . . , N j ≤ 2N z be the eigenvalues obtained from the solution of Eq. (17) and w j be the corresponding eigenvectors. The infinite and very large in magnitude α j are discarded. The solution is found by determining the j such that α j = c + can be identified (to within numerical error). There are two possible cases, either (i) c + is the greatest element of α j such that c + / ∈ [U min , U max ], or (ii) c + ∈ [U min , U max ] (a critical layer has been encountered).
The former case is trivial: merely pick the largest eigenvalue that is not in [U min , U max ]. The latter case may be handled using the methodology described in Sect. 3.4 below.
Although the collocation method converges exponentially fast with respect to increasing N z , the effects of roundoff error and ill-conditioning in the differentiation matrices soon start to dominate [11]. As shown later in Sect. 5, in double precision, the relative error in calculated c is at best around 10 −10 for this problem. Techniques such as domain decomposition may be used to reduce the error further but since these techniques are relatively well known, implementation details have been omitted.

Solving with Critical Layers in the Forward Problem
In most cases, one of the α j corresponds to c + and another to c − . The remaining α j eigenvalues are sampled from the essential spectrum and shall be described as spurious eigenvalues to match the terminology used in [11, ch. 7] and [45]. The underlying eigenvalue problem is quadratic and, as might be expected, the spurious eigenvalues occasionally manifest as complex conjugate pairs due to collisions with other spurious eigenvalues. In pathological situations, an eigenvalue collision between c + and a spurious eigenvalue will create a complex conjugate pair causing the sought c + to disappear from the solution set. This is dealt with in more detail in Sects. 3.4.1 and 3.4.2 below.
A similar approach to filtering of spurious eigenvalues by looking for singular behaviour in the eigenfunction away from the endpoints was used in [45]. However, it does not perform robustly when the order of discretization is too large or if the sought c + is too close to a spurious eigenvalue.

Eigenvector Filtering
The sought α j is usually within a good error margin of the exact c + but identifying this from the spurious eigenvalues is nontrivial. The eigenvectors corresponding to these spurious eigenvalues tend to exhibit features that permit identification aside from singular behaviour in the interior of the interval. These include 1. nonzero near z = −h; 2. nonzero imaginary component; 3. for suitably large k, the decay near the surface is not sufficiently 'fast'; 4. dependence on the order of discretization; and, 5. singular behaviour in the interior of the domain (detected using MATLAB's findpeaks()).
The approach taken is to successively eliminate eigenvectors from the set {w j } j until the only eigenvectors remaining are those corresponding to the sought α j . Computationally cheap methods are used as a first-pass filter (everything in the above list except peak-finding). If there are still spurious eigenvalues remaining then the using the peak-finding strategy will give the sought solution.
This process is demonstrated using the negative of the Columbia River profile, −U CR and N z = 80: Fig. 4a shows an example of the full set of eigenvectors for all finite eigenvalues at a fixed k; Fig. 4b shows the remaining eigenvectors after the initial filtering step; finally, Fig. 4c shows the sought eigenvector after all filtering steps are completed.

Eigenvalue Collisions
Although the spurious eigenvalues are sampled from the essential spectrum, when solving the finite-dimensional Eq. (17), this distinction is essentially lost. Eigenvalues can collide forming complex conjugate pairs. In particular, the sought c + eigenvalue may collide with a spurious eigenvalue causing c + to disappear from the result set and so the method in Sect. 3.4.1 would be ineffectual. By detecting these points and temporarily adjusting the size of the discretization so that the spurious eigenvalues are shifted, the sought eigenvalue reappears (similar in spirit to [11, sec. 7.5]). The detection is done by looking for complex conjugate pairs and matching-up against the eigenvector filtering, i.e. there should be an eigenvector corresponding to c + , if it is missing then a collision should be expected. This is demonstrated in Fig. 5. a b c Fig. 4 Eigenvector filtering for −U CR profile at k ≈ 85. a Eigenvectors for all finite eigenvalues. b Eigenvectors remaining after first-pass filtering. c Eigenvectors after filtering using peak-detection

Review of Loisel-Maxwell Path-Following Method for the Field of Values
In [29], the authors describe a path-following method to calculate the field of values boundary of a matrix, which we now briefly summarize. It concerns the solution of a parametrized Hermitian eigenvalue problem (which bounds the projection of the field of values onto the real axis), Here, and in the remainder of the paper, the overdot notation is used to indicate derivatives with respect to the problem parameter. This is to emphasize the parameter-varying or "time-varying" nature of the problems. Note that Eq. (18) is well-defined except perhaps for a finite number of τ j due to elementwise analyticity of the elements of H(e iτ A) and the analyticity, up to ordering, of the eigenvalue and eigenvectors. Differentiating Eq. (18) gives Since the system is under-determined, an additional constraint that u(τ ) must be perpendicular to its (elementwise) derivative is included: u(τ ) * u (τ ) = 0. In [29], this was derived by differentiating the algebraic normalization constraint u(τ ) * u(τ ) = 1 and fixing the 'phase-factor'. In retrospect, this appears to have been the same rationale applied in [49, p. 57] and similar to the normalization considerations in [9]. However, the normalization constraint is violated in the numerical integration and must be manually enforced at each step. In this sense, the derivation of the constraint does not offer a satisfactory explanation of its action. An anonymous referee points out that upon application of the Fredholm alternative, it is possible to express the derivative of the eigenvalue aṡ where (·, ·) is the usual inner product. Note that this expression is none other than the Hadamard first variation formula from [49, p. 57, eqn. 1.72].
Including the constraint gives the fully determined system, which can be rewritten in matrix form, The system described by Eq. (22) can be solved for [u(τ ) * λ (τ ) * ] * using a linear solver and used as the F(·) function for a Runge-Kutta numerical integrator, which then generates control points along the curve. The authors use the Dormand-Prince RK5(4)7M method [15, p. 23] and interpolation method of Shampine [42, p. 148]. The near-interpolant solution from this method is 5th-order accurate. There is comprehensive discussion presented of sharp points and flat segments on the boundary, and relevant event processing to detect these exceptional points in the path-following algorithm.

Path-Following Method for Forward Reduced Problem
We now extend the same process to the quadratic eigenvalue problem posed in Sect. 3.3. Recall Eq. (17), which upon differentiation (indicated with overdot) with respect to k gives We further impose that w(k) * ẇ (k) = 0. Writing in matrix form, This is the general form in which the structure is clear. In the subsections below, we perform the same derivation but include the specific expressions for the radial and angular paths including boundary conditions. The computational approach taken is analogous to [29]: an initial eigenpair {c 0 , w 0 } is calculated using CL-c. Then using Eq. (25) to solve for [ẇ(k) * ċ (k) * ] * , numerical integration can proceed along the curve in both directions. Hermite interpolation can then be used to query at arbitrary k.
Care must be taken to ensure that M remains nonsingular throughout the parameter ranges required. This may be done computationally at each Runge-Kutta step by, for example, checking the reciprocal condition number. Due to the bordered matrix structure, M only becomes singular under certain degenerate circumstances such as encountering nonsimple eigenvalues. For background on use of bordered matrices in numerical continuation, see, for example, [25] or [10]. From knowledge of the problem at hand, we know that the only region in which there may be nonsimple eigenvalues is for k such that c(k) ∈ [U min , U max ] (whereupon a critical layer is encountered). This is discussed further in Sect. 4.2.4. For other problems, more careful handling of points of higher algebraic multiplicity would likely be required.

System of Equations Along Radial Slice at Fixed Â (PF-R-r-c)
For PF-R-r-c, we fix angle θ = θ 0 and parametrize by k. Thus, we are in the setting of the reduced problem with the constant velocity profile being the relevant reduced velocity profile, U k (z). Writing Eq. (7a) in matrix form with c as eigenvalue and explicit dependence on parameter k, For notational succinctness, we use the shorthand c = c(k) and w = w(k). Differentiating Eq. (26) with respect to k (indicated by an overdot) gives The free-surface condition can be written as Differentiating Eq. (28) with respect to k, Upon rearranging terms, we define the block matrices: so that the system of ODEs can now be written in matrix form as In the same manner as in Sect. 3.1, the bottom boundary condition is imposed by eliminating the corresponding row and columns [52, ch. 7], only the free-surface is included explicitly.

System of Equations Along Angular Circle at Fixed k (PF-R-a-c)
For PF-R-r-c, the angle θ and hence the velocity profile was held constant. For PF-Ra-c, we instead hold k constant and seek to use a θ angular dependence. Therefore, we must also specify the parametrization of the velocity profile.
So that in matrix form and, upon differentiation with respect to θ (indicated with an overdot), Our starting point is the same, we use Eq. (26) but instead hold k constant and take the derivative with respect to θ . Temporarily adopting the abbreviated notation U int = U int (θ ), w = w(θ ), and c = c(θ ): As before, the free-surface condition is Eq. (28), which we take the derivative of with respect to θ using the shorthandu 1 =u 1 (θ ), In a similar manner as before, we define the block matrices: so that the system of ODEs can now be written in matrix form as . (35) Note that the P and Q matrix have the same structure as in Eq. (31), it is R that changes.

Path-Following Algorithm Specification for Reduced Problem
We describe the algorithm for PF-R-r-c, the algorithm for PF-R-a-c follows in the obvious manner. Using the definitions of P, Q, R from Eq. (30a,30b,30c) define matrix and vector functions, Given a candidate v(k) := [ w(k) * c(k) * ] * , define the Runge-Kutta F(·) function as F(·) can be easily obtained using a linear solver such as LU decomposition. The inverse M −1 is never calculated explicitly. The algorithm requires an initial v 0 = v(k 0 ) calculated using CL-c. As in [29], the Dormand-Prince RK5(4)7M method [15, p. 23] and Hermite interpolation strategy of Shampine [42, p. 148] are used. We use automatic step-size control as described in [22, p. 167]. For an interval [k ( j) , k ( j+1) ] with midpoint k (mid) , the integrator produces ). Thus, after numerical integration, a solution set of v ( j) ,v ( j) , and v ( j+1/2) is obtained upon which piecewise Hermite interpolation can be performed. If both c(k) and the eigenvector w(k) is required then interpolation is over N z length vectors; if only c(k) is required then interpolation is only one-dimensional. Example output is shown in Fig. 6 (one-dimensional output).

Nonsimple Eigenvalues and Critical Layer Processing for PF-R-{r,a}-c
As alluded to in Sect. 4.2, for c / ∈ [U min , U max ] then we know that for this problem the eigenvalues remain simple. However, there is essential spectrum in the interval a b For the path-following method, if c ∈ [U min , U max ] then M may be become nearly singular and it is possible to encounter nonsimple eigenvalues. Furthermore, the eigenvalues in this interval are sampled from the essential spectrum so they can be located very closely together; standard techniques to reliably process through this region would require a very small step-size indeed. Therefore, if within an integration step it is found that c ∈ [U min , U max ] then the path-following algorithm uses CL-c with the critical layer filtering described in Sect. 3.4 to calculate v(k) at each stage of the Runge-Kutta integration and then a linear solve to calculate the derivative as usual. This is obviously slower but a large step-size may be used and it ensures that the path-following algorithm tracks the sought c + (k) curve.

PFmp-R-{r,a}-c: Improving Accuracy for PF-R-{r,a}-c
As shall be described in Sect. 5, the error in the CL-methods is determined almost entirely by roundoff error incurred during the solution of the quadratic eigenvalue problem in double precision. In our numerical tests, the path-following algorithm essentially maintains the same error as is present in the initial v 0 . By calculating v(k 0 ) in high-precision arithmetic then executing the path-following schemes in double precision as normal, an improvement in accuracy of two to three orders of magnitude is obtained. This is discussed further in Sect. 5.1.

Path-Following Method for Forward General Problem (PF-G-c)
The PF-R-r-c and PF-R-a-c algorithms can be combined to create an efficient algorithm that can process scattered data query points, which we describe below.
(i) First, PF-R-a-c is executed at some nominal k = k 0 and interpolation points at angles {θ ( j) } J j=1 are calculated. See Fig. 7a, b. a c e b d f Fig. 7 Steps of PF-G. Intervals for k x and k y were chosen for purposes for visual clarity only; in practice, these are chosen as is appropriate for the problem. a Planar plot of PF-R-a used at radius k 0 , interpolated at angles θ ( j) . Blue circles are control points, red asterisks are interpolation points. Angles indicated in dotted grey. b 3d plot as panel (a). c Planar plot of PF-R-r used along each θ ( j) . Blue circles are control points, red asterisks are interpolation points (the k (i) ). d 3d plot as panel (c). e Planar plot showing interpolation of query point. Magenta crosses indicate the k q radius on each radial slice at angles θ (l) , θ (l+1) . Red asterisk indicates the query point (k q , θ q ). f 3d plot as panel (e) (ii) The results from step i. are used as the initial v 0 values for PF-R-r-c calculating radially along each θ ( j) . The curve on each radial slice is then interpolated at predefined {k (i) } I i=1 points. The control points for each radial slice are then replaced with the control points at these fixed k (i) (we do not calculate new midpoint values). So there is now a two-dimensional polar grid at angles θ ( j) and radii k (i) . See arrangement in Fig. 7c, d. (iii) For an arbitrary query point (k q , θ q ), the nearest angles θ (l) , θ (l+1) and radii k (m) , k (m+1) are identified. The interpolant on radial slices at angles θ (l) , θ (l+1) are calculated at radius k q . Equation (35) is then used to calculate the angular derivatives. Finally, (cubic) interpolation is performed in an angular direction for angle θ q to obtain the solution. See Fig. 7e, f.
Note that after steps i. and ii. are calculated once, only step iii. needs to be performed for further query points, in a similar manner as PF-R-{r,a}-c.
There is a loss in accuracy because of the required use of cubic interpolation-due to not having the midpoint-rather than the 4th-order interpolation used in PF-R-{r,a}c. However, for these purposes, it is not particularly significant. For clarity, we omit further analysis of PF-G: it is broadly similar to PF-R-r-c and does not add anything to the discussion.

Convergence and Error Estimates
It is well known that for sufficiently smooth solutions, spectral methods converge exponentially fast or with 'spectral accuracy' [11, ch. 1,2]. However, roundoff error poses a significant challenge for collocation methods due to the interaction of illconditioned matrices with commonly used double precision calculations [6,11]. We adopt a heuristic strategy to estimate the accuracy of each algorithm.

Dependence of Eigenvalue Accuracy on Order N z
To determine accuracy depending on N z , we first calculate a reference dispersion relation This is done in high-precision arithmetic, using the Advanpix library [1], for N ref = 384; this size of matrix exceeds what would be used in practice.
We calculate the relative normwise error in a candidate dispersion relation R cand (with k This is done for for the U T velocity profile using the CL-c, PF-R-r-c, PFmp-R-r-c, and DIM algorithms. The CL and PF methods reduce error with spectral accuracy until roundoff error starts to dominate. The high-precision initial calculation for the PFmp algorithm avoids this roundoff error and it can be seen that the path-following method a b Fig. 8 Plots of error in eigenvalue computations and backwards error + condition number estimates. a Loglog plot of normwise relative error in candidate algorithms depending on N z . The collocation method and path-following algorithm reduce error with spectral accuracy until roundoff error starts to dominate around N z ≈ 65. Path following+MP which uses a high-precision initial result maintains broadly the same low error despite the actual path-following calculations being performed in double precision (as is expected, the tolerance for the adaptive step-size control in the Dormand-Prince integration had to be set to around 10 −15 to achieve this). DIM reduces error as predicted, as O(N −2 z ). b Backwards error, condition number, and error estimate (backwards error × condition number). The time series for the quadratic eigenproblem solves for CL shown in blue, the linear solves for PF in magenta. From the condition numbers, shown with dotted lines, it can be seen that CL solves are of the same order of magnitude. The backwards error, shown in dashed lines, shows the linear solves in PF are appreciably more backwards stable. The error estimate clearly favours PF itself retains this improved accuracy even in double precision. DIM is included for indicative purposes. See Fig. 8a.
A possible explanation for this can be found in comparison of the backwards error and conditioning of the quadratic eigenvalue solve used for the CL-schemes compared to the linear solves predominantly used in PF-schemes. Although it is not a direct comparison-the linear solves are used to calculate a derivative, not the value itselfit may lend some insight. The backward error for the linear solve can be calculated with [23, eqn. 1.2] and the condition number in the usual manner: where w l is a corresponding left-eigenvector. Noting the usual inequality [51, eqn. For a range of N z , we calculate the · ∞ of the backwards errors η L , η Q and condition numbers κ L , κ Q over the k vector. This then permits calculating the product from Eq. (41). This is shown in Fig. 8b. The condition numbers are of the same magnitude, κ L ≈ κ Q , but the backwards error for the linear solves in PF-R-r-c clearly smaller, η L η Q . Although a direct comparison cannot be made, this suggests the path-following method has favourable numerical properties.

k-Dependent Convergence
As can be seen from Fig. 9a, the eigenvectors become numerically singular at the surface as k increases, implying that increasingly many basis polynomials are required to approximate the solution. This can be tested using a similar algorithm as in [5] to determine when the Chebyshev series has converged. Specifically, we calculate an envelope then use a histogram to locate the plateau convergence region. We then determine the required N z to reach convergence for a range of k values, as shown in Fig. 9b. For shorter wavelengths, much higher N z is required to reach convergence and so requiring more computational resources. This problem can be entirely ameliorated, as described in Sect. 6.

Adaptive Depth and Partition of Unity
It is clear from the results in Sect. 5.2 that as k increases, the required N z becomes infeasibly large due to the singular behaviour of the eigenfunction. This can be avoided by using a smaller h so that h 1 for higher k, on the following rationale. We expect that the eigenfunction decays roughly as e kz . Therefore, we can estimate the depth below which the eigenfunction is effectively zero, for numerical purposes. Let δ be the tolerance below which numerical values are considered zero, e.g. the "machine epsilon". Let h δ (k) := min{1, − log(δ)/k} be an estimate of the depth, taking into account the finite depth, at which the eigenfunction decays below tolerance δ for a given k. The CL-r scheme can be adapted for each calculated k. For a given k, we can set h = h δ (k). The calculated eigenvalue for the phase velocity will be correct automatically. The eigenvector may be mapped back onto the original interval on any suitably large set of z points chosen on the [−1, 0] interval using barycentric interpolation [7]; the eigenvector will be zero for z < −h δ (k).
This procedure becomes less obvious when considering the PF-r scheme because it would require mapping the entire system at each Runge-Kutta step. To avoid this, we instead split the k domain into several, partially overlapping, subintervals for which the chosen depth is suitable for all k in that subinterval. The path-following algorithm is then used on each subinterval independently with the appropriate depth. To combine the subintervals and avoid loss of smoothness in the computed dispersion relation, a partition of unity method is used on the overlaps.
We seek a scheme to choose subintervals I b ] and corresponding depths h ( j) that is both simple and easy to implement. For some k ( j) ∈ I where C min , C max are constants controlling the proportion of the [−h ( j) , 0] interval that −h δ (k ( j) ) should be within. In our computations, we found that C min = 0.3 and C max = 0.8 worked well. The subintervals and associated depths are then calculated as This generates intervals I ( j) k in such a manner that k , i.e. there is some overlap in the intervals. We use the partition of unity method described in [2] to join the subintervals. This is demonstrated in Fig. 10.

Performance Analysis
We compare the performance of the path-following algorithm against the collocation scheme CL-c from Sect. 3.3 and also some methods that were referenced in Sect. 1.2: DIM [27], the naive shooting method (SH), and the high-order interpolant method (CHEB) using Chebfun [17] with collocation for pointwise sampling.

Asymptotic Performance Analysis
There are two variables which control the expected computation time for the candidate algorithms: the number of z evaluation points, N z , and the number of query points N q . Since N z determines accuracy and is dependent on algorithm choice, we assume N z is set appropriately for each algorithm to achieve similar accuracy. Therefore, our primary concern shall be how the algorithms scale with N q .
DIM will incur a fixed per-point computational cost that depends on N z , which we denote σ DIM (N z ). So, the expected cost is O(σ DIM (N z )N q ). Similarly, the collocation algorithm incurs a fixed per-point computational cost which also depends on N z , σ CL (N z ). So, the expected cost is O(σ CL (N z )N q ). A naive shooting method would also have an expected cost of the form O(σ SH (N z )N q ). These estimates are valid for both the reduced and general problems.
In contrast, the reduced path-following algorithm incurs an initial computational cost dependent on N z , σ PF-NI (N z ), whereafter there is a very light-weight perpoint cost, σ PF-Q σ PF-NI (N z ). Therefore, the expected computational cost is O(σ PF-NI (N z ) + σ PF-Q N q ). The general path-following algorithm is similar with only the coefficients changed.
Finally, constructing a high-order interpolant using Chebfun with the collocation method as the pointwise sampling would incur an initial setup cost dependent on both N z and required number of sampling points N s , σ CHEB-S (N z , N s ); thereafter, there is a very low per-point cost, σ CHEB-Q σ CHEB-S (N z , N s ). So the expected computation This is summarized in the following table, assuming the eigenvector output is not required: It immediately becomes clear that if σ PF-NI (N z ) is not too large and σ PF-Q is sufficiently small then as N q increases, the path-following algorithms are the most efficient.
Analogous behaviour can be expected from the high-order interpolant with Chebfun but numerical experiments show that the initial setup cost σ CHEB-S (N z , N s ) is prohibitively high. This is discussed at further length in Sect. 7.2 below.

Numerical Experiments to Determine Algorithm Performance
The asymptotic complexity claims are confirmed by practical testing. For clarity, we only test with the reduced problem. By measuring the time taken for each algorithm to compute the dispersion relation for differing N q , we can determine the computational complexity in relation to N q as shown in the log-log plot, Fig. 11. Each algorithm can be executed with different parameter choices that influence accuracy. As such, we calibrated algorithms to produce output at three different accuracies (measured as relative normwise error using Eq. (38)): ≈ 10 −4 ('low' accuracy), ≈ 10 −7 ('medium' accuracy), and ≈ 10 −10 ('high' accuracy). The naive shooting algorithm (SH) is tested only at the medium and high accuracy levels. The high-order interpolant method (CHEB) is only performed for the high accuracy level.
As seen in the results, the path-following algorithm is asymptotically at least two orders of magnitude faster than DIM, SH, and the collocation scheme. The breakeven point in N q at which the path-following scheme becomes faster than DIM is N q ≈ 1000 for the low accuracy solution; for the medium and high accuracy levels, the path-following method is always faster.
Comparing the path-following algorithm against the only other scheme that uses interpolation, CHEB, it can be seen that they share similar characteristics. However, the initial setup cost for the path-following scheme is much lower than for the highorder interpolant scheme. This is because, in this setup, a quadratic eigenvalue problem must be solved for each point on the curve used to construct the high-order interpolant. Each quadratic eigenproblem solve incurs a QZ decomposition on a size 2N z matrix. The path-following algorithm avoids most of this cost by swapping the eigenproblem solves for linear solves on size N z matrices, which is considerably faster.
There is an issue with the interpolation for the path-following scheme that can be seen after around 10 5 points: the computation time starts to increase linearly. This is due to having to determine what the relevant control points are to use for the Hermite interpolation for a given query point. In the implementation used, this was done using MATLAB's histcounts() function. If some ordering of the query points were assumed, this cost could be eliminated but we do not assume order in the query points.

Guidance on Optimal Parameter Choices
Optimal parameter choices are predicated on two key properties: the required accuracy and the anticipated number of query points.
As can be observed from Fig. 11, the path-following algorithm is most effective when higher accuracy and at least a moderate number of query points are required. The nominal setup cost caused by the initial quadratic eigenproblem solve and numerical integration is dependent on the order of differentiation matrix used, N z . So this should Fig. 11 Performance plot for reduced problem: log-log plot of computational time against number of query points. The algorithms were tested at various relative normwise accuracies: ε ≈ 10 −4 ('low' accuracy), ε ≈ 10 −7 ('medium' accuracy), and ε ≈ 10 −10 ('high' accuracy). Note that the DIM, shooting (SH), and collocation algorithms are clearly linear in complexity with respect to N q . The high-order interpolant using Chebfun with collocation sampling (CHEB) is almost constant time but has a large initial computational cost. The path-following algorithm is also almost constant time but has a much lower initial computational cost due to swapping quadratic eigenproblem solves for linear solves be kept at the lowest value possible that maintains required accuracy. We found N z between 48 and 64 is optimal for the cases we tested. Furthermore, using N z too high risks roundoff error causing deleterious effects, cf. Fig. 8a.
The Dormand-Prince integrator requires a tolerance for the adaptive step-size control. We suggest that 10 −11 is the smallest value to use when the initial eigenproblem solve is performed in double precision. If the initial eigenvalue solve can be performed more accurately, for example in high-precision arithmetic, then the tolerance can be set around 10 −15 . In any case, if using a smaller N z then the tolerance should be adjusted to match the accuracy from the collocation solution.

Conclusions
By considering the boundary value eigenproblem posed by the Rayleigh stability equation with linearized free-surface boundary condition as parameterized by wave number k, we can adapt the path-following scheme in [29] to efficiently calculate the dispersion relation at high accuracy. This efficiency is achieved by first exchanging many expensive QZ decompositions on a size 2N matrix for one QZ decomposition and some linear solves on a size N matrix; second, we 'front load' the computational The accuracy tests in Sect. 5 suggest that the path-following algorithm can maintain the accuracy of the initial eigenpair v 0 and appears to be numerically more stable than the QZ decomposition used to obtain the initial eigenpair.
The algorithm is extended to permit calculation in the two-dimensional k-plane with scattered data. A method for processing through critical layers is given.
In other tests, not included here, it is clear the same approach works well for other problems from physics and engineering, assuming the problem is parametrized by a real scalar. Additional difficulties arise when there are exceptional points or bifurcations in the solution curve, or if the ODEs become stiff. These challenges may form the basis of future work.
The MATLAB library used to perform the calculations used in this paper is maintained at [31].

Example Usage
To demonstrate a possible use of the path-following algorithm, we calculate the result of a simple initial value problem of the form in [56, sec. 11.2] using the shear profile U(z) = (U T (z), 0) with F 2 = 0.15. Upon a sheared current, the anisotropic dispersion relation results in a ring wave pattern which is visibly asymmetric, predicted for a linear current [19] and recently observed experimentally for the first time [46].
Solving for ζ + and ζ − : In our example, we impose an initial Gaussian of the formζ 0 (x) = e − x 2 2σ 2 − y 2 2σ 2 andζ 0 (x) = 0. Specifically, we choose σ 2 ≈ 5/2 and plot for t = 16; this is shown in Fig. 12. The dispersion relation, ω ± (k) is calculated using the path-following algorithm, which allows us to pick a large number of grid points in x coordinates and still compute the result quickly.