Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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.

share|improve this question
    
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. –  J. M. 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. –  J. M. Aug 13 '11 at 7:04
    
The plot really helps! –  Alice Aug 18 '11 at 21:08
add comment

4 Answers

[ 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)]);

enter image description here

share|improve this answer
    
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
add comment

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.

share|improve this answer
    
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
add comment
up vote 1 down vote accepted

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.

share|improve this answer
    
I'm still interested in any other methods someone may come up with to obtain these curves. –  becko Aug 20 '11 at 6:39
add comment

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<AGROW>
            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

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.