Computación numérica de la de Rayleigh-Cordero curvas

Question

La de Rayleigh-Cordero ecuaciones:

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

(dos ecuaciones, una de ellas con el +1 exponente y el otro con el exponente -1) donde

$$p^2=\frac{\omega ^2}{c_L^2}-k^2$$ y $$q^2=\frac{\omega ^2}{c_T^2}-k^2$$

se muestran en las consideraciones físicas de la goma oscilaciones de placas sólidas. Aquí, $c_L$, $c_T$ y $d$ son constantes positivas. Estas ecuaciones determinan para cada valor positivo de la $\omega$ un conjunto discreto de real "valores propios" de $k$. Mi problema es el numérico cálculo de estos valores y, en particular, para la obtención de las curvas de la visualización de estos autovalores. ¿Qué tipo de método numérico puedo utilizar con este problema? Gracias.

Edit: Usando los valores numéricos $d=1$, $c_L=1.98$, $c_T=1$, las parcelas debe ser algo como esto (negro curvas corresponden a la exponente -1, azul curvas para el +1 exponente; el eje horizontal es $\omega$ y el eje vertical es $k$):

Answer 1

[ EDIT: se incluyen tanto $+$ $-$ curvas, intercambiar $k$ $\omega$ ejes como por su imagen ]

Aquí es la trama de $d=1$, $c_L = 1.98$, $c_T = 1$, $\omega$ de 0 a 14. Tenga en cuenta que tenemos $\omega \ge c_L k$ $q$ a ser real, así que me tomé $k$$14/c_L$. El Arce de comandos:

nca:= 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)}); nca:= eval(nca,{d=1; cl=1.98,ct=1}); con(parcelas): cols:= [azul,negro]: display([seq(implicitplot(nca[i],omega=0..14, k= 0 .. omega/1.98 - .01, grid=[50,50], gridrefine=3, crossingrefine=3, signchange=false, color=cols[i]),i=1..2)]);

Answer 2

Los métodos estándar para resolver numéricamente no-ecuaciones polinómicas debe trabajar para encontrar $k$ en un intervalo de tiempo dado para un valor dado de a $\omega$. En Maple me gustaría utilizar el fsolve comando para que. Para trazar las soluciones de los intervalos de $\omega$ $k$ me gustaría utilizar el implicitplot comando.

Answer 3

La de Rayleigh-Cordero ecuaciones:

$$\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 $p$ became imaginary. The trick is to divide by $p$ or $p$ 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 $p$. El suministro de estos a la ContourPlot comando en Mathematica he obtenido las curvas

para el +1 exponente, y

para el exponente -1.

Answer 4

Hay una descripción del algoritmo en el libro de las Ondas Ultrasónicas en Medios Sólidos por Joseph L. Rosa:

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

Y aquí están mis dos implementaciones en MATLAB:

1. El uso de la trama implícita:

   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. El uso de fzero función:

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

Donde la función 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) = [];

Otra idea es la trama de superficie 3D:

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