Modos de vibración de una cuerda no homogénea

Cuerda formada por dos segmentos

Las ecuaciones diferenciales del movimiento ondulatorio, son respectivamente

v 1 2 2 ψ 1 x 2 = 2 ψ 1 t 2 v 2 2 2 ψ 2 x 2 = 2 ψ 2 t 2   

Ψ es el desplazamiento trasversal de un punto x de la cuerda en el instante t.

v1 y v2 son las velocidades de propagación de las ondas trasversales en las partes izquierda y derecha de la cuerda sometida a la tensión T

v 1 = T μ 1 v 2 = T μ 2

Siguiendo los mismos pasos que para una cuerda homogénea sujeta por ambos extremos, ensayamos la solución

Ψ(x,t)=y(x)sin(ωt)

que describe una onda estacionaria. Cada punto de la cuerda vibra con una amplitud y(x) y con frecuencia angular ω.

Introduciendo esta expresión en cada una de las dos ecuaciones diferenciales del movimiento ondulatorio, obtenemos las ecuaciones diferenciales similares a las de un M.A.S.

d 2 y 1 d x 2 + ω 2 v 1 2 y 1 =0 d 2 y 2 d x 2 + ω 2 v 2 2 y 2 =0

Las soluciones de estas ecuaciones diferenciales son

