# 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, $c_L$, $c_T$ and $d$ are positive constants. 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$):

-
I posted a somewhat related question at physics.stackexchange.com/questions/10849/…;, just in case you might want a bit more info on the physics of the problem. –  becko Aug 5 '11 at 20:53
To be explicit... $k$ is the only unknown, or both $\omega$ and $k$? Newton-Raphson is standard for solving (simultaneous) transcendental equations; your problem here is coming up with initial approximations that Newton-Raphson can polish to a usable answer. –  Ｊ. Ｍ. Aug 7 '11 at 5:53
@J.M.: $k$ is the only unknown, in the sense that $\omega$ goes from 0 to infinity and, for each value of $\omega$, my objective is to find the corresponding values of $k$. Hope it's clearer now. And yes, if I use a Newton-Rhapson method then the problem reduces to finding first good approximations to $k$. –  becko Aug 9 '11 at 1:45
I'm working on it; but here's an idea I'm tossing out in case somebody does a better job than me: try replacing the tangent functions in the equation(s) with approximations of it (e.g. series or the Padé approximants); you'll obtain a polynomial equation that is hopefully slightly more tractable to solve, and then use those solutions to kickstart Newton-Raphson or some other iterative scheme. –  Ｊ. Ｍ. Aug 13 '11 at 7:04
The plot really helps! –  Alice Aug 18 '11 at 21:08

[ EDIT: included both $+$ and $-$ curves, interchanged $k$ and $\omega$ axes as per your image ]

Here is the plot for $d=1$, $c_L = 1.98$, $c_T = 1$, $\omega$ from 0 to 14. Note that we need $\omega \ge c_L k$ for $q$ to be real, so I took $k$ up to $14/c_L$. The Maple commands were:

eqs:= eval([tan(p*d)/tan(q*d) + 4*k^2*p*q/(k^2-q^2)^2, tan(p*d)/tan(q*d) + (4*k^2*p*q/(k^2-q^2)^2)^(-1)], {p=sqrt(omega^2/cl^2-k^2), q=sqrt(omega^2/ct^2 - k^2)}); eqs:= eval(eqs,{d=1,cl=1.98,ct=1}); with(plots): cols:= [blue,black]: display([seq(implicitplot(eqs[i],omega=0..14, k= 0 .. omega/1.98 - .01, grid=[50,50], gridrefine=3, crossingrefine=3, signchange=false, colour=cols[i]),i=1..2)]);

-
I edited my question to add an image of what the plot's supposed to look like. I'm looking at what you posted now. –  becko Aug 18 '11 at 5:35
These Maple curves don't look very accurate either. I'm writting an answer now using a suggestion that I got at stackoverflow.com/questions/7088632/… –  becko Aug 20 '11 at 6:19

The standard methods for numerically solving non-polynomial equations should work to find $k$ in a given interval for a given value of $\omega$. In Maple I would use the fsolve command for that. To plot the solutions given intervals for $\omega$ and $k$ I would use the implicitplot command.

-
I should have clarified that I'm trying to understand the numerical method used, and perhaps even program it myself. That being said, I have tried to use Mathematica, in particular the ContourPlot command, which I pressume is the equivalent of implicitplot in Maple. The curves I have obtained are very inaccurate. –  becko Aug 6 '11 at 0:47
It may help to write tan as sin/cos, and take the numerator of the difference of the two sides when placed over a common denominator, thus avoiding trouble with singularities. Note also that you probably want to restrict your attention to the region where $p$ and $q$ are real. If you provide typical values for the parameters, I can show you what Maple's plot looks like. –  Robert Israel Aug 7 '11 at 18:38
The following are typical numerical values for the constants: $d=1$, $c_L=1.98$ and $c_T=1$. Here, $k$ could go from $0$ to $14$ and $\omega$ from $0$ to $6$. And yes, the interesting regions are those where and are real. In Mathematica the curves I obtain have jumps and irregularities and even some zeros are missed completely. Let's see what the Maple plots look like. –  becko Aug 13 '11 at 15:12

The Rayleigh-Lamb equations:

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

are equivalent to the equations (as Robert Israel pointed out in a comment above)

$$\left(k^2-q^2\right)^2\sin pd \cos qd+4k^2pq \cos pd \sin qd=0$$

when the exponent is +1, and

$$\left(k^2-q^2\right)^2\cos pd \sin qd+4k^2pq \sin pd \cos qd=0$$

when the exponent is -1. Mathematica had trouble with the plots because $p$ or $q$ became imaginary. The trick is to divide by $p$ or $q$ in convenience.

