Propagación de las ondas en un medio dispersivo

En primer lugar, vamos a obtener el resultado de la integral

I= exp( u 2 2 σ 2 +iuy ) du= exp( ( u 2 σ i yσ 2 ) 2 y 2 σ 2 2 ) du= exp( y 2 σ 2 2 ) exp( ( u 2 σ i yσ 2 ) 2 ) du

Haciendo el cambio de variable

x= u 2 σ i yσ 2 dx= 1 2 σ du

Obtenemos

I= 2 σexp( y 2 σ 2 2 ) exp( x 2 ) dx

Sabiendo el resultado de la integral (véase la página anterior)

exp( x 2 ) dx= π

Obtenemos el resultado

I= 2π σexp( y 2 σ 2 2 )

que utlizaremos más adelante.

Transformada de Fourier

Ya hemos estudiado la transformada de Fourier de la función coseno.

f(x)=cos( k 0 x)= 1 2 ( exp(i k 0 x)+exp(i k 0 x) ) F(k)= f(x)exp(ikx)dx= 1 2 ( exp( i(k k 0 )x )dx+ exp( i(k+ k 0 )x )dx )= π( δ(k k 0 )+δ(k+ k 0 ) )

Dos funciones delta de Dirac situadas en +k0 y en -k0.

>> syms k0 x;
>> fx=cos(k0*x);
>> Fk=fourier(fx)
Fk =pi*(dirac(k0 - w) + dirac(k0 + w))

Vamos a obtener la transformada de Fourier de la función

f(x)=exp( ( x x 0 ) 2 2 σ x 2 )cos( k 0 x ) = 1 2 exp( ( x x 0 ) 2 2 σ x 2 )( exp(i k 0 x)+exp(i k 0 x) )

En primer lugar, obtenemos la transformada de Fourier F1(k) de la función f1(x)

f 1 (x)=exp( i k 0 x ( x x 0 ) 2 2 σ x 2 ) F 1 (k)= f(x)exp(ikx)dx = exp( ikx+i k 0 x ) exp( ( x x 0 ) 2 2 σ x 2 )dx F 1 (k)=exp( i(k k 0 ) x 0 ) exp( i( k 0 k)(x x 0 ) ( x x 0 ) 2 2 σ x 2 ) dx

Tenemos una integral del tipo descrito al principio de esta página con

u=x x 0 du=dx y= k 0 kσ= σ x F 1 (k)=exp( i(k k 0 ) x 0 ) exp( iyu u 2 2 σ 2 ) du F 1 (k)= 2π σ x exp( i k 0 x 0 )exp( ik x 0 )exp( ( k 0 k ) 2 σ x 2 2 )

De modo análogo, obtenemos la transformada de Fourier de la función f2(x)

f 2 (x)=exp( i k 0 x ( x x 0 ) 2 2 σ x 2 ) F 2 (k)= f(x)exp(ikx)dx = exp( ikxi k 0 x ) exp( ( x x 0 ) 2 2 σ x 2 )dx F 2 (k)=exp( i(k+ k 0 ) x 0 ) exp( i( k 0 k)(x x 0 ) ( x x 0 ) 2 2 σ x 2 ) dx u=x x 0 du=dx y=( k 0 +k)σ= σ x F 2 (k)= 2π σ x exp( i k 0 x 0 )exp( ik x 0 )exp( ( k 0 +k ) 2 σ x 2 2 )

La transformada de Fourier de la función f(x)= f1(x)+ f2(x) es F(k)=(F1(k)+F2(k))/2

x0=0;
k0=2;
s2=5;
x=-10:0.05:10;
fx=exp(-(x-x0).^2/(2*s2)).*cos(k0*x);
subplot(2,1,1)
plot(x,fx,'b')
grid on
xlabel('x')
ylabel('f(x)')
title('Función')

k=-5:0.05:5;
Fk=sqrt(2*pi*s2)*exp(-1i*k*x0).*(exp(1i*k0*x0)*exp(-(k-k0).^2*s2/2)...
+exp(-1i*k0*x0)*exp(-(k+k0).^2*s2/2))/2;
subplot(2,1,2)
plot(k,abs(Fk),'r')
grid on
xlabel('k')
ylabel('F(k)')
title('Transformada')