y1(x)=A1sin (k1x)+B1cos(k1x
y2
(x)=A2sin (k2x)+B2cos(k2x

k1=ω/v1 y k2=ω/v2  son los números de onda, respectivos

Los coeficientes A1, B1, A2 y B2, vienen determinados por la continuidad de la función de onda en x=0 y por las condiciones de contorno en los extremos fijos x=0 y x=L.

Condiciones de contorno

El extremo izquierdo de la cuerda está sujeto en x=0 y el extremo derecho en x=L. Se tiene que cumplir que

{ y 1 (0)=0, B 1 =0 y 2 (L)=0, A 2 sin( k 2 L)+ B 2 cos( k 2 L)=0

Continuidad de la función que describe la onda estacionaria

En el punto x=a, las dos mitades de la cuerda se unen. La función que describe la onda estacionaria tiene que ser continua y también su derivada primera.

y 1 (a)= y 2 (a) A 1 sin( k 1 a )= A 2 sin( k 2 a )+ B 2 cos( k 2 a ) y 1 x | a = y 2 x | a k 1 A 1 sin( k 1 a )= k 2 ( A 2 cos( k 2 a ) B 2 sin( k 2 a ) )

Tenemos un sistema homogéneo de tres ecuaciones con tres incógnitas, A1, A2 y B2, despejamos B2 en la primera ecuación

B 2 = A 2 sin( k 2 L) cos( k 2 L)

Sustituimos en las otras dos

{ A 1 sin( k 1 a )= A 2 sin( k 2 a ) A 2 sin( k 2 L) cos( k 2 L) cos( k 2 a ) k 1 A 1 sin( k 1 a )= k 2 ( A 2 cos( k 2 a )+ A 2 sin( k 2 L) cos( k 2 L) sin( k 2 a ) )

Utilizamos la fórmula, cos(A-B)=cosA·cosB+sinA·sinB, dividimos las dos ecuaciones

1 k 1 sin( k 1 a ) cos( k 1 a ) = 1 k 2 sin( k 2 (aL) ) cos( k 2 (aL) )

Frecuencias de los modos de vibración de la cuerda

Para calcular las frecuencias de los modos de vibración de la cuerda se resuelve la ecuación trascendente en ω

k 2 sin( k 1 a )cos( k 2 (La) )+ k 1 cos( k 1 a )sin( k 2 (La) )=0

Llamemos r al cociente de velocidades de propagación en las porciones de cuerda, v2=r·v1 o bien, k2=k1/r

k 1 r sin( k 1 a )cos( k 1 r (La) )+ k 1 cos( k 1 a )sin( k 1 r (La) )=0

Sea la longitud de la cuerda L=1 m, la proción izquierda de la cuerda tiene una longitud a=L/3. Las densidades lineales de la cuerda son tales que la relación de las velocidades de propagación entre las dos porciones de la cuerda es, v2=2v1

Como veremos en el ejemplo, más abajo, en este caso, no es necesario resolver la ecuación transcendente por procedimientos numéricos.

function cuerda3
    L=1; %longitud de la cuerda
    a=L/3;  %porción izquierda
    r=2; %v2=r·v1, k2=k1/r
    f=@(x) x.*sin(a*x).*cos((L-a)*x/r)/r+x.*cos(a*x).*sin(x*(L-a)/r);
    fplot(f,[0,20]);
    grid on
    xlabel('x')
    ylabel('f(x)')
    title('Raíces')
    
    x=linspace(1,20,20);
    disp('frecuencias (Hz)')
    raiz=raices(f,x);
    disp(raiz/(2*pi))

    function r = raices(f, x)
        y=f(x);
        indices=find(y(1:end-1).*y(2:end)<0);
        r=zeros(1,length(indices));
        for k=1:length(indices)
            r(k)=fzero(f, [x(indices(k)), x(indices(k)+1)]);
        end
    end
end

Las frecuencias ω=v1k1, f=ω/2π de los primeros modos de vibración son

frecuencias (Hz)
    0.7500    1.5000    2.2500    3.0000

Amplitudes de los modos de vibración

Despejamos A2 y B2 en función de A1, tomando las dos primeras ecuaciones del sistema homogéneo de tres ecuaciones con tres incógnitas

A 2 = A 1 cos( k 2 L) sin( k 1 a) sin( k 2 (La) ) B 2 = A 2 sin( k 2 L) cos( k 2 L) = A 1 sin( k 2 L) sin( k 1 a) sin( k 2 (La) )

Puede ocurrir una indeterminación 0/0, cuando k1a=nπ y a la vez, k2(L-a)=mπ. La indeterminación se resuelve aplicando la regla de de L'Hôpital tal como veremos en el ejemplo.

A 2 = A 1 cos( k 2 L) lim xa sin( k 1 x) sin( k 2 (xL)) = A 1 cos( k 2 L) k 1 k 2 cos( k 1 a) cos( k 2 (La) ) B 2 = A 1 sin( k 2 L) k 1 k 2 cos( k 1 a) cos( k 2 (La) )

Para el modo de vibración de frecuencia angular ω=f, con k1=ω/v1 y k2=ω/v2

{ y 1 (x)= A 1 sin( k 1 x ),0x<a y 2 (x)= A 2 sin( k 2 x )+ B 2 cos( k 2 x ),ax<L

A1 es un factor de escala. Para que todos los modos de vibración se representen a la misma escala, se calcula A1 haciendo que el área bajo la curva del cuadrado de la amplitud y(x) multiplicada por la densidad en cada tramo sea la unidad

μ 1 0 a y 1 2 (x)dx + μ 2 a L y 2 2 (x)dx =1 μ 1 0 a A 1 2 sin 2 ( k 1 x)dx + μ 2 a L ( A 2 sin( k 2 x)+ B 2 cos( k 2 x) ) 2 dx =1

La integral definida de la izquierda es el área sombreada en color azul y la integral de la derecha es el área sombreada en color rojo.

Relaciones de ortogonalidad

Las relaciones de ortogonalidad de Φn(x) para una cuerda no homogénea de densidad μ(x)

0 L μ(x) Φ m (x) Φ n (x)dx =0,mn μ 1 0 a y 1m (x) y 1n (x)dx + μ 2 a L y 2m (x) y 2n (x)dx=0, mn

m y n son dos los modos de vibración distintos. Véase el ejemplo más abajo

Ejemplo

Sea una cuerda no homogénea, de longitud L=1 m. La parte izquierda a=L/3, de densidad μ1=4μ, la parte derecha de longitud L-a=2L/3 de densidad lineal μ2=μ

La velocidad de propagación de las ondas trasnversales en las cuerdas son

v 1 = T μ 1 = T 4μ v 2 = T μ 2 = T μ v 2 =2 v 1 k 1 = ω v 1 , k 2 = ω v 2 = k 1 2

Frecuencias de los modos de vibración

k 1 2 sin( k 1 L 3 )cos( k 1 2 2L 3 )+ k 1 cos( k 1 L 3 )sin( k 1 2 2L 3 )=0 3 2 k 1 sin( k 1 L 3 )cos( k 1 L 3 )=0 3 4 k 1 sin( 2 k 1 L 3 )=0 k 1 L= 3nπ 2 ω n =nπ 3 v 1 2L ,n=1,2,3...

Amplitudes

Comprobamos las relaciones de ortogonalidad, sea m=1 y n=2

L=1;
a=L/3;
mu=4; %densidad parte izquierda, mu=1, derecha
k1=3*pi/(2*L);  %modo de vibración n=1
k2=k1/sqrt(mu);
a2=-cos(k2*L)*sin(k1*a)/sin(k2*(L-a));
b2=sin(k2*L)*sin(k1*a)/sin(k2*(L-a));
y1=@(x) sin(k1*x);
y2=@(x) a2*sin(k2*x)+b2*cos(k2*x);

K1=6*pi/(2*L);  %modo de vibración n=2; 
K2=K1/2;
Y1=@(x) sin(K1*x);
A2=cos(K2*L)*K1*cos(K1*a)/(K2*cos(K2*(L-a)));
B2=-sin(K2*L)*K1*cos(K1*a)/(K2*cos(K2*(L-a)));
Y2=@(x) A2*sin(K2*x)+B2*cos(K2*x);

z1=@(x) y1(x).*Y1(x);
z2=@(x) y2(x).*Y2(x);
res=mu*integral(z1,0,a)+integral(z2,a,L);
disp(res)

Comprobamos las relaciones de ortogonalidad, la integral vale cero. Los coeficientes son los que hemos obtenido con n=1 y n=2

    0
>> a2,b2
a2 =    0.7071
b2 =    0.7071
>> A2,B2
A2 =  -3.6739e-16
B2 =     2 

Representamos los primeros modos de vibración

L=1; %longitud de la cuerda
a=L/3; %porción izquierda
mu=4; %densidad parte izquierda, mu=1, derecha
hold on
for n=1:4
    k1=3*n*pi/(2*L); 
    k2=k1/sqrt(mu); %Relación de velocidades de propagación, v2=2*v1
    A2=-sin(k1*a)*cos(k2*L)/sin(k2*(L-a));
    B2=sin(k2*L)*sin(k1*a)/sin(k2*(L-a));
    if mod(k1*a/pi,1)<1e-6 && mod(k2*(L-a)/pi,1)<1e-6 %indeterminación 0/0
        A2=k1*cos(k2*L)*cos(k1*a)/(k2*cos(k2*(L-a)));
        B2=-k1*sin(k2*L)*cos(k1*a)/(k2*cos(k2*(L-a)));
    end
    y1=@(x) sin(k1*x);
    y2=@(x) A2*sin(k2*x)+B2*cos(k2*x);
    z1=@(x) y1(x).^2;
    z2=@(x) y2(x).^2;
    area=mu*integral(z1,0,a)+integral(z2,a,L);
    A1=1/sqrt(area);
    f=@(x) A1*(x<a).*y1(x)+(x>=a).*y2(x)*A1;
    fplot(f,[0,L],'displayName',num2str(k1/(2*pi)))
end
hold off
grid on
legend('-DynamicLegend','location','southwest')
xlabel('x')
ylabel('y')
title('Modos de vibración')

Cambiamos el valor de a=L/2. Resolvemos la ecuación transcendente para calcular las frecuencias de los primeros modos de vibración.

function cuerda4
    L=1; %longitud de la cuerda
    a=L/2; %porción izquierda
    mu=4; %densidad parte izquierda, mu=1, derecha
    r=sqrt(mu); %v2=r·v1, k2=k1/r
    %raices, k1
    f=@(x) x.*sin(a*x).*cos((L-a)*x/r)/r+x.*cos(a*x).*sin(x*(L-a)/r);
    x=linspace(1,15,20);
    raiz=raices(f,x);
    
    hold on
    for n=1:length(raiz)
        k1=raiz(n);
        k2=k1/r;
        A2=sin(k1*a)*cos(k2*L)/sin(k2*(a-L));
        B2=-sin(k2*L)*sin(k1*a)/sin(k2*(a-L));
        if mod(k1*a/pi,1)<1e-6 && mod(k2*(L-a)/pi,1)<1e-6 %indeterminación 0/0
            A2=k1*cos(k2*L)*cos(k1*a)/(k2*cos(k2*(a-L)));
            B2=-k1*sin(k2*L)*cos(k1*a)/(k2*cos(k2*(a-L)));
        end
        y1=@(x) sin(k1*x);
        y2=@(x) A2*sin(k2*x)+B2*cos(k2*x);
        z1=@(x) y1(x).^2;
        z2=@(x) y2(x).^2;
        area=mu*integral(z1,0,a)+integral(z2,a,L);
        A1=1/sqrt(area);
        f=@(x) A1*(x<a).*y1(x)+(x>=a).*y2(x)*A1;
        fplot(f,[0,L],'displayName',num2str(k1/(2*pi)))
    end
    hold off
    grid on
    legend('-DynamicLegend','location','southwest')
    xlabel('x')
    ylabel('y')
    title('Modos de vibración')

    function r = raices(f, x)
        y=f(x);
        indices=find(y(1:end-1).*y(2:end)<0);
        r=zeros(1,length(indices));
        for k=1:length(indices)
            r(k)=fzero(f, [x(indices(k)), x(indices(k)+1)]);
        end
    end
end

Comprobamos la relación de ortogonalidad con n=1 y m=2.

L=1; %longitud de la cuerda
a=L/2; %porción izquierda
mu=4; %densidad parte izquierda, mu=1, derecha
r=sqrt(mu); %v2=r·v1, k2=k1/r
%raices, k1
f=@(x) x.*sin(a*x).*cos((L-a)*x/r)/r+x.*cos(a*x).*sin(x*(L-a)/r);
x=linspace(1,40,20);
raiz=raices(f,x);

k1=raiz(1);  %modo de vibración n=1
k2=k1/r;
a2=-cos(k2*L)*sin(k1*a)/sin(k2*(L-a));
b2=sin(k2*L)*sin(k1*a)/sin(k2*(L-a));
y1=@(x) sin(k1*x);
y2=@(x) a2*sin(k2*x)+b2*cos(k2*x);

K1=raiz(2);  %modo de vibración n=2; 
K2=K1/r;
Y1=@(x) sin(K1*x);
A2=-cos(K2*L)*sin(K1*a)/sin(K2*(L-a));
B2=sin(K2*L)*sin(K1*a)/sin(K2*(L-a));

%A2=cos(K2*L)*K1*cos(K1*a)/(K2*cos(K2*(L-a))); %para n=3, 6, ...
%B2=-sin(K2*L)*K1*cos(K1*a)/(K2*cos(K2*(L-a)));
Y2=@(x) A2*sin(K2*x)+B2*cos(K2*x);

z1=@(x) y1(x).*Y1(x);
z2=@(x) y2(x).*Y2(x);
res=mu*integral(z1,0,a)+integral(z2,a,L);
disp(res)
   1.1102e-16

En este caso, a=L/2, para n=3, 6, 9... hay que utilizar las amplitudes obtenidas aplicando la regla de L'Hôpital (para resolver la indeterminación 0/0)

Cuerda formada por tres segmentos

Estudiamos los modos de vibración producidas en una cuerda no homogénea sujeta por ambos extremos. La cuerda de longitud 2L está formada por tres partes una de longitud L-a, otra de longitud 2a y otra de longitud L-a de distinta densidad lineal μ1, μ2, y μ3, respectivamente. Las ecuaciones diferenciales del movimiento ondulatorio, son respectivamente

v 1 2 2 ψ 1 x 2 = 2 ψ 1 t 2 , v 2 2 2 ψ 2 x 2 = 2 ψ 2 t 2  , v 3 2 2 ψ 3 x 2 = 2 ψ 3 t 2

v1, v2 y v3 son las velocidades de propagación de las ondas trasversales en las tres partes de la cuerda cuya tensión es T

v 1 = T μ 1 , v 2 = T μ 2 , v 3 = T μ 3

Siguiendo los mismos pasos que para una cuerda homogénea sujeta por ambos extremos, ensayamos la solución

Ψ(x,t)=y(x)sin(ωt)

que describe una onda estacionaria. Cada punto de la cuerda vibra con una amplitud y(x) y con frecuencia angular ω.

Introduciendo esta expresión en cada una de las tres ecuaciones diferenciales del movimiento ondulatorio, obtenemos las ecuaciones diferenciales similares a las de un M.A.S.

d 2 y 1 d x 2 + ω 2 v 1 2 y 1 =0, d 2 y 2 d x 2 + ω 2 v 2 2 y 2 =0, d 2 y 3 d x 2 + ω 2 v 2 2 y 3 =0

Las soluciones de estas ecuaciones diferenciales son la parte real de

{ y 1 (x)= A 1 exp( i k 1 x )+ B 1 exp( i k 1 x ), k 1 = ω v 1 y 2 (x)= A 2 exp( i k 2 x )+ B 2 exp( i k 2 x ), k 2 = ω v 2 y 3 (x)= A 3 exp( i k 3 x )+ B 3 exp( i k 3 x ), k 3 = ω v 3

Los coeficientes A1, B1, A2, B2, A3 y B3 vienen determinados por la continuidad de la función de onda en x=-a, x=a y por las condiciones de contorno en los extremos fijos x=-L y x=L.

Continuidad de la función que describe la onda estacionaria

En x=-a y x=a, la función que describe la onda estacionaria tiene que ser continua y también su derivada primera.

x=a,{ y 1 (a)= y 1 (a) A 1 exp( i k 1 a )+ B 1 exp( i k 1 a )= A 2 exp( i k 2 a )+ B 2 exp( i k 2 a ) d y 1 dx | x=a = d y 2 dx | x=a i k 1 ( A 1 exp( i k 1 a ) B 1 exp( i k 1 a ) )=i k 2 ( A 2 exp( i k 2 a ) B 2 exp( i k 2 a ) ) x=a,{ y 2 (a)= y 3 (a) A 2 exp( i k 2 a )+ B 2 exp( i k 2 a )= A 3 exp( i k 3 a )+ B 3 exp( i k 3 a ) d y 2 dx | x=a = d y 3 dx | x=a i k 2 ( A 2 exp( i k 2 a ) B 2 exp( i k 2 a ) )=i k 3 ( A 3 exp( i k 3 a ) B 3 exp( i k 3 a ) )

En forma matricial

( e i k 1 a e i k 1 a k 1 e i k 1 a k 1 e i k 1 a )( A 1 B 1 )=( e i k 2 a e i k 2 a k 2 e i k 2 a k 2 e i k 2 a )( A 2 B 2 ) ( e i k 2 a e i k 2 a k 2 e i k 2 a k 2 e i k 2 a )( A 2 B 2 )=( e i k 3 a e i k 3 a k 3 e i k 3 a k 3 e i k 3 a )( A 3 B 3 )

Relacionamos los coeficientes A1, B1, con A3 y B3

M k 1 ,a ( A 1 B 1 )= M k 2 ,a ( A 2 B 2 ),( A 2 B 2 )= M k 2 ,a 1 M k 1 ,a ( A 1 B 1 ) M k 2 ,a ( A 2 B 2 )= M k 3 ,a ( A 3 B 3 ), ( A 3 B 3 )= M k 3 ,a 1 M k 2 ,a ( A 2 B 2 )= M k 3 ,a 1 M k 2 ,a M k 2 ,a 1 M k 1 ,a ( A 1 B 1 )

Matriz inversa de una matriz 2×2

( a 11 a 12 a 21 a 22 ) 1 = 1 | a 11 a 12 a 21 a 22 | ( a 22 a 12 a 21 a 11 ) M k 1 = ( e ikx e ikx k e ikx k e ikx ) 1 = 1 2k ( k e ikx e ikx k e ikx e ikx )= 1 2k ( k e ikx e ikx k e ikx e ikx )

Condiciones de contorno

El extremo izquierdo x=-L de la cuerda está sujeto y el extremo derecho x=-L, también. Se tiene que cumplir que

{ y 1 (L)= A 1 exp( i k 1 L )+ B 1 exp( i k 1 L )=0, A 1 = B 1 exp( 2i k 1 L ) y 3 (L)= A 3 exp( i k 3 L )+ B 3 exp( i k 3 L )=0, B 3 = A 3 exp( 2i k 3 L )

Las frecuencias de los modos de vibración se obtienen

( A 3 B 3 )=( m 11 m 12 m 21 m 22 )( A 1 B 1 ),m= M k 3 ,a 1 M k 2 ,a M k 2 ,a 1 M k 1 ,a ( A 3 A 3 exp( 2i k 3 L ) )=( m 11 m 12 m 21 m 22 )( B 1 exp( 2i k 1 L ) B 1 ) { A 3 = m 11 B 1 e 2i k 1 L + m 12 B 1 A 3 e 2i k 3 L = m 21 B 1 e 2i k 1 L + m 22 B 1 ( m 11 B 1 e 2i k 1 L + m 12 B 1 ) e 2i k 3 L = m 21 B 1 e 2i k 1 L + m 22 B 1 m 11 e 2i( k 1 + k 3 )L m 12 e 2i k 3 L + m 21 e 2i k 1 L m 22 =0

Frecuencias de los modos de vibración

Las frecuencias ω de los modos de vibración son aquellas que hacen que la parte real e imaginaria de la expresión

m 11 e 2i( k 1 + k 3 )L m 12 e 2i k 3 L + m 21 e 2i k 1 L m 22 =0

sean nulas simultáneamente

Sea una cuerda hecha de tres segmentos cuyas velocidades de propagación son v1=1, v2=2, v3=1, respectivamente la mitad de la longitud de la cuerda es L=3 y a=1. Representamos la parte real e imaginaria de dicha expresión en el intervalo de frecuencias ω de 0.5 a 2 y marcamos sus intersecciones en el ej X

function cuerda
    v1=1; %velocidades de propagación
    v2=2;
    v3=1;
    a=1;
    L=3;
    %matriz
    f=@(k,x) [exp(1i*k*x), exp(-1i*k*x); k*exp(1i*k*x), -k*exp(-1i*k*x)];
    %matriz inversa
    g=@(k,x) [k*exp(-1i*k*x), exp(-1i*k*x); k*exp(1i*k*x), -exp(1i*k*x)]/(2*k);
    hold on
    fplot(@frecuencias_r,[0.5,2])
    fplot(@frecuencias_i,[0.5,2])
    hold off
    grid on
    xlabel('\omega')
    ylabel('f(\omega)')
    title('Parte real e imaginaria')
  
  function res=frecuencias_r(w) 
        m=g(w/v3,a)*f(w/v2,a)*g(w/v2,-a)*f(w/v1,-a);
        temp=m(1,1)*exp(2*1i*(1/v1+1/v3)*w*L)-m(1,2)*exp(2*1i*w*L/v3)+
m(2,1)*exp(2*1i*w*L/v1)-m(2,2);
        res=real(temp);
    end
    function res=frecuencias_i(w) 
        m=g(w/v3,a)*f(w/v2,a)*g(w/v2,-a)*f(w/v1,-a);
        temp=m(1,1)*exp(2*1i*(1/v1+1/v3)*w*L)-m(1,2)*exp(2*1i*w*L/v3)+
m(2,1)*exp(2*1i*w*L/v1)-m(2,2);
        res=imag(temp);
    end
  end

Los ceros de la parte real e imaginaria coinciden para las frecuencias ω=0.69 y 1.22, aproximadamente

Utilizamos la función fzero de MATLAB para calcularlas en el intervalo en el que la función cambia de signo

function cuerda1
    v1=1; %velocidades de propagación
    v2=2;
    v3=1;
    a=1;
    L=3;
     %matriz
    f=@(k,x) [exp(1i*k*x), exp(-1i*k*x); k*exp(1i*k*x), -k*exp(-1i*k*x)];
    %matriz inversa
    g=@(k,x) [k*exp(-1i*k*x), exp(-1i*k*x); k*exp(1i*k*x), -exp(1i*k*x)]/(2*k);
    
    %raíces
    xb=buscar_intervalos(@frecuencias_r,0.5,4, 40);
    nb=size(xb);
    for s=1:nb(1)
        frec_r(s)=fzero(@frecuencias_r,[xb(s,1) xb(s,2)]);
    end
    disp(frec_r) 
    xb=buscar_intervalos(@frecuencias_i,0.5,4, 40);
    nb=size(xb);
    for s=1:nb(1)
        frec_i(s)=fzero(@frecuencias_i,[xb(s,1) xb(s,2)]);
    end
    disp(frec_i)

    function res=frecuencias_r(w) 
        m=g(w/v3,a)*f(w/v2,a)*g(w/v2,-a)*f(w/v1,-a);
        temp=m(1,1)*exp(2*1i*(1/v1+1/v3)*w*L)-m(1,2)*exp(2*1i*w*L/v3)
+m(2,1)*exp(2*1i*w*L/v1)-m(2,2);
        res=real(temp);
    end
    function res=frecuencias_i(w) 
        m=g(w/v3,a)*f(w/v2,a)*g(w/v2,-a)*f(w/v1,-a);
        temp=m(1,1)*exp(2*1i*(1/v1+1/v3)*w*L)-m(1,2)*exp(2*1i*w*L/v3)
+m(2,1)*exp(2*1i*w*L/v1)-m(2,2);
        res=imag(temp);
    end
    
    function xb = buscar_intervalos(f,a,b,n)
        x = a:(b-a)/n:b;
        j = 0; 
        y1=f(x(1));
        xb=[];
        for i = 1:length(x)-1
            y2=f(x(i+1));
            if sign(y1) ~= sign(y2) 
                j = j + 1;
                xb(j,1) = x(i);
                xb(j,2) = x(i+1);
            end
            y1=y2;
        end
        if isempty(xb) 
            disp('no se han encontrado cambios de signo')
        else
            disp(['número de intervalos:' int2str(j)]) 
        end
    end
end
número de intervalos:11
0.5236  0.6957  1.0472  1.1216  1.5708  2.0200  2.0944  2.4459  2.6180  
3.6652  3.8373
número de intervalos:13
0.6957  0.7854  1.1216  1.3090  1.8326  2.0200  2.3562  2.4459   2.8798  3.1416  
3.4034  3.8373  3.9270

Coinciden las frecuencias, ω=0.6957, 1.1216, 2.0200, 2.4459, 3.8373

Amplitud de la onda estacionaria

Representamos la amplitud de la onda estacionaria y1(x) en el intervalo (-L,-a), y2(x) en el intervalo (-a,-a), e y3(x) en el intervalo (-a, L)

{ y 1 (x)= A 1 exp( i k 1 x )+ B 1 exp( i k 1 x ), k 1 = ω v 1 y 2 (x)= A 2 exp( i k 2 x )+ B 2 exp( i k 2 x ), k 2 = ω v 2 y 3 (x)= A 3 exp( i k 3 x )+ B 3 exp( i k 3 x ), k 3 = ω v 3

Asignamos A3=1. Calculamos los coeficientes A1, B1, A2, B2

( A 2 B 2 )= M k 2 ,a 1 M k 3 ,a ( A 3 B 3 ), B 3 = A 3 exp( 2i k 3 L ) ( A 1 B 1 )= M k 1 ,a 1 M k 2 ,a ( A 2 B 2 )

Para que todos los modos de vibración se representen a la misma escala, se calcula A3 haciendo que el área bajo la curva del cuadrado de la amplitud y(x) multiplicada por la densidad en cada tramo sea la unidad

μ 1 L a | y 1 (x) | 2 dx + μ 2 a a | y 2 (x) | 2 dx + μ 3 a L | y 3 (x) | 2 dx =1,v= T μ 1 v 1 2 L a | y 1 (x) | 2 dx + 1 v 2 2 a a | y 2 (x) | 2 dx + 1 v 3 2 a L | y 3 (x) | 2 dx =1

La integral de cada tramo se calcula de la siguiente forma

x 1 x 2 ( A e ikx +B e ikx ) ( A * e ikx + B * e ikx )dx= x 1 x 2 ( A A * +B B * +A B * e 2ikx + A * B e 2ikx ) dx= ( A A * +B B * )( x 2 x 1 )+ 1 2ik ( A B * ( e 2ik x 2 e 2ik x 1 ) A * B( e 2ik x 2 e 2ik x 1 ) )

donde el símbolo *, indica complejo conjugado

Representamos la amplitud de la onda estacionaria y(x) para las tres primeras frecuencias, ω=0.6957, 1.1216, 2.0200

function cuerda2
    v1=1; %velocidades de propagación
    v2=2;
    v3=1;
    a=1;
    L=3;
    %matriz
    f=@(k,x) [exp(1i*k*x), exp(-1i*k*x); k*exp(1i*k*x), -k*exp(-1i*k*x)];
    %matriz inversa
    g=@(k,x) [k*exp(-1i*k*x), exp(-1i*k*x); k*exp(1i*k*x), -exp(1i*k*x)]/(2*k); 

    m=1;
    for w=[0.6957, 1.1216, 2.0200]
        A3=1;
        B3=-A3*exp(2*1i*w*L/v3);
        V=[A3;B3];
        V=g(w/v2,a)*f(w/v3,a)*V;
        A2=V(1); B2=V(2);
        V=g(w/v1,-a)*f(w/v2,-a)*V;
        A1=V(1); B1=V(2);
        suma=integra(A1,B1,w/v1,-L,-a)/v1^2+integra(A2,B2,w/v2,-a,a)/v2^2+
integra(A3,B3,w/v3,a,L)/v3^2;
        subplot(3,1,m)
        m=m+1;
        hold on
        fc=@(x) real(A1*exp(1i*w*x/v1)+B1*exp(-1i*w*x/v1))/sqrt(suma);
        fplot(fc,[-L,-a])
        fc=@(x) real(A2*exp(1i*w*x/v2)+B2*exp(-1i*w*x/v2))/sqrt(suma);
        fplot(fc,[-a,a])
        fc=@(x) real(A3*exp(1i*w*x/v3)+B3*exp(-1i*w*x/v3))/sqrt(suma);
        fplot(fc,[a,L])
        hold off
        grid on
        xlabel('x')
        ylabel('y')
        title(sprintf('w=%1.3f',w))
    end

    function res=integra(A, B, k, x1,x2)
        res=(A*conj(A)+B*conj(B))*(x2-x1)+(A*conj(B)*(exp(2*1i*k*x2)-
exp(2*1i*k*x1))-conj(A)*B*(exp(-2*1i*k*x2)-exp(-2*1i*k*x1)))/(2*1i*k);
    end
end

Cuerda asimétrica

Las densidades o las velocidades de propagación de los segmentos extremos son iguales v1=v3, por lo que y(x) es simétrica. Consideremos el caso en el que v1=1, v2=2, y v3=3

Representamos la parte real e imaginaria en el intervalo de frecuencias ω de 0.5 a 2 y marcamos las intersecciones en el eje X

Los ceros de la parte real e imaginaria coinciden para las frecuencias ω=0.81 y 1.71, aproximadamente

Las calculamos

número de intervalos:7
1.5708  1.7163  2.3562  2.5094  3.1416  3.5726  3.9270
número de intervalos:8
0.7978  1.1781  1.7163  1.9635  2.5094  2.7489  3.5343  3.5726

Coinciden las frecuencias, ω=1.7163, 2.5094, 3.5726

Como vemos en la figura anterior, para la frecuencia ω=0.7978, la parte real (color azul) es tangente al eje X, el procedimiento numérico es incapaz de localizar esta raíz

Representamos la amplitud de la onda estacionaria y(x) para las tres primeras frecuencias, ω=0.7978, 1.7163, 2.5094

Referencias

Clendenning L. M. A laboratory approach to eigenvalue problem. Am. J. Phys. 36 (10) October 1968, pp. 879-881

Lim Yung-kuo. Problems and Solutions on Mechanics. World Scientific (1994). Problem 1236, pp. 389-393

G. Rodríguez Zurita, Ramón Alvarado Bustos, Rubén Alvarado Bustos, L. E. Zavala Ramírez. Modos de oscilación en cuerdas homogéneas por tercios. Revista Mexicana de Física, 48 (5) Octubre 2002. pp. 463–474.