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 se analiza 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.