En la parte superior de la figura vemos la función coseno, cuya envolvente (amplitud) es la función de Gauss.

No hay cambio apreciable en el valor absoluto de F(k) cuando se modifica la posición del centro de la función de Gauss, x0 a un valor distinto de cero. La transformada de Fourier de esta función consiste en dos picos centrados en ±k0.

Si incrementamos el parámetro σ2 a una valor grande por ejemplo, 100, vemos que la amplitud de la función f(x) es casi constante, y su transformada de Fourier son dos picos puntiagudos centrados en ±k0 que se van aproximando a una función delta de Dirac a medida que se incrementa σ2

Movimiento ondulatorio armónico

En la página web titulada "Descripción de la propagación" vemos en una animación la propagación de un pulso, cuya forma inicial viene descrito por la función f(x), a lo largo del eje X hacia la derecha con velocidad v y sin distorsionarse (sin cambiar de forma).

Ψ =f(x-vt) describe la propagación de una perturbación representada por la función f(x), sin distorsión, a la largo del eje X, hacia la derecha, con velocidad v.

2 Ψ t 2 = v 2 2 Ψ x 2

Estudiamos un caso particular importante, aquél en el que la función f(x) es una función armónica (seno o coseno).

Ψ(x,t)=Ψ0·sin k(x-vt)

Las características de esta función de dos variables, son las siguientes:

  1. La función seno es periódica y se repite cuando el argumento se incrementa en 2π . La función Ψ(x, t) se repite cuando x se incrementa en 2π/k.

  2. k( x+ 2π k vt )=kxkvt+2π

    Se trata de una función periódica, de periodo espacial o longitud de onda λ =/k. La magnitud k se denomina número de onda.

  3. Cuando se propaga un movimiento ondulatorio armónico, un punto x del medio describe un Movimiento Armónico Simple de amplitud Ψ0 y frecuencia angular ω =kv.

  4. Ψ(x,t)=Ψ0·sin (kx-ω t)

    El periodo de la oscilación es P= y la frecuencia  f =1/P.

  5. La igualdad ω =kv, nos permite relacionar el periodo espacial o longitud de onda λ y el periodo de la oscilación P de un punto del medio.

  6. ω=kv 2π P = 2π λ vλ=vP 

    La longitud de onda λ está relacionada con la frecuencia f de la forma λ=v/f. Para una velocidad de propagación v, cuanto mayor es la longitud de onda menor es la frecuencia y viceversa.

x=-6:0.05:6;
k=2*pi/2;
w=2*pi/1;

for i=1:5
    subplot(5,1,i)
    t=(i-1)/4;
    y=sin(k*x-w*t);
    plot(x,y,'r');
    grid on
    ylim([-1.1,1.1])
    str=sprintf('t=%1.2f',t);
    title(str)
end

En la figura vemos que la longitud de onda es λ=2. Un punto del medio, por ejemplo x=0, describe un MAS cuyo periodo es P=1. Un pico (señalado con una flecha) se desplaza una longitud de onda durante un periodo de oscilación.

Superposición de ondas armónicas

Superponemos dos ondas armónicas de la misma amplitud A, la primera onda se propaga con velocidad v1=ω1/k1 y la segunda con velocidad v2=ω2/k2

Ψ 1 (x,t)=Asin( k 1 x ω 1 t) Ψ 2 (x,t)=Asin( k 2 x ω 2 t) Ψ 1 (x,t)+ Ψ 2 (x,t)=2Acos( k 1 k 2 2 x ω 1 ω 2 2 t )sin( k 1 + k 2 2 x ω 1 + ω 2 2 t )

Supongamos que k1 y k2 son números de onda cercanos, por ejemplo k1=12 y k2=10. Una fotografía de la onda resultante en el instante t=0, sería la siguiente

k1=12;
k2=10;
x=0:0.01:6;
y=2*cos((k1-k2)*x/2).*sin((k1+k2)*x/2); %superposición de dos ondas armónicas
y1=2*cos((k1-k2)*x/2); %envolvente
plot(x,y,'b',x,y1,'r',x,-y1,'r');
ylim([-2.1,2.1])
xlabel('x');
ylabel('y')
title('Superposición de ondas')
grid on

