Numerical computation of the Rayleigh-Lamb curves

Question

The Rayleigh-Lamb equations:

$$\frac{\tan (pd)}{\tan (qd)}=-\left[\frac{4k^2pq}{\left(k^2-q^2\right)^2}\right]^{\pm 1}$$

(two equations, one with the +1 exponent and the other with the -1 exponent) where

$$p^2=\frac{\omega ^2}{c_L^2}-k^2$$ and $$q^2=\frac{\omega ^2}{c_T^2}-k^2$$

show up in physical considerations of the elastic oscillations of solid plates. Here, $d$ is the thickness of a elastic plate, $c_L$ the velocity of longitudinal waves, and $c_T$ the velocity of transverse waves. These equations determine for each positive value of $\omega$ a discrete set of real "eigenvalues" for $k$. My problem is the numerical computation of these eigenvalues and, in particular, to obtain curves displaying these eigenvalues. What sort of numerical method can I use with this problem? Thanks.

Edit: Using the numerical values $d=1$, $c_L=1.98$, $c_T=1$, the plots should look something like this (black curves correspond to the -1 exponent, blue curves to the +1 exponent; the horizontal axis is $\omega$ and the vertical axis is $k$):

Answer 1

How about instead of finding the roots and then making the plots, you skip right to the plots using a Monte Carlo method?

Choose a random k, then choose a random ω and calculate the left-hand side and the right-hand side of the R-L equations. If the RHS is close enough to the LHS (you pick how close), put a point on the plot (blue or black, depending on which branch you used).

The more points you process, and the more stringent the condition of RHS=LHS equality you pick, the more accurate will the plot look.

A problem with this approach is that when you put a point on the plot, your algorithm doesn't know which branch it belongs to. But if it is the plot you are after, you will have no problem figuring it out by the eye when the calculation is complete.

To numerically read from that Monte Carlo plot, you can sort the solutions you find (sort them in k-space or in ω-space), and do some kind of search using interpolation.

Here is an approximate C code to explain what I mean:

void main()
{
srand(1);
const int N=1000000000;
const float eps=0.01;
const float cl=...;
const float ct=...;
const float kmax=20.0;
const float omegamax=20.0;
float k, omega, p, q, lhs, rhs;
for (int i=0; i<N; i++)
 {
  k=(float)rand()/(float)RAND_MAX*k_max;
  omega=(float)rand()/(float)RAND_MAX*omega_max;
  lhs=tan(...)/tan()...;
  rhs=(4.0*k*k...);
  if (fabs(lhs-rhs) < eps*fabs(lhs+rhs))
    printf("%f %f\n", k, omega);
 }
}

and that's it.