## Abstract

This paper presents a method of analysing the dispersion relation and field shape of any type of wave field for which the dispersion relation is transcendental. The method involves replacing each transcendental term in the dispersion relation by a finite-product polynomial. The finite products chosen must be consistent with the low-frequency, low-wavenumber limit; but the method is nevertheless accurate up to high frequencies and high wavenumbers. Full details of the method are presented for a non-trivial example, that of anti-symmetric elastic waves in a layer; the method gives a sequence of polynomial approximations to the dispersion relation of extraordinary accuracy over an enormous range of frequencies and wavenumbers. It is proved that the method is accurate because certain gamma-function expressions, which occur as ratios of transcendental terms to finite products, largely cancel out, nullifying Runge’s phenomenon. The polynomial approximations, which are unrelated to Taylor series, introduce no spurious branches into the dispersion relation, and are ideal for numerical computation. The method is potentially useful for a very wide range of problems in wave theory and stability theory.

## 1. Introduction

The aim of this paper is to present a practical method for obtaining detailed information about wave fields governed by a dispersion relation that contains transcendental terms. In principle, such information can be obtained by numerical computation, but in practice the need to account for multiple branches in a complex frequency or wavenumber space leads to difficulties in producing robust computer code. This problem arises particularly in the modelling of complex systems, where a wave field of one type can be strongly coupled to a wave field of a different type. Then it is desirable to have the simplest possible representation of the wave field in each component system.

Of all representations, polynomials are the simplest. Modern software packages make numerical work with polynomials very rapid, even when the order of the polynomials is high. But the polynomials used must be fit for their purpose. In particular, truncated Taylor series may have a limited range of accuracy or involve massive cancellation or introduce spurious roots. An alternative is a finite-product polynomial approximation which, in a bounded region, has the same roots as the original function. The difficulty here is Runge’s phenomenon (Trefethen 2000, p. 42), i.e. unwanted oscillations similar to, but worse than, those of Gibbs’ phenomenon in Fourier series. The use of Chebyshev polynomials removes this difficulty, but a disadvantage for some purposes is that such polynomials hide the analytical structure of the underlying problem.

The mathematical idea behind the present paper is that Runge’s phenomenon can be eliminated by a method based not on Chebyshev polynomials but on two formulae involving finite products. The formulae, though elementary, are omitted from some standard reference works, e.g. Abramowitz & Stegun (1965) and Olver (1974), but are given in, for example, Erdélyi *et al.* (1953, p. 4) and Magnus *et al.* (1966, p. 2). These formulae, expressed in terms of a factor *g*_{m}(*z*) defined by
1.1
are
1.2
1.3
Here *m*=0,1,2,…, with the convention that when *m*=0 the products are assumed equal to 1. We use the factorial notation, namely *z*!=*Γ*(*z*+1) for all *z*, where *Γ* denotes the gamma function. The factor *g*_{m} will be called the gamma-conversion factor, because by equations (1.2) and (1.3) it provides a method of conversion, involving gamma functions, between trigonometric and polynomial functions. Thus application of equations (1.2) and (1.3) may be called gamma conversion. In this paper, we proceed from sines and cosines to polynomials. In a different type of problem, equations (1.2) and (1.3) may be used to express a finite product as an amplitude-modulated sine or cosine; hence Runge’s phenomenon may be regarded as amplitude modulation governed by gamma functions.

The crucial fact underlying the finite-product method is that certain combinations of gamma-conversion factors cancel out almost exactly; hence in a very wide class of problems in waves and stability, these factors may be ignored from the outset, i.e. taken equal to 1, and all trigonometric or hyperbolic functions in the dispersion relation may be replaced by elementary finite products from the start. That is, Runge’s phenomenon is simply ignored, because it cancels out almost exactly at the end. The cancellation may be verified by numerical evaluation of the conversion factors or by Stirling’s approximation (Abramowitz & Stegun 1965, p. 257), which gives (for *m*≠0)
1.4
This approximation to *g*_{m}(*z*) is extremely accurate for |*z*| up to *m*−1/2.

The paper is organized as follows. Section 2 summarizes the required theory for a non-trivial example of wave propagation, that of anti-symmetric elastic waves in a layer. The result is a ‘two-tangent’ dispersion relation known to have a rich analytical structure. Section 3 introduces a family of finite-product approximations and identifies, by comparison with numerical computations, a subfamily of extremely high accuracy in representing the dispersion relation. Section 4 gives an analytical theory, based on the gamma-conversion formulae (1.2) and (1.3), relating the finite-product approximations to the exact dispersion relation and uses this theory to provide a refined accuracy analysis; a complete quantitative explanation is found for why the accuracy is so high. Section 5, the conclusion, explains why the method of finite-product approximation is promising for a variety of problems in wave theory and stability.

## 2. Rayleigh–Lamb theory

### (a) Governing equations