La curva en color rojo es la envolvente que varía como 2·cos(x) y la curva en color azul la superposición de las dos ondas en el instante t=0, que varía como 2·cos(x)·sin(11x). Dado que k1 k2 y ω1 ω2 la velocidad de la superposición (en color azul) de las dos ondas es v=ω1/k1 ω2/k2 y se denomina velocidad de fase. Sin embargo, la velocidad de la envolvente es Δωk=(ω1- ω2)/(k1-k2 ), se denomina velocidad de grupo.

Las definiciones de velocidad de fase vp y de grupo vg son, respectivamente

v p = ω k v g = dω dk

Ambas velocidades son funciones de k y son distintas, en general.

fichero = 'onda_1.gif';
hg=figure;
set(hg,'Position',[0,0,568,180]) %posición y tamaño de la ventana gráfica
k1=12; w1=21;
k2=10; w2=20;
x=0:0.02:20;
for t=0:0.05:2*pi
    %superposición de dos ondas armónicas
    y=2*cos((k1-k2)*x/2-(w1-w2)*t/2).*sin((k1+k2)*x/2-(w1+w2)*t/2);
    y1=2*cos((k1-k2)*x/2-(w1-w2)*t/2); %envolvente
    plot(x,y,'b',x,y1,'r',x,-y1,'r');
    ylim([-2.1,2.1])
    xlabel('x');
    ylabel('y')
    title('Superposición de dos ondas')
   
    %GIF animado
    frame=getframe;
    im = frame2im(frame);
    [imind,cm] = rgb2ind(im,256);
    if t==0
        imwrite(imind,cm,fichero,'gif','DelayTime',0,'loopcount',inf);
    else
        imwrite(imind,cm,fichero,'gif','DelayTime',0,'writemode','append');
    end
end

En esta animación, comparamos la velocidad de la envolvente (ω1- ω2)/(k1-k2 )=1/2 (en color rojo) con la velocidad de fase de las ondas (ω1+ω2)/(k1+k2 )=(21+20)/(12+10)≈2 (en color azul)

Ejemplos de medios dispersivos y no dispersivos

Velocidad de grupo

La ecuación de una onda armónica es

Ψ(x,t)= Ψ 0 sin( kxωt) )

donde k es el número de onda, ω es la frecuencia angular, y v=ω/k es la velocidad de propagación

La ecuación de una onda, es la superposición

Ψ(x,t)= a(k)exp( i(kxωt) ) dk

En un medio dispersivo la frecuencia angular ω depende del número de onda k, las ondas de diferentes longitudes de onda λ=2π/k se propgan con distinta velocidad.

En la sección anterior, hemos calculado la transformada de Fourier de la función.

Ψ ( x ,0 ) = exp ( i k 0 x ) · exp ( ( x x 0 ) 2 2 σ x 2 ) a ( k ) = A exp ( i k x 0 ( k k 0 ) 2 2 σ k 2 ) σ x = 1 σ k

Como hemos visto a(k) es una función que tiene un pico agudo centrado en k0, por tanto, a(k) es distinto de cero en las proximidades de k0 y es nula en el resto.

La fase φ(k)=kx-ωt es una función de k. Aproximamos φ(k) alrededor de k0 (donde a(k) es distinto de cero) tomando los primeros términos del desarollo en serie de esta función.

φ(k)= φ 0 + dφ dk | 0 ( k k 0 )+ 1 2 d 2 φ d k 2 | 0 ( k k 0 ) 2 +...

donde

φ 0 = k 0 x ω 0 t dφ dk | 0 =x dω dk | 0 t=x v g t d 2 φ d k 2 | 0 = d 2 ω d k 2 | 0 t=αt φ(k) k 0 x ω 0 t+(x v g t)( k k 0 ) 1 2 αt ( k k 0 ) 2

vg de denomina velocidad de grupo de ondas, mientas que el cociente ω/k se denomina velocidad de fase.

Calculamos la ecuación de la onda en cualquier instante t

