5
$\begingroup$

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

The curves are supposed to look like this.

  • 0
    The plot really helps!2011-08-18

4 Answers 4

1

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

http://books.google.de/books?id=DEtHDJJ-RS4C&pg=PA110&dq=&redir_esc=y#v=onepage&q&f=false

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') 

enter image description here

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

enter image description here

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)     shading interp     view(2)     caxis([-5e-10 5e-10]) 

enter image description here

4

[ 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(pd)/tan(qd) + 4*k^2*pq/(k^2-q^2)^2, tan(pd)/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)]);

enter image description here

  • 0
    These Maple curves don't look very accurate either. I'm writting an answer now using a suggestion that I got at http://stackoverflow.com/questions/7088632/how-to-obtain-accurate-plot-curves-in-mathematica2011-08-20
2

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.

  • 0
    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.2011-08-13
1

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

Curves for the +1 exponent

for the +1 exponent, and

Curves for the -1 exponent

for the -1 exponent.

  • 0
    I'm still interested in any other methods someone may come up with to obtain these curves.2011-08-20