An elastic medium of thickness *h* is assumed to occupy the layer , −*h*/2<*y*<*h*/2 and support elastic waves satisfying the linear equations (Achenbach 1973, pp. 59, 257)
2.1
Here (*u*(*x*,*y*,*t*),*v*(*x*,*y*,*t*)) are longitudinal and transverse displacements at position (*x*,*y*) and time *t*, subscripts *x*,*y* and *t* denote partial differentiation, *c*_{1} is the *P*-wave speed, i.e. compression-wave speed, and *c*_{2} is the *S*-wave speed, i.e. the shear-wave speed. The elastic medium has Young’s modulus *E*, Poisson’s ratio *ν*, undisturbed density *ρ* and ‘reference speed’ *c*_{0}=(*E*/*ρ*)^{1/2}. For an elastic medium in plane strain, the wave speeds are
2.2
For plane stress, the value of *c*_{1} here must be replaced by *c*_{0}/(1−*ν*^{2})^{1/2}. In numerical calculations, we shall take *ν*=0.3. Thus for plane strain, we have *c*_{1}=1.16*c*_{0} and *c*_{2}=0.62*c*_{0}, and for plane stress, we have *c*_{1}=1.05*c*_{0} and *c*_{2}=0.62*c*_{0}. The normal stresses *τ*_{xx} and *τ*_{yy} and shear stress *τ*_{xy} corresponding to displacements (*u*,*v*) are
2.3

### (b) The Rayleigh–Lamb dispersion relation

This paper gives details for anti-symmetric (‘bending-wave’) solutions of the above equations, i.e. waves such that (*u*,*v*) are (odd, even) in the transverse coordinate *y*. Accordingly, solutions with real frequency *ω* and longitudinal wavenumber *k* are sought in the form
2.4
where *l* is the transverse wavenumber, possibly complex, and the sum is over suitable solutions *l* of the equation obtained by substituting equation (2.4) into equation (2.1) and requiring that *A*(*l*) and *B*(*l*) are both not zero. This equation factorizes to give the result that *l* must satisfy one of the equations *l*^{2}=(*ω*/*c*_{1})^{2}−*k*^{2} or *l*^{2}=(*ω*/*c*_{2})^{2}−*k*^{2}, corresponding to *P*-waves and *S*-waves. The values of *A*(*l*) and *B*(*l*) are found from the traction-free boundary conditions, obtained by equating to zero at *y*=±*h*/2 the stresses *τ*_{yy},*τ*_{xy} in equation (2.3). In the dimensionless variables
2.5
in terms of which the elastic layer occupies the region |*Y* |≤1/2, the displacement field is
2.6
where
2.7
and
2.8
with
2.9
The stresses (*τ*_{xy},*τ*_{yy}) corresponding to equation (2.6) are
2.10
Throughout this paper, a subscript 1 refers to *P*-waves, and a subscript 2 refers to *S*-waves. For complex *K*, we write *K*=*K*_{r}+i*K*_{i}; if *K* is real by implication, we usually write *K* rather than *K*_{r}. Most formulae will be written in terms of *c*_{1} and *c*_{2}; but it is convenient to define *Ω* in terms of *c*_{0} (which does not depend on *ν*) for comparison of our plots with many of those in the engineering literature.

For plane strain with Poisson’s ratio *ν*=0.3, the real and imaginary branches of the Rayleigh–Lamb dispersion relation (2.7) are plotted for *Ω* and *K* or *K*_{i} up to 40 as the solid lines in figure 1*a*,*b*. The lower left corner of each plot is the Euler–Bernoulli region, for which a Taylor-series expansion of the dispersion relation (2.7) gives
2.11
Here *c*_{B}, the bending-wave characteristic speed, is defined by . For plane strain, *c*_{B}=*c*_{0}/{12(1−*ν*^{2})}^{1/2}, i.e. 0.303*c*_{0} when *ν*=0.3; for plane stress, . Timoshenko’s approximation (Achenbach 1973, p. 253; Howe 1998, p. 21), containing a parameter *κ* often chosen to be *π*^{2}/12, is
2.12

## 3. Finite products

### (a) Finite products and their numerical accuracy

We shall need two families of polynomials, *S*_{m}(*s*) and *C*_{m}(*s*), defined for *m*=0 by *S*_{0}(*s*)=1 and *C*_{0}(*s*)=1 for all *s*, and for *m*=1,2,… by
3.1
Here *S*_{m} is ‘sine-based’ and *C*_{m}(*s*) is ‘cosine-based’, because and as . A finite-product approximation to the dispersion relation (2.7) is
3.2
where *m*_{1},*m*_{2},*n*_{1} and *n*_{2} are at our disposal. On the left-hand side of equation (3.2), the numerator is an approximation to and the denominator is an approximation to . At fixed *L*_{1} and *L*_{2}, these approximations become progressively more accurate as *m*_{1},*m*_{2},*n*_{1} and *n*_{2} are increased. However, this fact by itself is of limited interest, because we aim to use equation (3.2) at fixed values of *m*_{1},*m*_{2},*n*_{1} and *n*_{2}. Then the crucial question is: for fixed *m*_{1},*m*_{2},*n*_{1} and *n*_{2}, in which regions of (*Ω*,*K*) space, where *Ω* or *K* may be complex, is equation (3.2) accurate? There are grounds here for pessimism. For example, if a plot of for 0≤*z*≤5*π* is superposed with plots of for various values of *m*, it is found that *m* has to be enormous before high accuracy is obtained for the whole range of *z* up to 5*π*. This is a manifestation of Runge’s phenomenon.

It is therefore a surprise that there are small values of *m*_{1},*m*_{2},*n*_{1} and *n*_{2} for which equation (3.2) captures whole branches of the dispersion relation with great precision up to arbitrarily large values of *Ω* and *K*. The requirements are (i) *m*_{1}=*m*_{2} (=*m*, say) and *n*_{1}=*n*_{2} (=*n*, say), to capture the behaviour at small *Ω* and *K* and (ii) *n*=*m*+1, to capture the behaviour at large *Ω* and *K*. This gives a discrete one-parameter family of polynomial approximations to the dispersion relation, labelled by (*m*,*n*)=(0,1),(1,2),(2,3),…. For brevity, we refer to the (0,1) finite-product approximation, the (1,2) finite-product approximation, etc. On cross-multiplying equation (3.2), we obtain a polynomial equation in and , and hence a polynomial equation in *Ω*^{2} and *K*^{2}, on using the definitions of and in equation (2.9). It is therefore a trivial matter computationally to find, for example, all the real and complex roots *K* for given *Ω*, using a polynomial root-finding subroutine in any standard numerical software system, for example, Matlab. It will be shown in equations (3.3)–(3.6) that the (*m*,*n*) finite-product approximation is of degree *m*+*n*+1 in *Ω*^{2} and *K*^{2}.

For plane strain with Poisson’s ratio *ν*=0.3, the real and imaginary branches of the finite-product approximations (*m*,*n*)=(0,1),(1,2),(2,3), (8,9) for real *Ω* are plotted as dotted lines in figure 2*a*–*h*; the exact dispersion relation is superposed as solid lines. The (1,2) approximation, in figure 2*c*,*d*, is accurate (e.g. has error at most 1%) in almost all of a square region defined by *Ω* up to 6 and *K* or *K*_{i} up to 6; moreover, the real branch through the first symbol ∇ on the *Ω*-axis in figure 2*c* is captured accurately over its whole infinite length, i.e. up to arbitrarily large (*Ω*,*K*), not shown in the figure. The (2,3) approximation, in figure 2*e*,*f*, is accurate in almost all of a square region defined by *Ω* up to 10 and *K* or *K*_{i} up to 10, and the two real branches through the first two symbols ∇ on the *Ω*-axis in figure 2*e* are captured accurately over their whole infinite length. The (8,9) approximation, in figure 2*g*,*h*, is accurate in almost all of a square region defined by *Ω* up to 35 and *K* or *K*_{i} up to 35, and the eight real branches through the first five symbols ∇ and the first three symbols * on the *Ω*-axis in figure 2*g* are captured accurately over their whole infinite length. In general, the (*m*,*m*+1) approximation captures *m* branches accurately over their whole infinite length. This is a remarkable achievement for a polynomial approximation to a transcendental equation. In many applied problems, dimensionless frequencies as high as *Ω*=6 are already very high indeed. Therefore, the (1,2) approximation is ideal for many practical purposes, and in particular offers a massive improvement over attempts to extend Timoshenko-type approximations (2.12) to ever higher frequencies.

If extreme accuracy (e.g. error at most 0.01%) is needed for very large *Ω* and *K*, the value of *m* required is not too high. For example, if *m*=20, the exact and approximate dispersion curves are indistinguishable by eye in a square region defined by *Ω* up to about 30 and *K* or *K*_{i} up to about 30; if *m*=30, the curves are indistinguishable for *Ω* up to about 60 and *K* or *K*_{i} up to about 60. Even for these values of *m*, the numerical computation, involving only polynomials, is almost instantaneous on any computer. Therefore, in practice, unlimited numerical accuracy is obtainable at negligible cost from approximations of the type (3.2). A remarkable feature of these polynomial approximations to a transcendental dispersion relation is that they introduce no spurious roots or branches; the only error is that, for any given *m*, some of the branches become inaccurate at high enough *Ω* or *K*. This is in sharp contrast to, for example, Taylor series truncated at high order.

The method presented in this section is sufficient for the numerical analysis of a very wide range of dispersion relations if the sole aim is ease of computation. The rest of the paper gives a detailed mathematical analysis of why an equation of such a simple type as (3.2), in which a Chebyshev polynomial is nowhere to be seen, is able to overwhelm Runge’s phenomenon. The analysis yields, as a by-product, accurate mathematical formulae for the difference between the exact and approximate curves plotted in figure 2.

### (b) The low-frequency limit

In the finite-product approximation (3.2), we have insisted that *m*_{1}=*m*_{2}=*m* and *n*_{1}=*n*_{2}=*n*. Such a choice is not obvious, because each original tangent term in the dispersion relation (2.7) is given different numbers of factors in its sine numerator and cosine denominator. Nevertheless, the choice *m*_{1}=*m*_{2} and *n*_{1}=*n*_{2} is essential, because it is the only way to eliminate the term in *k*^{2} from the low-frequency limit of equation (3.2); this limit must reduce to a balance between terms in *ω*^{2} and *k*^{4}, to be consistent with the low-frequency limit of equation (2.7), i.e. the Euler–Bernoulli bending-wave limit.

To prove the above remark, note from equation (2.9) that and contain (*k**h*)^{2} with identical coefficients, −1. Hence in the ratio , the factors in equation (3.1)_{1} must pair off exactly; i.e. each factor 1−*s*/(*m*′*π*)^{2} for a given *m*′ must either occur in the numerator with and also in the denominator with or occur in neither the numerator nor the denominator. This forces *m*_{1}=*m*_{2}. The same argument applied to forces *n*_{1}=*n*_{2}.

### (c) The structure of the (*m*,*n*) polynomial and of the field

We henceforth assume that *m*_{1}=*m*_{2}=*m* and *n*_{1}=*n*_{2}=*n*. On clearing the finite-product approximation (3.2) of fractions, and inserting the definitions (2.9) of , and , the result is a polynomial in *ω*^{2} and *k*^{2}. We shall refer to such a polynomial by its degree in *ω*^{2} and *k*^{2} rather than in *ω* and *k*. Then at first sight, the degree of the polynomial obtained from equation (3.2) is *m*+*n*+2. However, the polynomial is divisible by . To see this, we subtract 1 from each side of equation (3.2) and use the identities
3.3
and
3.4
Here , where *c*_{B} is the bending-wave characteristic speed as defined after equation (2.11), so that the Euler–Bernoulli limit is (*ω**h*/*c*_{B})^{2} = (*k**h*)^{4}. Then equation (3.2) becomes
3.5
where
3.6
Here the factor divides the numerator, so that *I*_{m,n} is a polynomial of degree *m*+*n*−1 in and , or in *ω*^{2} and *k*^{2}. Hence equation (3.5) is of degree *m*+*n*+1 in *ω*^{2} and *k*^{2}; it is the (*m*,*n*) finite-product approximation to the exact dispersion relation (2.7).

Corresponding to the (*m*,*n*) approximation are the exact displacement field
3.7
and
3.8
and the exact stress field
3.9
and
3.10
Here a factor is omitted. These expressions are obtained from displacement (2.6) and stress (2.10) using the gamma-conversion formulae (1.2) and (1.3) and finite products (3.1). The expressions may be simplified by means of Stirling’s approximation (1.4); since the elastic layer occupies the region |*Y* |≤1/2, the remark after equation (1.4) shows that the approximation is extremely accurate for the larger of |*L*_{1}| and |*L*_{2}| up to the smaller of 2*π*(*m*−1/2) and 2*π*(*n*−1). The corresponding regions in the (*Ω*,*K*) plane, for *K* real or imaginary, may be determined from simple plots of constant |*L*_{1}| and |*L*_{2}|, which are hyperbolae or ellipses. Expressions (3.7)–(3.10) make explicit the way in which the gamma-conversion factors *g*_{m} and *g*_{n−1/2} nullify Runge’s phenomenon in the finite-product polynomials *S*_{m} and *C*_{n}.

### (d) The high-frequency limit

If we retain only the terms of highest degree in equation (3.5), and divide by a common coefficient, the result is
3.11
This equation determines the asymptotes of the branches of equation (3.5) at high frequencies and wavenumbers. The particular cases *n*=*m* and *n*=*m*+1 of equation (3.11) are
3.12
3.13
More generally, if *n*=*m*±*r*, with *r*≥1, then equation (3.11) is
3.14
and
3.15
Here the power of outside the braces, i.e. *n* in equations (3.12) and (3.14) and *m* in equation (3.13) and (3.15), is the number of branches of equation (3.5) which tend at large *ω* and *k* to the *P*-wave line *ω**h*/*c*_{1}=*k**h*; the power of outside the braces, i.e. *n* in equation (3.12) and (3.14) and *m*+1 in equations (3.13) and (3.15), is the number of branches of equation (3.5) which tend to the *S*-wave line *ω**h*/*c*_{2}=*k**h*. The number of real factors of the terms in braces determines the number of the remaining branches of equation (3.5) which have real asymptotes. This last number varies irregularly with *r*, as do the directions of the corresponding asymptotes, because these features depend on the accidental occurrence of real roots of the expressions in braces in equations (3.14) and (3.15) regarded as functions of *ω* and *k*.

By contrast, for the exact dispersion relation (2.7), plotted for real *ω* and *k* in figure 1*a*, all of the branches except the branch through the origin tend to the *S*-wave line *ω**h*/*c*_{2}=*k**h* at large *ω* and *k*, outside the region plotted. The *S*-wave line is the dashed line OS in the figure, where O will always denote the origin. Moreover, at large but not very large values of *ω* and *k*, all of these branches are close to the *P*-wave line *ω**h*/*c*_{1}=*k**h*, in the inflection region of the real (*ω*,*k*) plane where the branches change from being concave downwards to concave upwards. The *P*-wave line is the dashed line OP in the figure.

We may therefore expect that for a finite-product approximation (3.5) of given degree *m*+*n*+1 in *ω*^{2} and *k*^{2}, the greatest accuracy at medium and large *ω* and *k* is attained by maximizing the number of factors of and in the terms of highest degree, i.e. by choosing *n*=*m* or *n*=*m*+1 to attain the form (3.12) or (3.13). Here, an ‘accidental’ advantage of equation (3.13) emerges: the value of is close to *c*_{1}, and so gives a ‘bonus’ factor approximately equal to ; then equation (3.13) is factored completely into terms advantageous for numerical accuracy and is superior to equation (3.12), which has a spurious asymptote parallel to the *k*-axis. The expressions after (2.2) and (2.11) show that takes the value (1−2*ν*)^{1/2}/(1−*ν*) for plane strain and (1−*ν*^{2})^{1/2} for plane stress, i.e. 0.90 and 0.95 for Poisson’s ratio *ν*=0.3.

This analysis of the accuracy of the finite-product approximation (3.5) as a function of *m* and *n* is confirmed by computer plots for the range 0≤*m*≤10 and 0≤*n*≤10. The plots reveal the extraordinary accuracy consequent on choosing the highest-order terms to be of the form (3.13), i.e. on choosing *n*=*m*+1; plots for (*m*,*n*)=(0,1), (1,2), (2,3), (8,9) are shown in figure 2. The splitting of branches with the same asymptote into distinct well-defined branches, before the limiting behaviour at large *ω* and *k* is attained, is determined by the terms just below the highest degree in equation (3.5). The extremely high accuracy displayed in figure 2 suggests that for *n*=*m*+1 these splitting terms are very close to the corresponding splitting terms which could be computed from the exact dispersion relation (2.7).

Another remarkable feature of the finite-product approximations for *n*=*m*+1 is that they introduce no spurious branches anywhere in the complex *K* plane within a large ‘box’ containing the origin. This is in contrast to high-order Taylor-series approximations, which are beset by spurious branches in such a box. Thus for (*m*,*n*)=(1,2), within a three-dimensional box defined approximately by *Ω* up to 8 and *K*_{r}, *K*_{i} up to 16, the finite-product approximation captures all the branches of the dispersion relation, i.e. all wholly real, wholly imaginary and complex branches, including all bifurcations, and introduces no extra branches whatever. Hence within the box the topological structure of the branches is captured perfectly, and the accuracy is everywhere as high as revealed in the box cross sections shown for real and imaginary *K* in figure 2*c*,*d*. Box sizes for other (*m*,*n*) with *n*=*m*+1 can readily be determined; for example, for (*m*,*n*)=(0,1), corresponding to a Timoshenko-type approximation, a suitable box is *Ω*, *K*_{r}, *K*_{i} up to 4, and for (*m*,*n*)=(8,9), a suitable box is *Ω* up to 25 and *K*_{r}, *K*_{i} up to 40.

Extremely fine features of the exact dispersion curve can be resolved by the *n*=*m*+1 finite-product approximations without *m* needing to be very large. For example, as Poisson’s ratio *ν* is varied, branches can meet and connect up differently, changing their topology (e.g. Freedman 1990). The finest details of such a reconnection can be tracked with *m* at most about 30, and the polynomial calculations required take negligible time on any computer.

### (e) The Rayleigh-wave branch

For the exact dispersion relation (2.7), the branch through the origin tends at large *ω* and *k* to the Rayleigh-wave line *ω*=*c*_{R}*k*, where (Achenbach 1973, p. 189)
3.16
For typical Poisson’s ratios, *c*_{R} is about 10 per cent less than *c*_{2}; for example, *c*_{R}/*c*_{2}=0.93 for plane strain and *c*_{R}/*c*_{2}=0.92 for plane stress when *ν*=0.3. It follows from §3d that the finite-product approximation (3.5) does not capture the high-frequency limit of the exact Rayleigh-wave branch as accurately as it captures the high-frequency limit of the other branches. This is evident in figure 2*a*,*c*,*e*,*g*, which shows the gradual divergence of the exact and approximate Rayleigh-wave branches at high frequencies. Nevertheless, the branch through the origin is captured accurately to high, but not arbitrarily high, frequencies, well into the region where the displacement field is approximately a Rayleigh wave on each face of the elastic layer. In particular, the branch is captured accurately within the (*Ω*,*K*) boxes defined in the previous section. The possibility of ‘tuning’ the finite-product approximation by introducing an adjustable parameter into equation (3.5) is discussed in §4; the approximation may be tuned to give the Rayleigh-wave asymptote exactly.

### (f) The grid

Finite-product approximations agree with the exact dispersion relation on a grid of points. As this is a most important aspect of their accuracy, a detailed description of the grid will now be presented for the dispersion relation (2.7). The methods that follow apply generally; it will normally be worthwhile to determine the grid in complete detail.

Corresponding to the (*m*,*n*) finite-product approximation is the exact equation
3.17
where
3.18
and the gamma-conversion factors *g*_{m} and *g*_{n−1/2} are defined via equation (1.1). Equation (3.17) is the exact dispersion relation (2.7) written in gamma-converted form by means of mathematical identities (1.2) and (1.3) and the definitions (3.1) of *S*_{m}, *C*_{n}. Comparison of the exact equation (3.17) with the finite-product approximation (3.2) shows that the sole effect of the approximation is to replace *α* by 1. It follows from the structure of equation (3.17) that solutions of the finite-product equation agree with those of the exact equation at all points (*ω*,*k*) for which, simultaneously,
3.19
Inspection of the factors of *S*_{m} and *C*_{n} determined by equation (3.1) shows that in equation (3.19) there are nine cases to consider, obtained from the intersections of any one of (i) , (ii) , *m*′=1,…,*m*, and (iii) , *n*′=1,…,*n* with any one of (iv) *K*=0, (v) , *n*′=1,…,*n*, and (vi) , *m*′=1,…,*m*. Hence a prominent role is played by the ‘grid of bounds’ (Mindlin 1960) defined by the curves and for *r*=1,2,… and *r*′=1,2,…. The definitions of and in equation (2.9) show that for real *ω* these curves are hyperbolae when *k* is real, and ellipses when *k* is imaginary. Examination of the nine cases reveals that the solutions of equation (3.19) fall into six categories (the symbols on diagrams are indicated):

(i) the Euler–Bernoulli point (), at which (

*ω*,*k*)=(0,0);(ii) the

*P*-wave cut-on points (*), also called thickness-stretch resonances, at which (*ω**h*/*c*_{1},*k**h*)=(2*m*′*π*,0),*m*′=1,…,*m*;(iii) the

*S*-wave cut-on points (), also called thickness-shear resonances, at which (*ω**h*/*c*_{2},*k**h*)=((2*n*′−1)*π*,0),*n*′=1,…,*n*;(iv) the Lamé points (), at which , where and ;

(v) the

*m*(*m*−1)/2 even interior points (°), specified parametrically by*m*_{1}=1,…,*m*;*m*_{2}=1,…,*m*;*m*_{1}<*m*_{2}, at which 3.20 where , , and*k*may be real or imaginary; and(vi) The

*n*(*n*−1)/2 odd interior points (), specified parametrically by*n*_{1}=1,…,*n*;*n*_{2}=1,…,*n*;*n*_{1}<*n*_{2}, at which (with and ) 3.21 where , , and*k*may be real or imaginary.

These six categories are mutually exclusive and exhaustive, so that the total number of solutions of equation (3.19) is (*m*+1)(*m*+2)/2+*n*(*n*+1)/2, which for the preferred choice *n*=*m*+1 simplifies to *n*(*n*+1). Thus for (*m*,*n*)=(0,1), (1,2), (2,3), (8,9), the number of solutions is 2,6,12,90. The corresponding points, referred to as grid points, are indicated in all plots by their symbols.

The *m* Lamé points () lie on the straight line , i.e. on . This is the Lamé line, drawn as the dotted line *O**L* in figure 1*a*. At the Lamé points, the exact and approximate dispersion curves touch, i.e. have not only a common point but also a common tangent; this is a consequence of the term in both the exact and approximate equations. Since the (*m*,*n*) approximation captures *m* Lamé points (for any *n*), the tangency is part of the explanation for the remarkable accuracy of finite-product approximations at high frequencies and wavenumbers. The hyperbolae , also share this tangent at the Lamé points; arcs of these hyperbolae above the *P*-wave line OP are drawn as dash-dotted curves in figure 1*a*.

The wavenumbers *k* of the *m*(*m*−1)/2 even interior points (°) are real for *m*_{2}/*m*_{1}≥*c*_{1}/*c*_{2} and imaginary for *m*_{2}/*m*_{1}<*c*_{1}/*c*_{2}, by equation (3.20). Hence the exact numbers of real and imaginary wavenumbers vary somewhat irregularly with *c*_{1}/*c*_{2}, depending on the number of integer lattice points on each side of a line of slope *c*_{1}/*c*_{2}. A number-theoretic argument, based on placing a square box of known area around each lattice point, shows that the number of even interior points with real wavenumber is approximately (1/2)(*c*_{2}/*c*_{1}){*m*−(*c*_{1}−*c*_{2})/(2*c*_{2})}^{2}, so that the rest of the *m*(*m*−1)/2 even interior points have imaginary wavenumber. This approximation becomes less accurate as *c*_{1}/*c*_{2} becomes closer to 1, because in the limit *c*_{1}/*c*_{2}→1 (which is possible for plane stress but not plane strain) all the even interior points have real wavenumbers, by equation (3.20), but the approximate formula then gives *m*^{2}/2 not *m*(*m*−1)/2. Similarly, the wavenumbers of the *n*(*n*−1)/2 odd interior points () are real for (*n*_{2}−1/2)/(*n*_{1}−1/2)≥*c*_{1}/*c*_{2} and imaginary for (*n*_{2}−1/2)/(*n*_{1}−1/2)<*c*_{1}/*c*_{2}, by equation (3.21). A lattice-point argument shows that the number of real points is approximately (*c*_{2}/*c*_{1})*n*^{2}/2, and the number of imaginary points is approximately *n*(*n*−1)/2−(*c*_{2}/*c*_{1})*n*^{2}/2. Again, this approximation loses some accuracy as *c*_{1}/*c*_{2} approaches 1, because in this limit all the odd interior points have real wavenumbers, by equation (3.21). For plane strain with Poisson’s ratio *ν*=0.3, i.e. *c*_{1}/*c*_{2}=1.87, the number of (real, imaginary) even interior points for (*m*,*n*)=(0,1), (1,2), (2,3), (8,9) is (0,0), (0,0), (1,0), (16,12) from accurate computations, and (0.05,−0.05), (0.09,−0.09), (0.7,0.3), (15.3,12.7) from the approximate formula just given; for the odd interior points, the corresponding numbers are (0,0), (1,0), (2,1), (20,16) from accurate computations and (0.3,−0.3), (1.07,−0.07), (2.4,0.6), (21.6,14.4) from the approximate formula.

### (g) The grid quartic and self-similarity

In plots for real *k*, e.g. figure 1*a*, the interior grid points are seen to lie on smooth arcs. From equations (3.20) and (3.21), these arcs lie on the family of quartic curves
3.22
labelled by a positive integer *s* for which *s*=*m*_{1}+*m*_{2}=*n*_{1}+*n*_{2}−1. Each quartic curve has a real arc occupying the region
3.23
The arc is tangential to the *P*-wave line , i.e. *ω**h*/*c*_{1}=*k**h*, at the point (*ω*,*k*) for which
3.24
The left-hand part of an arc, i.e. the part to the left of the tangent point (3.24), forms a skeleton for the part of an exact branch which lies below the *P*-wave line. The right-hand part of an arc, i.e. the part to the right of the tangent point (3.24), does not correspond to a single branch, but crosses successive branches at grid points. The geometrical self-similarity of the arcs with different *s*, evident in figure 1*a*, is a consequence of the fact that in equation (3.22) the quantities *ω*, *k* and *s* occur only in the combinations *ω*/*s* and *k*/*s*. A suitable name for the curve (3.22) is the grid quartic. This quartic provides a broad-brush description of the solutions of the dispersion relation below the *P*-wave line *ω**h*/*c*_{1}=*k**h*. The fact that a quartic can provide such a description of the solutions of a transcendental equation is a further indication of the cancelling-out of Runge’s phenomenon.

In plots for imaginary *k*, e.g. figure 2*b*, the grid points lie on arcs obtained from equation (3.22) by putting *k*=i*k*_{i}, so that the left-hand side becomes −(*k*_{i}*h*/(2*π**s*))^{2}. There is then a bounded arc in the region *ω**h*/(2*π**s*)≤*c*_{1}*c*_{2}/(*c*_{1} + *c*_{2}) and an unbounded arc in *ω**h*/(2*π**s*)≥*c*_{1}*c*_{2}/(*c*_{1}−*c*_{2}). The grid points are at intersections of ellipses, not shown in the figure. Each unbounded arc forms a skeleton for an exact imaginary branch; self-similarity is apparent. On each unbounded arc, there is an inflection point (indicated by +), which explains the arc’s near straightness. The inflection point and tangent there may be calculated numerically by solving a cubic equation in (*ω**h*/(2*π**s*))^{2}; alternatively, a perturbation analysis in the small parameter *η*={(*c*_{1}−*c*_{2})/(*c*_{1}+*c*_{2})}^{2} gives the inflection point extremely accurately as
3.25
at which
3.26
On a plot, the inflection point and tangent line given by equations (3.25) and (3.26) are indistinguishable from their numerically computed values. By self-similarity, the inflection points for different arcs all lie on a straight line through the origin, indicated in the figure by the straight dotted line OI.

## 4. Accuracy analysis

The (*m*,*n*) finite-product approximation is exact at values of (*ω*,*k*) for which *α*=1. Hence a first attempt at estimating the accuracy of the approximation at arbitrary (*ω*,*k*) is to plot contours of *α* in the (*ω*,*k*) plane and expect that the region of accuracy is where *α* is close to 1. However, this approach is too pessimistic. Typical contours of *α* are shown in figure 3*a*, for the (1,2) finite-product approximation. The contours are roughly parallel to the *k* axis, and even where *α* is as low as 0.5 the (1,2) approximation is good. The same is true for any (*m*,*m*+1) approximation.

Another approach is to exploit the divisibility of equation (3.17) by , as in equations (3.3)–(3.6). This gives the exact dispersion relation in the form
4.1
where
4.2
Thus, the (*m*,*n*) approximation (3.5) is exact where *ε*=0, so that a possible accuracy criterion is |*ε*|≪1. This criterion is an improvement over the criterion based on *α*, as may be seen from the contours of *ε* shown in figure 3*b*. The shape of the contours correctly reveals the high accuracy of the approximation near the *Ω*-axis and also high up on the left-most branches. Nevertheless, the numerical values of |*ε*| are still pessimistically high. For example, in figure 3*b*, the region of the (*ω*,*k*) plane between the contours *ε*=−0.2 and *ε*=−0.4 contains segments of high accuracy which are not indicated by the values of *ε*.

The feature missing from accuracy criteria based on *α* or *ε* is that on the grid of points identified in §3*f* the approximate equations are exact, irrespective of the values of *α* or *ε*. Since the grid points are a distance ‘order-one’ apart in the dimensionless variables *Ω*=*ω**h*/*c*_{0} and *K*=*k**h*, it follows that an accuracy criterion should be based on the angle in the (*Ω*,*K*) plane between the tangents to the exact and approximate solutions at the grid points. An analytical advantage of this approach is that at the grid points many formulae for derivatives simplify by L’Hôpital’s rule, because the grid points are at zeros of the polynomials *S*_{m} and *C*_{n}. Let us denote solutions of the exact equations (3.17) or (4.1) by a subscript ‘ex’ and solutions of the corresponding finite-product approximations (in which we have put *α*=1 or *ε*=0) by ‘fp’. Then equation (3.17) gives for the even interior grid points
4.3
and for the odd interior grid points
4.4
The definitions of and in equation (2.9) give
4.5
so that
4.6
With the notation *K*_{Ω} for d*K*/d*Ω* on a dispersion curve, we obtain at even interior grid points
4.7
and at odd interior grid points
4.8
The difference between an approximate and an exact expression will be indicated by *δ*, so that the left-hand sides here may be written *δ*(*K*_{Ω})/*K*_{Ω}. To convert from slopes to angles in the (*Ω*,*K*) plane, we put and *δ**θ*=*θ*_{fp}−*θ*_{ex}. Hence (with *δ**θ* in radians)
4.9
The factor is tame: it never exceeds 1/2, a value attained when *K*_{Ω}=1, and it tends to zero as *K*_{Ω} when *K*_{Ω}→0 and as 1/*K*_{Ω} when . Hence equation (4.9), in conjunction with expressions (4.7) and (4.8), is a suitable tool for analysing the accuracy of a finite-product approximation. A check of formulae (4.7)–(4.9) is that they give *δ**θ*>0 at even interior grid points and *δ**θ*<0 at odd interior grid points, because *α*<1 at grid points. This is consistent with figure 2*a*,*c*,*e*,*g*, which shows that the approximate dispersion curves are steeper than the exact curves at the even points, but shallower than the exact curves at the odd points. The numerator of equation (4.7) gives *δ**θ*=0 on the *Ω*-axis (*K*=0), the Lamé line () and the *P*-wave line () and the numerator of equation (4.8) gives *δ**θ*=0 on the *Ω*-axis, the Lamé line and the *S*-wave line (). The occurrence of rather than in these numerators keeps *δ**θ* small in a broad region around the Lamé line. Similarly, the occurrence of *K*^{2} rather than *K* in the numerators keeps *δ**θ* small in a broad region around the *Ω*-axis.

Contours of *δ**θ* (in degrees) for (*m*,*n*)=(1,2) are shown in figure 3*c* for odd interior grid points; the corresponding contours for even interior grid points are similar. The contour values are *δ**θ*=0^{°}, −0.1^{°}, −1^{°}, −3^{°}, −5^{°}. The region of the (*Ω*,*K*) plane for which |*δ**θ*| is less than any specified value is apparent and can be seen to correspond to a region of closeness of the exact and approximate dispersion curves, thus confirming the validity of the accuracy criterion based on tangents at grid points. The large size of the regions for which |*δ**θ*|≤1^{°} or |*δ**θ*|≤3^{°}, for example, even for (*m*,*n*) as low as (1,2), demonstrates the high accuracy and wide range of validity of finite-product approximations. All formulae in this section apply whether *K* is wholly real or wholly imaginary. Thus analogues of figure 3 may be constructed for the curves in the (*Ω*,*K*_{i}) plane shown in figures 1*b* and 2*b*,*d*,*f*,*h*.

The uppermost curve in figure 3 is the branch connecting the Euler–Bernoulli bending-wave region, near the origin, to the Rayleigh-wave region, higher up. An accuracy criterion based on *δ**θ*, i.e. on tangents at grid points, cannot be expected to apply to the Rayleigh-wave region, because the only grid point on this branch is the Euler–Bernoulli point, at the origin. The contours of *α* plotted in figure 3*a* show that for the Rayleigh-wave region an appropriate error criterion is |*α*−1|≪1. Physically, this different criterion is required because a Rayleigh wave is entirely a boundary-layer phenomenon, whereas the waves corresponding to the high regions of all the other real branches are body shear waves perturbed only slightly by their boundary-layer adjustment region.

The finite-product approximations obtained in this paper contain no adjustable parameters: given any choice of *m* and *n*, all subsequent formulae are definite. Nevertheless, finite-product approximations can be ‘fine-tuned’, if desired, to give enhanced accuracy in chosen regions of the (*Ω*,*K*) plane, for example, to match exactly with the Rayleigh wave at high frequencies. If such fine tuning is not to destroy the agreement of the exact and approximate dispersion curves at grid points, the best approach is to replace the factor on the right-hand side of equation (3.5) by , so that the right-hand side of equation (3.5) becomes . Then *κ* may be chosen to give, for example, the correct Rayleigh-wave phase speed, or, for (*m*,*n*)=(0,1), the Timoshenko equation (2.12). The above replacement is not the same as replacing by on the right-hand side of equation (3.2); in the transition from equation (3.2) to equation (3.5), it is necessary to preserve the exact cancellation of the factor . Thus equation (3.5) is a more fundamental type of equation than equation (3.2).

## 5. Conclusions and further work

The results in this paper demonstrate that a finite-product polynomial approximation to a transcendental dispersion relation can be extraordinarily accurate, and, moreover, be amenable to an analytically based accuracy analysis founded on gamma functions. Thus, a finite-product method is a useful addition to directly numerical approaches (e.g. Pagneux & Maurel 2001) or Padé-approximant methods (e.g. Poncelet *et al.* 2006). An accurate polynomial approximation is particularly useful when the dispersion equation needs to be solved for a large number of complex roots; polynomial root finders are now trivial to use, whereas the task of computing all the complex roots of an arbitrary equation is still far from automated. The authors have verified, both for equation (2.7) and for more complicated dispersion relations, that the finite-product method finds all the complex branches of a dispersion relation just as accurately as it finds the real and imaginary branches; not a single spurious branch, of any type, is generated in this process, and the resulting three-dimensional line plots give a complete picture of all connections and bifurcations in three-dimensional space.

Many extensions of the results presented in this paper are possible, for example to other types of mode, to initial-value problems (Kaplunov *et al.* 2006) and to resonance (Pagneux 2006; Zernov *et al.* 2006). Other applications are to layers which may be fluid-loaded (Smith 2007; Sorokin 2007) or pre-stressed, composite or homogenized (Pichugin *et al.* 2008), anisotropic (Shuvalov & Poncelet 2008) or unstable or a selection of these (Fu & Rogerson 1994; Rogerson & Sandiford 1997; Sorokin 2004; Shuvalov *et al.* 2006). These papers contain a very wide range of dispersion relations which, however complicated, can be written in a form containing linear combinations of ratios of tangents, and so are natural generalizations of equation (2.7). The finite-product method applies to cylindrical geometry, with or without the complications just listed, because an infinite-product representation of Bessel functions is available (Abramowitz & Stegun 1965, p. 370), involving Bessel-function zeros in the factors. Theoretical developments are also possible, especially in relation to wave hierarchies (Whitham 1974, p. 353). The (0,1) finite-product approximation is explicitly of Whitham’s type, and satisfies the criterion Whitham gives for stability. Hence the question arises of whether equation (3.2) or (3.5), perhaps with some restriction to allowed values of *m* and *n*, can be embedded within a Scott-type hierarchical theory (e.g. Scott 2008).

## Acknowledgements

This work has been supported by a Joint Project Grant from the Royal Society. The authors thank J. D. Kaplunov and A. V. Pichugin for helpful remarks.

## Footnotes

- Received May 11, 2009.
- Accepted September 23, 2009.

- © 2009 The Royal Society