Ψ(x,t)=A exp( ik x 0 ( k k 0 ) 2 2 σ k 2 +i( k 0 x ω 0 t+(x v g t)(k k 0 ) 1 2 αt ( k k 0 ) 2 ) ) dk Ψ(x,t)=Aexp( i( k 0 x ω 0 t k 0 x 0 ) ) exp( i(k k 0 ) x 0 ( k k 0 ) 2 2 σ k 2 +i( (x v g t)(k k 0 ) 1 2 αt ( k k 0 ) 2 ) ) dk Ψ(x,t)=Aexp( i( k 0 x ω 0 t k 0 x 0 ) ) exp( i(x x 0 v g t)(k k 0 )(iαt σ k 2 +1) ( k k 0 ) 2 2 σ k 2 ) dk

Hacemos el cambio de variable y tenemos una integral del tipo descrito al principio de esta página

u=k k 0 ,du=dk y=x x 0 v g t, 1 σ 2 = iαt σ k 2 +1 σ k 2 =iαt+ 1 σ k 2 Ψ(x,t)=Aexp( i( k 0 x ω 0 t k 0 x 0 ) ) exp( yu u 2 2 σ 2 ) du Ψ(x,t)= 2π σ k 2 iαt σ k 2 +1 Aexp( i( k 0 x ω 0 t k 0 x 0 ) )exp( ( x x 0 v g t ) 2 σ k 2 2(iαt σ k 2 +1) ) Ψ(x,t)= 2π σ k 2 iαt σ k 2 +1 Aexp( i( k 0 x ω 0 t k 0 x 0 ) )exp( ( x x 0 v g t ) 2 (1iα σ k 2 t) 2( α 2 σ k 2 t 2 + 1 σ k 2 ) )

Propagación de un pulso sin distorsión, α=0, velocidad de grupo vg=2. El punto de color azul marca el centro del grupo que se mueve con velocidad constante vg.

x0=0; %posición inicial de la función de Gauss
k0=2; %número de onda
sx2=5; %extensión de la función de Gauss (cuadrado de sigma)
w0=4*k0; %velocidad de fase w0/k0
vg=2; %velocidad de grupo
alfa=0; %dispersión
sk2=1/sx2;
x=-8:0.02:30;

for i=1:5
    subplot(5,1,i)
    t=(i-1)*3;
    y=sqrt(2*pi*sk2/(1i*alfa*t*sk2+1))*exp(1i*(k0*x-w0*t-k0*x0))...
.*exp(-(x-x0-vg*t).^2*(1-1i*alfa*sk2*t)/(2*(alfa^2*sk2*t^2+1/sk2)));
    hold on
    plot(x,real(y),'r'); %movimiento ondulatorio
    plot(vg*t,0,'bo','markersize',2,'markerfacecolor','b') %centro del grupo
    hold off
    grid on
    ylim([-1.1,1.1])
    str=sprintf('t=%1.2f',t);
    title(str)
end

Propagación del pulso con distorsión α=1, velocidad de grupo vg=2. El punto de color azul marca el centro del grupo.

Animamos la segunda figura

fichero = 'onda_2.gif';
hg=figure;
%posición y tamaño de la ventana gráfica
set(hg,'Position',[0,200,690,180]) 

x0=0; %posición inicial de la función de Gauss
k0=2; %número de onda
sx2=5; %extensión de la función de Gauss (cuadrado de sigma)
w0=4*k0; %velocidad de fase w0/k0
vg=2; %velocidad de grupo
alfa=1; %dispersión
sk2=1/sx2;
x=-8:0.02:30;
for t=0:0.1:10;
    y=sqrt(2*pi*sk2/(1i*alfa*t*sk2+1))*exp(1i*(k0*x-w0*t-k0*x0))...
.*exp(-(x-x0-vg*t).^2*(1-1i*alfa*sk2*t)/(2*(alfa^2*sk2*t^2+1/sk2)));
    plot(x,real(y),'r')
    ylim([-1.1,1.1])
    grid on
    
    %GIF animado
    frame=getframe;
    im = frame2im(frame);
    [imind,cm] = rgb2ind(im,256);
    if t==0
        imwrite(imind,cm,fichero,'gif','DelayTime',0,'loopcount',inf);
    else
        imwrite(imind,cm,fichero,'gif','DelayTime',0,'writemode','append');
    end
end

Ejemplo

La ecuación diferencial que describe la vibración de una barra elástica

λ 2 4 y x 4 + 2 y t 2 =0

