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
    I posted a somewhat related question at http://physics.stackexchange.com/questions/10849/rayleigh-lamb-dispersion-curves, just in case you might want a bit more info on the physics of the problem.2011-08-05
  • 0
    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.2011-08-07
  • 0
    @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$.2011-08-09
  • 0
    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.2011-08-13
  • 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
    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.2011-08-18
  • 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
    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.2011-08-06
  • 0
    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.2011-08-07
  • 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