Numerical computation of the Rayleigh-Lamb curves

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$):

-
once you insert $\omega$ you have an equation $f(k)=0$ that you want to solve for $k$. This is a 1 d non-linear root finding problem (as you can see from reading a few introductory paragraphs from each section in the book numerical recipes (I would not use their code though). As it is 1d you can bracket the root to search for it using secant or other method. – Alice Aug 12 '11 at 11:44
can you solve for the eigenvalues analytically? and then just use the software to plot the graph? : ) – Timtam Aug 13 '11 at 12:39
@Timtam: Solving for the eigenvalues analytically is part of the problem – becko Aug 18 '11 at 2:06
@Alice: The problem is there is more than one root. Blindly bracketing for roots can miss some of them. – becko Aug 18 '11 at 2:06
@becko. Only three free parameters, $d$, $\omega/c_T$, $\omega/c_L$ that will specify the roots. For most situations $d$ likely tells you the spacing with roots every $k^2 d$ mod $2 \pi$. Once you bracket the roots you can run a root finder. If you plot the function for various regimes of your three parameters you can check that you have appropriately bracketed the roots. – Alice Aug 18 '11 at 12:49

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.

-
+1 That's a good idea certainly worth exploring. I'll keep you posted. – becko Aug 19 '11 at 16:22
I am glad it is along the lines of what you wanted. For a lightweight problem like this, brute force Monte Carlo approach may pay off. The calculation will not be that long, especially if you parallelize it, but you will not need to put a lot of time and effort into the search for a more elegant solution. – drlemon Aug 19 '11 at 16:41