Buscamos una solución de la forma

y(x,t)=Aexp( i( kxωt ) ) ω=λ k 2

Obtenemos una relación entre la frecuencia angular ω y el número de onda k que no es lineal (como en un medio no dispersivo)

La solución de la ecuación diferencial es la superposición

y(x,t)= A(k)exp( i( kxω(k)t ) )

Consideremos el caso, en el que la amplitud A(k) es no nula en un pequeño intervalo

A(k)={ A 0 , k 0 Δk<k< k 0 +Δk 0,el resto

El desplazamiento de los puntos de la barra es

y(x,t)= A 0 k 0 Δk k 0 +Δk cos( kxλ k 2 t )dk

Esta integral no es sencilla. El resultado que proporciona Math Symbolic de MATLAB es el siguiente. Se ha sustituido el parámetro λ por la variable a.

>> syms x k a t;
>> z=int(cos(k*x-a*k^2*t),k)
z =pi^(1/2)*cos(x^2/(4*a*t))*fresnelc((2^(1/2)*(x/2 - 
a*k*t)*(-1/(a*t))^(1/2))/pi^(1/2))*(-1/(2*a*t))^(1/2) - 
pi^(1/2)*sin(x^2/(4*a*t))*fresnels((2^(1/2)*(x/2 - a*k*t)*
(-1/(a*t))^(1/2))/pi^(1/2))*(-1/(2*a*t))^(1/2)

Mediante el comando

>> latex(z)

La respuesta ans se transforma en código Latex que copiamos y pegamos en el programa MathType, que a su vez, lo transforma en código MathML

π cos( x 2 4λt )C( 2 ( x 2 λkt ) 1 λt π ) 1 2λt π sin( x 2 4λt )S( 2 ( x 2 λkt ) 1 λt π ) 1 2λt

Teniendo en cuenta que

z =i z C(ix)=iC(x) S(ix)=iS(x)

donde

S(z)= 0 z sin( π t 2 2 )dt ,C(z)= 0 z cos( π t 2 2 )dt

son las integrales de Fresnel

El desplazamiento de los puntos de la barra elástica es

y(x,t)= A 0 π 2λt { cos( x 2 4λt )C( x2λkt 2πλt )+sin( x 2 4λt )S( x2λkt 2πλt ) } k 0 Δk k 0 +Δk = A 0 π 2λt { cos( x 2 4λt ){ C( x2λ( k 0 +Δk )t 2πλt )C( x2λ( k 0 Δk )t 2πλt ) }+sin( x 2 4λt ){ S( x2λ( k 0 +Δk )t 2πλt )S( x2λ( k 0 Δk )t 2πλt ) } }

Representamos esta función en los instantes t=0.02, 5.02 y 10.02, para los valores de los parámetros

syms t x lambda k0 dk;
y=-sqrt(pi/(2*lambda*t))*(cos(x.^2/(4*lambda*t)).*
(fresnelc((x-2*lambda*t*(k0+dk))/sqrt(2*pi*lambda*t))-
fresnelc((x-2*lambda*t*(k0-dk))/sqrt(2*pi*lambda*t)))+
sin(x.^2/(4*lambda*t)).*(fresnels((x-2*lambda*t*(k0+dk))/
sqrt(2*pi*lambda*t))-fresnels((x-2*lambda*t*(k0-dk))/sqrt(2*pi*lambda*t))));

subplot(3,1,1)
y1=subs(y,{t,lambda,dk, k0},{0.02,10,0.05,1});
fplot(y1,[0,400])
grid on
xlabel('x')
ylabel('y(x)')
title('t=0.02')

subplot(3,1,2)
y1=subs(y,{t,lambda,dk, k0},{5.02,10,0.05,1});
fplot(y1,[0,400])
grid on
xlabel('x')
ylabel('y(x)')
title('t=5.02')

subplot(3,1,3)
y1=subs(y,{t,lambda,dk, k0},{10.02,10,0.05,1});
fplot(y1,[0,400])
grid on
xlabel('x')
ylabel('y(x)')
title('t=10.02')

Referencias

David Romero-Abad, Roberto Suárez-Córdova. Energy transported by mechanical waves - two interesting examples: waves on a string and waves on a bar. No publicado