Ecuación de ondas (Raúl, Sofía, Jaime)
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Ecuación de ondas |
| Asignatura | EDP |
| Curso | 2023-24 |
| Autores | Raúl Ortega
Sofía Gómez Jaime Sáenz de Miera |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
En este artículo vamos a estudiar la ecuación de ondas mediante el estudio del movimiento de una cuerda. Para ello, vamos a ver cómo varía la solución en función de las condiciones iniciales...
Contenido
1 Ecuación de ondas
Vamos a considerar una cuerda vibrante en el intervalo [0,1], con densidad d y tensión [math] \tau _0[/math] constante, tal que la velocidad de propagación es: [math] c=\frac{\tau_0}{d}=1[/math]. Además, vamos a suponer que la cuerda está fija en los extremos con [math] u_0(x)[/math] y [math] u_1(x)[/math] su posición e impulso iniciales respectivamente. Teniendo todo esto en cuenta y la ecuación de ondas, obtenemos el siguiente sistema:
[math] \begin{cases} u_{tt}-u_{xx}=0,\\ u(0,t)=u(1,t)=0,\\ u(x,0)=u_0(x),\\ u_t(x,0)=u_1(x). \end{cases} [/math]
Donde [math]u(x)[/math] representa la posición de cada punto de la cuerda en cada momento.
1.1 Resolución
Para resolver este sistema, vamos a utilizar el método de separación de variables, el cual consiste en suponer que [math] u(x,t)=T(t)X(x)[/math], lo cual nos da lugar a resolver dos problemas por separado, uno que depende de t y otro que depende de x. Al hacer esto, obtenemos la siguiente solución:
siendo [math] a_k=\frac{\int_0^1 sin(k\pi x)u_0(x)dx}{\int_0^1 sin^2(k \pi x)dx}[/math] y [math] b_k=\frac{1}{k\pi}\frac{\int_0^1 sin(k\pi x)u_1(x)dx}{\int_0^1 sin^2(k \pi x)dx}[/math]
Ahora, para entender mejor esta solución, vamos a tomar como datos iniciales [math]u_0(x)=e^{-100(x-\frac{1}{2})^2}[/math] y [math] u_1(x)=0[/math], lo cual nos da que la solución se expresa de la siguiente manera:
Para representarla vamos a calcular las integrales necesarias para obtener los [math] a_k [/math] mediante el método del trapecio y a dibujar los primeros 50 términos de la serie.
%Definimos las funciones
u0=@(x) exp(-100*(x-0.5).^2); %Condición inicial
uk=@(x,t,k) sin(k.*pi.*x).*cos(k*pi.*t); %Solución del problema
%Aplicamos la regla del trapecio
a=0; b=1; %Extremos del intervalo
x=linspace(a,b,100); %Vector de espacio
t=linspace(a,2,100); %Vector de tiempo
N=length(x)-1; %Número de puntos
h=(b-a)/N;
w=ones(N+1,1); %Vector de pesos
w(1)=1/2; w(N+1)=1/2;
f_n=zeros(length(x),length(t));
%Suma de los 10 primeros términos en todo el vector de tiempo
for i=1:length(t)
for k=1:50
g=u0(x)*sin(k*pi*x)';
a_k=2*h*w'*g; %Resultado de la regla del trapecio (c_k)
f_n(i,:)= f_n(i,:)+a_k.*uk(x,t(i),k);
end
end
%Hacemos un video representando la solución para distintos tiempos
pelicula=VideoWriter("cuerda",'MPEG-4');
pelicula.FrameRate=10;
open(pelicula)
for i=1:length(t)
f=f_n(i,:);
figura=figure(1);
plot(x,f,color='b',LineWidth=1.5)
title(['t = ' num2str(t(i))]);
axis([0 1 -1 1])
imagen=getframe(figura);
writeVideo(pelicula,imagen);
end
close(pelicula)
Vemos que los extremos de la cuerda se mantienen siempre en cero, por lo que la solución satisface las condiciones frontera en todo momento. Además la cuerda empieza dividiéndose en dos y cada parte viaja hacia un lado. Cuando estas partes llegan al final del intervalo se reflejan y cambian el sentido de su desplazamiento.
Ahora si observamos la solución en un intervalo de tiempo más largo [math]t\in[0,4][/math] vemos que la solución hace el mismo recorrido dos veces por lo que podemos concluir que la solución es periódica de periodo 2.
1.2 Onda que viaja en un solo sentido
En este caso tomamos como condiciones iniciales: [math]u_0(x)=e^{-100(x-\frac{1}{2})^2}[/math] y [math] u_1(x)=200(x-\frac{1}{2})e^{-100(x-\frac{1}{2})^2}[/math] que es la derivada de [math]u_0(x)[/math] cambiada de signo. Al igual que en el caso anterior para representarlo vamos a realizar las integrales correspondientes para calcular los coeficientes de la solución mediante el método del trapecio y a representar los primeros 50 términos.
%Definimos las funciones
u0=@(x) exp(-100*(x-0.5).^2); %Condición inicial
u1=@(x) -(100-200.*x).*exp(-100*(x-0.5).^2); %Impulso inicial
uk_cos=@(x,t,k) sin(k.*pi.*x).*cos(k*pi.*t); %Primera parte de la solución del problema
uk_sen=@(x,t,k) sin(k.*pi.*x).*sin(k*pi.*t); %Segunda parte de la solución del problema
%Aplicamos la regla del trapecio
a=0; b=1; %Extremos del intervalo
x=linspace(a,b,100); %Vector de espacio
t=linspace(a,2,100); %Vector de tiempo
N=length(x)-1; %Número de puntos
h=(b-a)/N;
w=ones(N+1,1); %Vector de pesos
w(1)=1/2; w(N+1)=1/2;
f_n=zeros(length(x),length(t));
%Suma de los 10 primeros términos en todo el vector de tiempo
for i=1:length(t)
for k=1:50
g1=u0(x)*sin(k*pi*x)';
g2=-u1(x)*sin(k*pi*x)'/(k*pi);
a_k=2*h*w'*g1; %Resultado de la regla del trapecio (c_k)
b_k=2*h*w'*g2;
f_n(i,:)= f_n(i,:)+a_k.*uk_cos(x,t(i),k)+b_k.*uk_sen(x,t(i),k);
end
end
%Hacemos un video representando la solución para distintos tiempos
pelicula=VideoWriter("cuerda",'MPEG-4');
pelicula.FrameRate=10;
open(pelicula)
for i=1:length(t)
f=f_n(i,:);
figura=figure(1);
plot(x,f)
axis([0 1 -2 2])
imagen=getframe(figura);
writeVideo(pelicula,imagen);
end
close(pelicula)
En este caso vemos que efectivamente la cuerda viaja en un solo sentido y sigue manteniendo los extremos fijos en cero. Además, como el intervalo en el que viaja la cuerda tiene longitud uno y lo recorre dos veces en el intervalo de tiempo [math][0,2][/math], podemos deducir que la velocidad a la que viaja la curva es 1. Esto coincide con la constante [math]c[/math] que hemos elegido la cual representa esta velocidad.
También podemos ver que cuando la cuerda llega a un extremo sigue el mismo comportamiento que en el caso anterior, se refleja y cambia de sentido pero manteniendo la misma velocidad.
1.3 Condiciones Neumann
En los dos casos anteriores hemos considerado condiciones frontera de tipo Dirichlet, es decir imponíamos que la solución en la frontera fuera siempre cero. Ahora vamos a poner condiciones frontera de tipo Neumann imponiendo que la derivada espacial de la solución en los extremos sea siempre cero: [math]u_x(0,t)=u_x(1,t)=0[/math].
Para resolverlo procedemos igual que en los apartados anteriores y representamos la solución con el siguiente código.
%Definimos las funciones
u0=@(x) exp(-100*(x-0.5).^2); %Condición inicial
u1=@(x) -(100-200.*x).*exp(-100*(x-0.5).^2); %Impulso inicial
uk_cos=@(x,t,k) cos(k.*pi.*x).*cos(k*pi.*t); %Solución del problema
uk_sen=@(x,t,k) cos(k.*pi.*x).*sin(k*pi.*t); %Solución del problema
%Aplicamos la regla del trapecio
a=0; b=1; %Extremos del intervalo
x=linspace(a,b,100); %Vector de espacio
t=linspace(a,2,100); %Vector de tiempo
N=length(x)-1; %Número de puntos
h=(b-a)/N;
w=ones(N+1,1); %Vector de pesos
w(1)=1/2; w(N+1)=1/2;
f_n=zeros(length(x),length(t));
%Suma de los 100 primeros términos en todo el vector de tiempo
for i=1:length(t)
for k=1:100
g1=u0(x)*cos(k*pi*x)';
g2=-u1(x)*cos(k*pi*x)'/(k*pi);
a_k=2*h*w'*g1; %Resultado de la regla del trapecio (c_k)
b_k=2*h*w'*g2;
f_n(i,:)= f_n(i,:)+a_k.*uk_cos(x,t(i),k)+b_k.*uk_sen(x,t(i),k);
end
f_n(i,:) = f_n(i,:) + h*w'*u0(x)' ;
end
pelicula=VideoWriter("cuerda",'MPEG-4');
pelicula.FrameRate=10;
open(pelicula)
for i=1:length(t)
f=f_n(i,:);
figura=figure(1);
plot(x,f,color='b',LineWidth=1.5)
axis([0 1 -2 2])
imagen=getframe(figura);
writeVideo(pelicula,imagen);
end
close(pelicula)
En este caso vemos como ya los extremos de la cuerda no están fijos sino que suben y bajan como si viajaran por un carril. En el vídeo se ve cómo en los extremos la solución hace un pico un poco raro, lo que podría dar lugar a pensar que la derivada en los extremos es distinta de cero y que esta solución no satisface las condiciones frontera, pero este comportamiento se debe a que estamos aproximando unas integrales con el método del trapecio y solo estamos pintando 100 términos de la solución.
2 Solución fundamental
En esta sección vamos a analizar la solución fundamental de la ecuación de ondas en dimensiones 1, 2 y 3. La solución fundamental es crucial para comprender cómo una perturbación inicial localizada se propaga en el espacio con el tiempo. Para la ecuación de ondas, la solución fundamental se obtiene resolviendo el siguiente sistema:
[math]\begin{cases} u_{tt}- c^2 \Delta u = 0,& x\in\mathbb{R}^n,t\gt0,\\ u(x,0)=0,u_t(x,0)=\delta(x),&x\in\mathbb{R}^n \end{cases}[/math]
donde \(\delta(x)\) es la delta de Dirac, que representa una perturbación puntual en el origen. La forma de la solución fundamental depende de la dimensión del espacio en el que se propaga la onda, por lo que la vamos a dibujarlas en dimensiones 1, 2 y 3 y, como todas son radiales, lo haremos en su variable radial.
2.1 Dimensión 1
En una dimensión, la solución fundamental está dada por [math]K_1(x, t) = \frac{1}{2c} [H(x + ct) - H(x - ct)][/math]. Donde \(H(s)\) es la función de Heaviside.
clear all
% Parámetros
c = 1; % Velocidad de propagación
t = linspace(0.01, 3, 100); % Tiempo
x = linspace(-3, 3, 1000); % Dominio espacial
% Función de Heaviside
H = @(s) double(s >= 0);
% Solución fundamental en dimensión 1
K1 = @(x, t) (1 / (2 * c)) * (H(x + c * t) - H(x - c * t));
%Lo representamos en un video
pelicula=VideoWriter("fundamental1",'MPEG-4');
pelicula.FrameRate=10;
open(pelicula)
for i=1:length(t)
figura=figure(1);
plot(x,K1(x,t(i)),color='b',LineWidth=1.5);
imagen=getframe(figura);
writeVideo(pelicula,imagen);
end
close(pelicula)
Esta solución describe una onda que se propaga en ambas direcciones a partir del punto de perturbación.
2.2 Dimensión 2
En dos dimensiones, la solución fundamental se expresa como: \(K_2(x, t) = \frac{1}{2\pi c \sqrt{c^2 t^2 - |x|^2}} \chi_{B(0,ct)}(x)\). Donde \(\chi_{B(0,ct)}(x)\) es la función característica de la bola de radio \(ct\) centrada en el origen. Debido a la singularidad en esta expresión, se introduce una regularización:
[math] K_2^\epsilon (x, t) = \frac{1}{\epsilon + 2\pi c \sqrt{c^2 t^2 - |x|^2}} \chi_{B(0,ct)}(x), [/math]
donde \(\epsilon\) es un pequeño valor positivo para evitar la singularidad, que vamos a tomar \(\epsilon=0.01\).
% Parámetros
epsilon = 0.01; % Regularización
c = 1; % Velocidad de propagación
t = linspace(0, 1, 1000); % Tiempo
r = linspace(0, 1, 1000); % Radio
% Solución fundamental regularizada en dimensión 2
K2 = @(r,t) 1./(epsilon+2*pi*c.*sqrt(c^2*t^2-r.^2)).*(r <= c * t);
pelicula=VideoWriter("fundamental2",'MPEG-4');
pelicula.FrameRate=100;
open(pelicula)
for i=1:length(t)
figura=figure(1);
plot(r,K2(r,t(i)),color='b',LineWidth=1.5);
axis([0 1 0 100])
title(['t = ' num2str(t(i))]);
imagen=getframe(figura);
writeVideo(pelicula,imagen);
end
close(pelicula)
En este caso también vemos cómo la solución también representa una onda que se propaga radialmente formando círculos concéntricos alrededor del origen.
2.3 Dimensión 3
En tres dimensiones, la solución fundamental es \(K_3(x, t) = \frac{\delta(|x| - ct)}{4\pi c |x|}\). Donde \(\delta(s)\) es la delta de Dirac. Para representar esta solución en gráficos, se utiliza una aproximación de la delta de Dirac mediante una función gaussiana: \(\delta(s) \approx \phi_k(s) = \sqrt{\frac{k}{\pi}} e^{-ks^2}\), donde \(k\) es un parámetro grande que vamos a tomar \(k=1000\).
% Parámetros
k = 1000; % Parámetro de aproximación
c = 1; % Velocidad de propagación
t = linspace(0.01, 3, 100); % Tiempo
r = linspace(0, 3, 1000); % Radio
% Aproximación de la delta de Dirac
phi_k = @(s) sqrt(k / pi) * exp(-k * s.^2);
% Solución fundamental en dimensión 3
K3 = @(r, t) (phi_k(r - c * t) ./ (4 * pi * c * r));
pelicula=VideoWriter("fundamental3",'MPEG-4');
pelicula.FrameRate=10;
open(pelicula)
for i=1:length(t)
figura=figure(1);
plot(r,K3(r,t(i)),color='b',LineWidth=1.5);
axis([0 3 0 20])
imagen=getframe(figura);
writeVideo(pelicula,imagen);
end
close(pelicula)
En el caso de dimensión 3 vemos que de nuevo la solución representa una onda que se propaga radialmente por culpa de una perturbación formando círculos concéntricos alrededor del origen. A diferencia del caso anterior, aquí se ve cómo la onda va perdiendo intensidad a medida que avanza el tiempo.
2.4 Principio de Huygens
El principio de Huygens es un concepto fundamental en la física de ondas que dice que cada punto de un frente de onda se puede considerar como una fuente secundaria de ondas esféricas. Estas ondas secundarias se superponen para formar el nuevo frente de onda a medida que avanza.
Las soluciones fundamentales de la ecuación de ondas que hemos visto en las secciones anteriores describen cómo una perturbación inicial localizada se propaga a través del espacio. Estas soluciones proporcionan una interpretación del principio de Huygens.
En una dimensión, la perturbación inicial se propaga en dos direcciones opuestas a velocidad c. Cada punto a lo largo de la línea puede ser visto como una fuente de ondas que se propagan hacia la izquierda y hacia la derecha, formando así el nuevo frente de onda. Este comportamiento lineal refleja cómo las ondas se suman en una dimensión para avanzar. En el caso de dos dimensiones, la perturbación inicial se expande radialmente formando círculos concéntricos y en el de tres dimensiones la perturbación inicial se propaga esféricamente.
2.5 algo
clear all
close all
c=1; % Velocidad de propagación.
t=[2]; % Valores de t
theta=linspace(0,2*pi); % intervalo de theta
r=linspace(0,1,50); % intervalo del radio
[TT,RR]=meshgrid(theta,r);
n=size(TT);
U=zeros(n); % Inicializamos la matriz U
XX=RR.*cos(TT);
YY=RR.*sin(TT);
% Solución fundamental en polares
K=@(x,y,r,t)(sqrt(r.^2+y.^2-2*r.*y.*cos(x))<c*t)./(0.01*(c^2*t.^2-r.^2- y.^2 + 2.*r.*y.*cos(x)<10^-16)+2*pi*c*sqrt(abs(c^2*t.^2- r.^2- y.^2 + 2.*r.*y.*cos(x))));
% Solución para cada valor de t
for j=1:length(t)
for i=1:length(r)
U(i,:)=integral2(@(x,y)K(x,y,r(i),t(j)), 0,2*pi, 0,1/2)*ones(1,length(theta));
end
colormap jet
figure(j)
surf(XX,YY,U)
shading interp
title("t="+num2str(t(j)))
end