Using the numerical values $d=1$, $c_L=1.98$, $c_T=1$, we divide equation for the +1 exponent by $p$ and the equation for the -1 exponent by $q$. Supplying these to the ContourPlot command in Mathematica I obtained the curves

for the +1 exponent, and

for the -1 exponent.

-
I'm still interested in any other methods someone may come up with to obtain these curves. –  becko Aug 20 '11 at 6:39

There is a description of algorithm in the book Ultrasonic Waves in Solid Media by Joseph L. Rose:

And here are my two implementations in MATLAB:

1. Using implicit plot:

close all
clear all
clc
cL = 1.98;
cT = 1;
d = 1;
p = @(k,w)sqrt(w.^2/cL^2-k.^2);
q = @(k,w)sqrt(w.^2/cT^2-k.^2);
symmetric = @(w,k)tan(q(k,w)*d)./q(k,w)+4*k.^2.*p(k,w).*tan(p(k,w)*d)./(q(k,w).^2-k.^2).^2;
asymmetric = @(w,k)q(k,w).*tan(q(k,w)*d)+(q(k,w).^2-k.^2).^2.*tan(p(k,w)*d)./(4*p(k,w).*k.^2);
h1 = ezplot(symmetric,[0 6 0 14]);
hold on
h2 = ezplot(asymmetric,[0 6 0 14]);
set(h1, 'Color', 'b');
set(h2, 'Color', 'k')
legend('symmetric','antisymmetric')
set(gca,'YGrid','on')


1. Using fzero function:

% function tries to find all roots of disspersion equations for Lamb waves
close all
clear all
clc
cL = 1.98;
cT = 1;
d = 1;
N = 1000;
omega = linspace(0,6,N);

for idx = N:-1:1
w = omega(idx);
p = @(k)sqrt(w.^2/cL^2-k.^2);
q = @(k)sqrt(w.^2/cT^2-k.^2);
symmetric = @(k)tan(q(k)*d)./q(k)+4*k.^2.*p(k).*tan(p(k)*d)./(q(k).^2-k.^2).^2;
asymmetric = @(k)q(k).*tan(q(k)*d)+(q(k).^2-k.^2).^2.*tan(p(k)*d)./(4*p(k).*k.^2);
try
lb = 0;
ub = 14;
bstep = 0.1;
tmps = findAllZeros(symmetric,lb,ub,bstep);
tmpa = findAllZeros(asymmetric,lb,ub,bstep);
result{idx,1} = [w tmps];
result{idx,2} = [w tmpa];
catch ME
disp(ME)
end
end

%%
figure
hold on
h1 = plot(NaN,NaN,'b.');
h2 = plot(NaN,NaN,'k.');
hold off
for idx=1:N
if numel(result{idx,1}) > 1
x1 = result{idx,1}(1);
y1 = result{idx,1}(2:end);
x1 = [x1+0*y1 get(h1,'xdata')];
y1 = [y1 get(h1,'ydata')];
set(h1,'xdata',x1,'ydata',y1)
end
if numel(result{idx,2}) > 1
x2 = result{idx,2}(1);
y2 = result{idx,2}(2:end);
x2 = [x2+0*y2 get(h2,'xdata')];
y2 = [y2 get(h2,'ydata')];
set(h2,'xdata',x2,'ydata',y2)
drawnow
end
end
xlim([0 6])
ylim([0 14])
set(gca,'YGrid','on')
xlabel('\omega')
ylabel('k')


Where function findAllZeros:

    function x = findAllZeros(fun,lb,ub,bstep)
% fun - handle to the function
% lb - a lower bound for the function.
% ub - an upper bound for the function.
% bstep - step for iteration

x = []; % Initializes x.

for i=lb:bstep:ub
if sign(fun(i-bstep))~=sign(fun(i+bstep))
tmp = fzero(fun, i);
if isreal(tmp) && abs(fun(tmp))<1 % eliminate complex values and discontinuities
x = [x tmp]; %#ok<AGROW>
end
end
end

% Make sure that there are no duplicates.
x = unique(x);
DUPE = (diff([x NaN]) < 1e-16) | isnan(x);
x(DUPE) = [];


Another idea is to plot 3D surface:

    close all
clear all
clc
figure
cL = 1.98;
cT = 1;
d = 1;
p = @(k,w)sqrt(w.^2/cL^2-k.^2);
q = @(k,w)sqrt(w.^2/cT^2-k.^2);
symmetric = @(w,k)tan(q(k,w)*d)./q(k,w)+4*k.^2.*p(k,w).*tan(p(k,w)*d)./(q(k,w).^2-k.^2).^2;
N = 1000;
[ww,kk] = meshgrid(linspace(0,6,N),linspace(0,14,N));
zz = symmetric(ww,kk);
surf(ww,kk,zz)