Diferencia entre revisiones de «Ecuación de Ondas, (equipo LUA)»

De MateWiki
Saltar a: navegación, buscar
(Se creó una página vacía)
 
 
Línea 1: Línea 1:
 +
{{TrabajoED | Ecuación de ondas. Grupo LUA | [[:Categoría:EDP|EDP]]|[[:Categoría:EDP23/24|2023-24]] | Luis Carreras Hoyos, Lucía Gil Ruiz y Alejandra Hernández Sieber}}
 +
==Introducción==
 +
La ecuación de ondas describe cómo se propaga una perturbación a lo largo del tiempo y el espacio en un medio dado. Su forma más común es una ecuación diferencial parcial de segundo orden, <math> u_{tt}-u_{xx} = 0</math>.
 +
En este trabajo, representaremos diferentes soluciones de la ecuación de ondas en una dimensión y trataremos la solución fundamental de la ecuación de ondas en dimensiones <math>1, 2</math> y <math> 3 </math>. La cual nos servirá para interpretar el principio de Huygens.
  
 +
 +
==Conceptos previos==
 +
 +
'''Velocidad de propagación <math>c</math>''': La velocidad de propagación de una onda se refiere a la velocidad a la cual las perturbaciones o variaciones en la onda se desplazan a través del medio. Esta constante depende de las propiedades del medio en el que se propaga la onda.
 +
 +
'''Principio de Huygens: ''' El principio de Huygens postula que todo punto alcanzado por una onda se comporta como un emisor de ondas. Basándose en este principio y utilizando un método geométrico Huygens explicó las propiedades de las ondas. (Reflexión, refracción, difracción e interferencias).
 +
 +
==Ecuación de ondas en una dimensión==
 +
 +
Consideramos una cuerda vibrante en el intervalo <math>[0,1]</math> con densidad <math>d</math> y tensión <math>\tau_0</math> constante de manera que la velocidad de propagación <math>c = \frac{\tau_0} {d} = 1</math>. Supondremos que la cuerda está fija en los extremos. Llamaremos <math>u_0(x)</math> y <math>u_1(x)</math> su posición e impulso iniciales respectivamente.
 +
 +
De esta manera el sistema de EDPs que modeliza el comportamiento de los desplazamientos transversales de la cuerda es:
 +
 +
<center><math> \begin{cases}
 +
u_{tt}-u_{xx} = 0~~~  \text{  con } x \in [0,1] \text{ y } t>0\\
 +
u(0,t)=u(1,t)=0\\
 +
u(x,0)=u_0(x) \\
 +
u_t(x,0)=u_1(x) \\
 +
\end{cases}</math></center>
 +
 +
Para resolver este sistema aplicaremos separación de variables. De esta manera, suponemos que la solución <math>u(x,t)</math> es de la forma <math>u(t,x)=X(x) \cdot T(t) </math>.  Por consiguiente, se obtienen dos problemas de valor inicial, uno asociado a la parte espacial de la solución <math>X(x)</math> y el otro asociado a la parte temporal <math>T(t)</math>. El primero de ellos será:
 +
 +
<center><math> \begin{cases}
 +
X^{’’}(x)+ \lambda \cdot X(x)~~~ \text{ en } x \in [0,1]\\
 +
X(0)=X(1)=0\\
 +
\end{cases}</math></center>
 +
 +
Y el segundo problema será:
 +
 +
<center><math> \begin{cases}
 +
T^{’’}(t)+ \lambda \cdot T(t)~~~ \text{ con } t>0\\
 +
\end{cases}</math></center>
 +
 +
Teniendo en cuenta los posibles signos del valor de <math> \lambda </math> se llega a que la solución general de <math>u(x,t)</math> es de la forma:
 +
 +
<center><math> u(x,t)=\sum_{k=1}^{\infty} c_k \cdot sen(k \pi x) \cdot cos(k  \pi t) + d_k \cdot sen(k \pi x) \cdot sen(k \pi t) </math></center>
 +
 +
Donde <math> c_k, d_k </math> son los coeficientes de la serie de Fourier que calcularemos imponiendo las condiciones iniciales.
 +
 +
Ahora, suponemos que tomamos como datos iniciales <math> u_{0}(x) = e^{ −100(x−1/2)^{2}} </math> y <math> u_{1}(x) = 0 </math>. Estos datos iniciales los incorporamos a nuestro programa de Matlab y dibujamos la solución en el intervalo de tiempo <math> t ∈ [0, 2] </math> tomando los primeros <math> 50 </math> términos de la serie.
 +
{{matlab|codigo= clear all
 +
close all
 +
 +
% CONDICIONES INICIALES
 +
u_0=@(x) exp(-100.*(x-1/2).^2);
 +
u_1=@(x) 0;
 +
 +
% TÉRMINO K-ÉSIMO DE LA SOLUCIÓN DEL PROBLEMA HABIENDO RESTADO LA SOLUCIÓN ESTACIONARIA
 +
u_k=@(x,t,k,coef_fourier_c,coef_fourier_d) coef_fourier_c.*sin(k.*pi.*x).*cos(k.*pi.*t)+coef_fourier_d.*sin(k.*pi.*x).*sin(k.*pi.*t);
 +
 +
% DATOS PARA LA DISCRETIZACIÓN DEL ESPACIO
 +
extremo_izquierdo=0;
 +
extremo_derecho=1;
 +
paso_espacio=0.001;
 +
 +
% DATOS PARA LA DISCRETIZACIÓN DEL TIEMPO
 +
tiempo_inicial=0;
 +
tiempo_final=2;
 +
paso_tiempo=0.001;
 +
 +
% DISCRETIZACIÓN DEL ESPACIO Y DEL TIEMPO
 +
[X,T]=meshgrid(extremo_izquierdo:paso_espacio:extremo_derecho,tiempo_inicial:paso_tiempo:tiempo_final); % matrices de mallado espacio-temporal
 +
 +
u=zeros(size(X)); % inicializamos la variable w (solución del problema habiendo restado la solución estacionaria)
 +
 +
% SOLUCIÓN EN EL INTERVALO DE TIEMPO [0,1] TOMANDO LOS k_end PRIMEROS TÉRMINOS DE LA SERIE
 +
k_end=50; % último valor de k del sumatorio, siendo k el índice del sumatorio de la solución
 +
 +
for ks=1:k_end % queremos resolver el problema para distintos valores de k
 +
    c_k=trapz(X(1,:),sin(ks.*pi.*X(1,:)).*u_0(X(1,:)))./trapz(X(1,:),(sin(ks.*pi.*X(1,:))).^2);
 +
    d_k=trapz(X(1,:),sin(ks.*pi.*X(1,:)).*u_1(X(1,:)))./trapz(X(1,:),(sin(ks.*pi.*X(1,:))).^2);
 +
    u=u+u_k(X,T,ks,c_k,d_k);
 +
end
 +
grafica=figure(1)
 +
mesh(X,T,u)
 +
colorbar
 +
colormap('jet')
 +
xlabel('Espacio')
 +
ylabel('Tiempo')
 +
foto=getframe(grafica);
 +
title('Representación de la solución')
 +
imwrite(foto.cdata,"Ej1_Ap3.png")
 +
 +
 +
% Crear un vídeo
 +
pelicula=VideoWriter('Ej1_Ap3_Cuerda','MPEG-4'); %creo el video
 +
pelicula.FrameRate=50; %controla velocidad
 +
open(pelicula); %abrimos el vídeo
 +
vector_tiempo=T(:,1)
 +
for t=1:2:length(vector_tiempo)
 +
    figura=figure(2)
 +
    plot(X(1,:),u(t,:),'Linewidth',1.5)
 +
    title(['Cuerda en t = ' num2str(vector_tiempo(t))]);
 +
    axis([0 1 -1 1]);
 +
 +
    % Añadir el frame al vídeo
 +
    imagen=getframe(figura);
 +
    writeVideo(pelicula,imagen); %inserto la imagen
 +
 
 +
end
 +
 +
% Cerrar el vídeo
 +
close(pelicula);
 +
 +
}}
 +
 +
[[Archivo: EQUIPOLUA Ej1 Ap3.png|600px|thumb|center| Representación de la solución obtenida para <math> u_{0}(x) = e^{ −100(x−1/2)^{2}} </math> y <math> u_{1}(x) = 0 </math> en el intervalo de tiempo <math> t \in [0,2] </math> tomando los <math>50</math> primeros términos de la serie.]]
 +
 +
[[Archivo: EQUIPO LUA Ej1 Ap3 Cuerda.gif|600px|thumb|center| Representación de la evolución de la solución obtenida para <math> u_{0}(x) = e^{ −100(x−1/2)^{2}} </math> y <math> u_{1}(x) = 0 </math> en el intervalo de tiempo <math> t \in [0,2] </math> tomando los <math>50</math> primeros términos de la serie.]]
 +
 +
Para el código, hemos comenzado definiendo las condiciones iniciales de la onda “u_0” y “u_1”. Luego, hemos definido una función “u_k”que representa el término k-ésimo de la solución del problema. A continuación, se discretiza el espacio y el tiempo para generar una malla espacio-temporal. Se inicializa la solución u y se calcula iterativamente la contribución de los primeros k términos a la solución. Después de obtener la solución en toda la malla, se grafica en forma de superficie tridimensional. Finalmente, se crea un video que muestra la evolución de la onda en el tiempo.
 +
 +
Podemos ver como la solución es periódica de tiempo <math> T=2 </math> . Al inicio, se tiene una onda con una oscilación de valor máximo <math>1 </math> en <math>x=0.5</math>.
 +
 +
Con el paso del tiempo, la onda se divide en dos y cada una se va trasladando a uno de los extremos de la cuerda. Cuando las ondas alcanzan los extremos se trasladan al eje <math> y </math> negativo, produciéndose el mismo movimiento de manera simétrica e inversa con respecto a este eje. 
 +
 +
Las ondas vuelven a la frontera, se trasladan otra vez al eje <math> y </math> positivo original y se realiza el mismo movimiento transformándose ambas ondas en la onda original. Vuelve al estado inicial en tiempo <math> t=2 </math>.
 +
 +
Además, hemos representado el movimiento de la cuerda en <math> 3D </math>. En el eje horizontal de la gráfica se representa tanto la posición espacial como temporal de la onda. Vemos como en tiempo <math>t=0</math> la cuerda se encuentra en la posición inicial del video. A medida que pasa el tiempo, su posición  máxima de <math>1</math> unidad en el centro del espacio se va trasladando conservando siempre el valor  <math>0</math> en los extremos de la cuerda.
 +
 +
 +
===Cambio en los valores iniciales: Una onda que viaja en un sentido===
 +
 +
Suponemos ahora que tomamos como datos iniciales los correspondientes a una onda que viaja en un solo sentido, es decir, la solución <math>u(x,t)</math> es de la forma <math> u(x,t)=f(x-t)</math>.
 +
 +
Para ello, tomamos <math> u(x,0)=u_0(x)=f(x) </math> y <math> u_t(x,0)=u_1(x)=-f’(x) </math>, siendo <math> u_0(x)</math> como en el caso anterior, es decir, <math> u_0(x)=f(x)= e^{ −100(x−1/2)^{2}}</math>.
 +
 +
Implementando la solución en Matlab hemos obtenido los siguientes resultados:
 +
 +
{{matlab|codigo= clear all
 +
close all
 +
 +
% CONDICIONES INICIALES
 +
u_0=@(x) exp(-100.*(x-1/2).^2);
 +
u_1=@(x) 200.*(x-1/2).*exp(-100.*(x-1/2).^2);
 +
 +
% TÉRMINO K-ÉSIMO DE LA SOLUCIÓN DEL PROBLEMA HABIENDO RESTADO LA SOLUCIÓN ESTACIONARIA
 +
u_k=@(x,t,k,coef_fourier_c,coef_fourier_d) coef_fourier_c.*sin(k.*pi.*x).*cos(k.*pi.*t)+coef_fourier_d/(k.*pi).*sin(k.*pi.*x).*sin(k.*pi.*t);
 +
 +
% DATOS PARA LA DISCRETIZACIÓN DEL ESPACIO
 +
extremo_izquierdo=0;
 +
extremo_derecho=1;
 +
paso_espacio=0.001;
 +
 +
% DATOS PARA LA DISCRETIZACIÓN DEL TIEMPO
 +
tiempo_inicial=0;
 +
tiempo_final=2;
 +
paso_tiempo=0.001;
 +
 +
% DISCRETIZACIÓN DEL ESPACIO Y DEL TIEMPO
 +
[X,T]=meshgrid(extremo_izquierdo:paso_espacio:extremo_derecho,tiempo_inicial:paso_tiempo:tiempo_final); % matrices de mallado espacio-temporal
 +
 +
u=zeros(size(X)); % inicializamos la variable w (solución del problema habiendo restado la solución estacionaria)
 +
 +
% SOLUCIÓN EN EL INTERVALO DE TIEMPO [0,1] TOMANDO LOS k_end PRIMEROS TÉRMINOS DE LA SERIE
 +
k_end=50; % último valor de k del sumatorio, siendo k el índice del sumatorio de la solución
 +
 +
for ks=1:k_end % queremos resolver el problema para distintos valores de k
 +
    c_k=trapz(X(1,:),sin(ks.*pi.*X(1,:)).*u_0(X(1,:)))./trapz(X(1,:),(sin(ks.*pi.*X(1,:))).^2);
 +
    d_k=trapz(X(1,:),sin(ks.*pi.*X(1,:)).*u_1(X(1,:)))./trapz(X(1,:),(sin(ks.*pi.*X(1,:))).^2);
 +
    u=u+u_k(X,T,ks,c_k,d_k);
 +
end
 +
grafica=figure(1)
 +
mesh(X,T,u)
 +
colorbar
 +
colormap('jet')
 +
xlabel('Espacio')
 +
ylabel('Tiempo')
 +
foto=getframe(grafica);
 +
title('Representación de la solución')
 +
imwrite(foto.cdata,"Ej1_Ap4.png")
 +
 +
 +
% Crear un vídeo
 +
pelicula=VideoWriter('Ej1_Ap4_Cuerda','MPEG-4'); %creo el video
 +
pelicula.FrameRate=50; %controla velocidad
 +
open(pelicula); %abrimos el vídeo
 +
vector_tiempo=T(:,1)
 +
for t=1:2:length(vector_tiempo)
 +
    figura=figure(2)
 +
    plot(X(1,:),u(t,:),'Linewidth',1.5)
 +
    title(['Cuerda en t = ' num2str(vector_tiempo(t))]);
 +
    axis([0 1 -8 8]);
 +
 +
    % Añadir el frame al vídeo
 +
    imagen=getframe(figura);
 +
    writeVideo(pelicula,imagen); %inserto la imagen
 +
 
 +
end
 +
 +
% Cerrar el vídeo
 +
close(pelicula);
 +
 +
}}
 +
 +
[[Archivo: EQUIPOLUA Ej1 Ap4.png|600px|thumb|center| Representación de la solución obtenida para <math> u_{0}(x) = e^{ −100(x−1/2)^{2}} </math> suponiendo que la onda viaja en un solo sentido.]]
 +
 +
[[Archivo: EQUIPO LUA Ej1 Ap4 Cuerda.gif|600px|thumb|center| Representación de la evolución de la solución obtenida para <math> u_{0}(x) = e^{ −100(x−1/2)^{2}} </math> suponiendo que la onda viaja en un solo sentido]]
 +
 +
En la gráfica, la solución muestra un comportamiento oscilatorio, típico de las ondas, con los puntos máximos y mínimos moviéndose a través del dominio espacial a medida que el tiempo avanza. Además, se puede apreciar como la solución viaja a velocidad constante con velocidad <math> c=1 ~ m/s</math>. Por otro lado, se puede observar como cuando la onda llega a la frontera rebota y se comienza a propagar con sentido y amplitud contrarios. Esto ocurre porque para mantener <math>u(0,t)=u(1,t)=0</math> , la única forma en que la onda puede seguir cumpliendo estas condiciones después de la reflexión es si la parte reflejada de la onda tiene la amplitud opuesta. Se produce ese comportamiento de manera constante en el tiempo. Hay que destacar que, en tiempo <math> t=2 </math> la onda ha regresado a su posición original.
 +
 +
===Cambio en condiciones frontera: Incorporación de condiciones de frontera tipo Neumann===
 +
 +
 +
En este último caso, vamos a realizar un cambio en las condiciones frontera del sistema. Incorporamos condiciones fronteras tipo Neumann para sustituir a las condiciones tipo Dirichlet presentes en el problema original. De esta manera, tendremos el siguiente problema:
 +
<center><math> \begin{cases}
 +
u_{tt}-u_{xx} = 0~~~  \text{  con } x \in [0,1] \text{ y } t>0\\
 +
u_x(0,t)=u_x(1,t)=0 \\
 +
u(x,0)=u_0(x) \\
 +
u_t(x,0)=u_1(x) \\
 +
 +
\end{cases}</math></center>
 +
 +
Implementamos la solución en Matlab y obtenemos el siguiente resultado:
 +
 +
{{matlab|codigo= clear all
 +
close all
 +
 +
% CONDICIONES NEUMANN
 +
u_0=@(x) exp(-100.*(x-1/2).^2);
 +
u_1=@(x) 200.*(x-1/2).*exp(-100.*(x-1/2).^2);
 +
 +
% TÉRMINO K-ÉSIMO DE LA SOLUCIÓN DEL PROBLEMA
 +
u_k=@(x,t,k,coef_fourier_c,coef_fourier_d) coef_fourier_c.*cos(k.*pi.*x).*cos(k.*pi.*t)+coef_fourier_d/(k.*pi).*cos(k.*pi.*x).*sin(k.*pi.*t);
 +
 +
% DATOS PARA LA DISCRETIZACIÓN DEL ESPACIO
 +
extremo_izquierdo=0;
 +
extremo_derecho=1;
 +
paso_espacio=0.001;
 +
 +
% DATOS PARA LA DISCRETIZACIÓN DEL TIEMPO
 +
tiempo_inicial=0;
 +
tiempo_final=2;
 +
paso_tiempo=0.001;
 +
 +
% DISCRETIZACIÓN DEL ESPACIO Y DEL TIEMPO
 +
[X,T]=meshgrid(extremo_izquierdo:paso_espacio:extremo_derecho,tiempo_inicial:paso_tiempo:tiempo_final); % matrices de mallado espacio-temporal
 +
 +
u=zeros(size(X)); % inicializamos la variable w (solución del problema habiendo restado la solución estacionaria)
 +
 +
% SOLUCIÓN EN EL INTERVALO DE TIEMPO [0,1] TOMANDO LOS k_end PRIMEROS TÉRMINOS DE LA SERIE
 +
k_end=50; % último valor de k del sumatorio, siendo k el índice del sumatorio de la solución
 +
 +
% Añadimos a la solución el término cuando k=0
 +
c_0=trapz(X(1,:),cos(0.*pi.*X(1,:)).*u_0(X(1,:)))./trapz(X(1,:),(cos(0.*pi.*X(1,:))).^2);
 +
u=u+c_0.*cos(0.*pi.*X).*cos(0.*pi.*T);
 +
 +
for ks=1:k_end % queremos resolver el problema para distintos valores de k
 +
    c_k=trapz(X(1,:),cos(ks.*pi.*X(1,:)).*u_0(X(1,:)))./trapz(X(1,:),(cos(ks.*pi.*X(1,:))).^2);
 +
    d_k=trapz(X(1,:),cos(ks.*pi.*X(1,:)).*u_1(X(1,:)))./trapz(X(1,:),(cos(ks.*pi.*X(1,:))).^2);
 +
    u=u+u_k(X,T,ks,c_k,d_k);
 +
end
 +
grafica=figure(1)
 +
mesh(X,T,u)
 +
colorbar
 +
colormap('jet')
 +
xlabel('Espacio')
 +
ylabel('Tiempo')
 +
title('Representación de la solución')
 +
foto=getframe(grafica);
 +
imwrite(foto.cdata,"Ej1_Ap5.png")
 +
 +
 +
% Crear un vídeo
 +
pelicula=VideoWriter('Ej1_Ap5_Cuerda','MPEG-4'); %creo el video
 +
pelicula.FrameRate=50; %controla velocidad
 +
open(pelicula); %abrimos el vídeo
 +
vector_tiempo=T(:,1)
 +
for t=1:2:length(vector_tiempo)
 +
    figura=figure(2)
 +
    plot(X(1,:),u(t,:),'Linewidth',1.5)
 +
    title(['Cuerda en t = ' num2str(vector_tiempo(t))]);
 +
    axis([0 1 -8 8]);
 +
 +
    % Añadir el frame al vídeo
 +
    imagen=getframe(figura);
 +
    writeVideo(pelicula,imagen); %inserto la imagen
 +
 
 +
end
 +
 +
% Cerrar el vídeo
 +
close(pelicula);
 +
 +
}}
 +
[[Archivo: EQUIPOLUA Ej1 Ap5.png|600px|thumb|center| Representación de la solución obtenida]]
 +
 +
[[Archivo: EQUIPO LUA Ej1 Ap5 Cuerda.gif|600px|thumb|center| Representación de la evolución de la solución obtenida]]
 +
 +
 +
 +
En este código, se sigue el mismo algoritmo que el anterior, pero cambiando las condiciones de contorno de Neumann en lugar de condiciones de contorno de Dirichlet.
 +
 +
El término k-ésimo de la solución “u_k” se redefine para adaptarse a estas nuevas condiciones de contorno. Nuestra base de Fourier para el problema anterior estaba formada por el conjunto <math> \left\{sin(k \pi x)\right\}</math>. Sin embargo, imponiendo las condiciones de Neuman nuestra nueva base es <math>\left\{1,cos(k \pi x) \right\} </math>. Se muestra gráficamente y se genera un video de la evolución temporal de la onda.
 +
 +
En la gráfica, la solución muestra un comportamiento oscilatorio, esta vez solo de valores máximos moviéndose a través del dominio espacial a medida que el tiempo avanza. La condición nueva <math>u_x(0,t) = u_x(1,t) = 0 </math> indica que la derivada espacial de  u  es cero en los extremos. 
 +
 +
Cabe destacar que no se produce el fenómeno de reflexión con respecto al eje x una vez la onda alcanza el extremo izquierdo de la cuerda, tal y como sucedía en el caso anterior. Esto se debe a que ya no tenemos las condiciones de Dirichlet que fijaban la condición frontera para el valor 0, en este nuevo caso la condición frontera puedo tomar valores distintos de 0.
 +
 +
Este comportamiento, lo podemos apreciar también en el vídeo, donde vemos que la onda se deplaza y alcanzan su máximo valor cuando dicho desplazamiento alcanza los extremos de la cuerda. En conreto,  la solución toma un valor de <math>2</math>.
 +
== Solución fundamental de ondas ==
 +
 +
La solución fundamental de ondas es la solución que se obtiene al dar un impulso inicial muy localizado en <math> x=0 </math>. A lo largo de esta subsección analizaremos la solución fundamental de la ecuación de ondas en dimensiones <math>1, 2</math> y <math>3</math>.
 +
 +
Matemáticamente, la solución fundamental es la solución del sistema relativo a la ecuación de ondas:
 +
 +
<center><math>\left \{ \begin{array}{ll} 
 +
u_{tt}-c^2\Delta u=0, \quad x \in \mathbb{R}^n, t>0, \\
 +
u(x,0)=0, u_t(x,0)=\delta (x), \quad x \in \mathbb{R}^n,
 +
\end{array} \right. 
 +
</math></center>
 +
 +
 +
donde <math> \delta(x) </math> es la delta de Dirac. Formalmente, se define como el siguiente límite:
 +
 +
<center><math>
 +
\delta (x) = \lim \limits_{r \rightarrow 0} \frac{1}{|B(0,r)|} \chi_{B(0,r)}(x) \sim \left \{ \begin{array}{ll}
 +
\infty \quad x=0, \\
 +
0 \quad x\neq 0,
 +
\end{array} \right. 
 +
</math></center>
 +
 +
donde <math>\chi_{B(0,r)}(x) </math> es la función característica de la bola centrada en <math> 0 </math> de radio <math> r </math> y <math>|B(0,r)|</math> es el volumen de la bola.
 +
 +
Nos centraremos ahora en la solución fundamental en dimensiones <math> 1,2 </math> y <math> 3 </math>, dado que son las que Podemos representar con facilidad. Sus expresiones asociadas son:
 +
 +
En dimensión <math>n=1</math>:
 +
 +
<center><math> K_1(x,t) = \frac{1}{2c}[H(x+ct)-H(x-ct)], </math></center>
 +
 +
Donde <math>H(s)</math> es la función de Heaviside,
 +
 +
<center><math> H(s) = \left \{ \begin{array}{ll}
 +
0, \quad x<0 \\
 +
1, \quad x \geq 0
 +
\end{array} \right. 
 +
</math></center>
 +
 +
En dimensión <math> n=2 </math>, la expresión de la solución fundamental es:
 +
 +
<center><math>
 +
K_2(x,t)=\frac{1}{2\pi c \sqrt{c^2t^2-|x|^2}} \chi_{B(0,ct)}(x),
 +
</math></center>
 +
 +
donde <math> \chi_{B(0,ct)}(x) </math> es la función característica de la bola de centro <math> 0 </math> y radio <math> ct </math>.
 +
 +
Por último, en dimensión <math> 3 </math>, la expresión es:
 +
 +
<center><math>
 +
K_3(x,t)=\frac{\delta (|x|-ct)}{4\pi c|x|}, t>0.
 +
 +
</math></center>
 +
 +
=== Solución fundamental de ondas en una dimensión ===
 +
 +
En primer lugar, estudiamos la solución fundamental de ondas en una dimensión. Recordemos que dicha expresión era:
 +
 +
<center><math> K_1(x,t) = \frac{1}{2c}[H(x+ct)-H(x-ct)], </math></center>
 +
 +
Donde <math>H(s)</math> es la función de Heaviside,
 +
 +
<center><math> H(s) = \left \{ \begin{array}{ll}
 +
0, \quad x<0 \\
 +
1, \quad x \geq 0
 +
\end{array} \right. 
 +
</math></center>
 +
 +
Gráficamente la solución es la siguiente:
 +
 +
[[Archivo:EQUIPOLUA Ej2 Ap1 n1.png|600px|thumb|center| Representación de la solución fundamental de ondas en dimensión <math> n=1 </math>]]
 +
 +
{{matlab|codigo= clear all
 +
close all
 +
 +
c=1;
 +
 +
% Definimos la solución fundamental para n=1
 +
k_1=@(x,t) 1/(2.*c).*((x+c.*t>=0)-(x-c.*t>=0));
 +
 +
% Datos inicial y final del espacio
 +
espacio_inicial=-100;
 +
espacio_final=100;
 +
 +
% Datos inicial y final del tiempo
 +
tiempo_inicial=0;
 +
tiempo_final=100;
 +
 +
% Discretización del espacio y del tiempo
 +
espacio=linspace(espacio_inicial,espacio_final,1000);
 +
tiempo=linspace(tiempo_inicial,tiempo_final,1000);
 +
 +
% Mallado
 +
[Espacio,Tiempo]=meshgrid(espacio,tiempo);
 +
 +
% Representamos la solución
 +
grafica=figure(1)
 +
mesh(Espacio,Tiempo,k_1(Espacio,Tiempo))
 +
colorbar
 +
colormap('jet')
 +
xlabel('Espacio')
 +
ylabel('Tiempo')
 +
title('Representación de la solución fundamental para n=1')
 +
foto=getframe(grafica);
 +
imwrite(foto.cdata,"Ej2_Ap1_n1.png")
 +
}}
 +
 +
 +
Podemos observar como la solución obtenida es radial, es decir, no depende de <math> x </math>, únicamente del valor absoluto de x, <math> |x| </math>. Esto puede verse fácilmente observando la paridad de la función en la representación. Analíticamente vemos como:
 +
 +
<center><math> K_1(-x,t) = \frac{1}{2c}[H(-x+ct)-H(-x-ct)] = \frac{1}{2c}[H(-(x-ct))-H(-(x+ct))] = \frac{1}{2c}[-H(x-ct)+H(x+ct)] = \frac{1}{2c}[H(x+ct)-H(x-ct)] = K_1(x,t)</math></center>
 +
Donde en la tercera igualdad se ha usado que, según la definición de la función de Heaviside, <math> H(-s)= 1-H(s) </math>.
 +
 +
Con respecto al código, primero se define la solución fundamental en una dimensión además de los datos necesarios: puntos del espacio y tiempo iniciales y finales. Con ello, se define una discretización espacio-temporal y su correspondiente mallado. Por último, representamos la solución gráficamente e introducimos un comando para que nos guarde directamente la foto en el ordenador.
 +
 +
=== Solución fundamental de ondas en dos dimensiones ===
 +
 +
En segundo lugar, estudiamos la solución fundamental de ondas en dos dimensiones. Su expresión es:
 +
 +
<center><math>
 +
K_2(x,t)=\frac{1}{2\pi c \sqrt{c^2t^2-|x|^2}} \chi_{B(0,ct)}(x),
 +
</math></center>
 +
 +
Podemos ver como la solución es radial ya que solo depende de la norma de <math>x</math>. Esta solución vista como solución radial simplifica considerablemente el problema, en particular permite la representación de la solución con el transcurso del tiempo.
 +
 +
Si implementamos nuestra solución en Matlab, tenemos que añadir la condición de <math>c=\epsilon=0.01</math> para evitar la singularidad.
 +
 +
[[Archivo: EQUIPOLUA Ej2 Ap1 n2.png|600px|thumb|center| Representación de la solución fundamental de ondas en dimensión <math> n=2 </math>]]
 +
 +
{{matlab|codigo= clear all
 +
close all
 +
 +
c=1;
 +
 +
% Definimos la solución fundamental para n=2
 +
k_2=@(r,t) (r<=c*t)./(0.01+2*pi*c.*sqrt(c^2*t.^2-r.^2));
 +
 +
% Datos inicial y final del espacio
 +
radio_inicial=0;
 +
radio_final=1;
 +
 +
% Datos inicial y final del tiempo
 +
tiempo_inicial=0;
 +
tiempo_final=1;
 +
 +
% Discretización del espacio y del tiempo
 +
radio=linspace(radio_inicial,radio_final,1000);
 +
tiempo=linspace(tiempo_inicial,tiempo_final,1000);
 +
 +
% Mallado
 +
[Radio,Tiempo]=meshgrid(radio,tiempo);
 +
 +
% Representamos la solución
 +
grafica=figure(1)
 +
mesh(Radio,Tiempo,k_2(Radio,Tiempo))
 +
colorbar
 +
colormap('jet')
 +
xlabel('Radio')
 +
ylabel('Tiempo')
 +
title('Representación de la solución fundamental para n=2')
 +
foto=getframe(grafica);
 +
imwrite(foto.cdata,"Ej2_Ap1_n2.png")
 +
}}
 +
 +
En esta imagen, apreciamos el comportamiento de la onda, el cual se inicia en el extremo izquierdo con un valor de 100 en <math>(x,t)=(0,0)</math>. Sin embargo, en cuanto aumenta el radio, el valor de la función cambia abruptamente a cero.  Con el paso del tiempo, este valor máximo se desplaza hacia el lado derecho hasta el tiempo <math>t=1</math> donde finalmente alcanza el extremo derecho.
 +
 +
En relación con el código empleado, es muy similar al de para una dimensión. La única diferencia radica en la definición de la solución fundamental. En este caso, se define en ella la función característica en la bola como <math>r</math>=<math>c</math>*<math>t</math> y se define además la regularización de la solución fundamental para evitar la singularidad en su representación.
 +
 +
 +
=== Solución fundamental de ondas en tres dimensiones ===
 +
 +
En este último caso, representamos la solución fundamental de ondas en tres dimensiones.  Nuevamente, recordamos que su expresión era:
 +
 +
<center><math>
 +
K_3(x,t)=\frac{\delta (|x|-ct)}{4\pi c|x|}, t>0 </math></center>
 +
 +
Para representar gráficamente la solución, hemos tenido que optar por emplear la aproximación de la delta de Dirac. Esto es dado que la delta de Dirac presenta una singularidad. Tomamos entonces:
 +
 +
<center><math> \delta (s) \sim \phi _k (s) = \sqrt{\frac{k}{\pi}}e^{-ks^2}, ~~~~~ k>>1, </math> </center>
 +
 +
Concretamente, hemos tomado <math>k=1000</math>. Además, podemos observar como se verifica que <math> \int_{-\infty}^{+\infty} \phi _k (s) ~ \text{ds}=1</math>, criterio que debe verificar la delta de Dirac.
 +
 +
La solución obtenida en este caso es la siguiente:
 +
 +
[[Archivo: EQUIPOLUA Ej2 Ap1 n3.png|600px|thumb|center| Representación de la solución fundamental de ondas en dimensión <math> n=3 </math>]]
 +
 +
{{matlab|codigo= clear all
 +
close all
 +
 +
c=1;
 +
 +
% Definimos la aproximación de la delta de Dirac
 +
aprox_dirac=@(k,s) sqrt(k/pi)*exp(-k*s.^2);
 +
 +
% Definimos la solución fundamental para n=3
 +
k_3=@(r,t) aprox_dirac(1000,r-c.*t)./(4*pi*c.*r);
 +
 +
% Datos inicial y final del espacio
 +
radio_inicial=0;
 +
radio_final=1;
 +
 +
% Datos inicial y final del tiempo
 +
tiempo_inicial=0;
 +
tiempo_final=1;
 +
 +
% Discretización del espacio y del tiempo
 +
radio=linspace(0,1,1000);
 +
tiempo=linspace(0,1,1000);
 +
 +
% Mallado
 +
[Radio,Tiempo]=meshgrid(radio,tiempo);
 +
 +
% Representamos la solución
 +
grafica=figure(1)
 +
mesh(Radio,Tiempo,k_3(Radio,Tiempo))
 +
colorbar
 +
colormap('jet')
 +
xlabel('Radio')
 +
ylabel('Tiempo')
 +
title('Representación de la solución fundamental para n=3')
 +
foto=getframe(grafica);
 +
imwrite(foto.cdata,"Ej2_Ap1_n3.png")
 +
 +
}}
 +
 +
Análogamente, la solución fundamental es radial, dependiendo únicamente de la norma de <math> \textbf{x} </math>.
 +
 +
Por otro lado, podemos observar como la solución fundamental tiene una singularidad en el punto <math>(\textbf{x},t)=(\textbf{0},0) </math> debido a la definición de la delta de Dirac. Al haber tomado la aproximación de la delta de Dirac, esta singularidad no se produce gráficamente, aunque sí que toma valores bastante grandes con respecto a los demás valores espacio-temporales del dominio representado.
 +
 +
Con respecto al código, de nuevo, el funcionamiento es análogo a los dos anteriores. En primer lugar, se define la aproximación de la delta de Dirac y la solución fundamental en dimensión 3, además de los datos necesarios: puntos del espacio y tiempo iniciales y finales. Con ello, se define una discretización espacio-temporal y el mallado a emplear. Por último, representamos la solución gráficamente.
 +
 +
===Solución del sistema usando la solución fundamental para n=2 ===
 +
 +
En esta última sección, se busca calcular y representar la solución del problema en dimensión 2
 +
<center><math>\left \{ \begin{array}{ll} 
 +
u_{tt}-c^2\Delta u=0, \quad x \in \mathbb{R}^n, t>0, \\
 +
u(x,0)=0, u_t(x,0)=h (x)=\chi_{B(0,1/2)}(x) , \quad x \in \mathbb{R}^n,
 +
\end{array} \right. 
 +
</math></center>
 +
que viene dada por la convolución
 +
<math> u(x,t)=\int_{\mathbb{R^2}} K_2(x-y,t) h(y) dy </math> .
 +
 +
En primer lugar, es importante destacar que dicha solución es radial en x pues tanto la solución fundamental en dos dimensiones como la función <math> h(x)=\chi_{B(0,1/2)}(x) </math> son funciones radiales y, con ello, la convolución también será una función radial.
 +
 +
Por lo tanto, si suponemos que <math> x=(r_1 \cdot cos(\theta_1), r_1 \cdot sin(\theta_1)) </math> al usar coordenadas polares, esta propiedad nos permite que a la hora de calcular la solución del problema planteado, basta con calcular la solución cuando <math> \theta_1=0 </math> y posteriormente tomar el mismo valor de la solución para cualquier valor <math> \theta_1 \in [0,2\pi] </math>. 
 +
 +
Retomando la integral con la que definimos la solución, vamos a calcularla usando coordenadas polares, de manera que sean <math> x=(r_1 cos(0), r_1 sin(0)) </math> e <math> y= (r_2 cos(\theta_2),r_2 sin(\theta_2)) </math>, podemos expresar <math> |x-y|^2=r_1^2+r_2^2-2 r_1 r_2 cos(\theta_2) </math>. Con ello, la solución del problema viene dada en coordenadas polares por la siguiente expresión:
 +
 +
<center> <math> u(r_1,t)=\int_{0}^{2 \pi} \int_{0}^{1/2} \frac{1_{B(0,ct)} (r_1 cos(\theta_1)-r_2 cos(\theta_2),-r_2 sin(\theta_2))}{2 \pi c \sqrt{c^2 t^2 – (r_1^2+r_2^2-2 r_1 r_2 cos(\theta_2))}} dr_2 d\theta_2  </math> </center>
 +
Teniendo esto en cuenta, hemos elaborado el siguiente código de Matlab que representa la solución para los tiempos <math> t=0,0.5,1,2 </math> y velocidad de propagación <math> c=1 </math>, sin olvidar que hay que deshacer los cambios de variables realizados.
 +
 +
{{matlab|codigo= clear all
 +
close all
 +
 +
warning('off','all')
 +
 +
c=1;
 +
 +
% Datos inicial y final del espacio
 +
radio_inicial=0;
 +
radio_final=1;
 +
 +
% Datos inicial y final de theta
 +
theta_inicial=0;
 +
theta_final=2*pi;
 +
 +
% Tiempos a estudiar
 +
tiempos=[0,0.5,1,2];
 +
 +
% Discretización del espacio y del tiempo
 +
theta_1=linspace(theta_inicial,theta_final,1000);
 +
radio=linspace(radio_inicial,radio_final,1000);
 +
 +
% Mallado
 +
[Theta,Radio]=meshgrid(theta_1,radio);
 +
 +
u=zeros(size(Theta));
 +
u_eps=zeros(size(Theta));
 +
 +
k_2=@(r_1,r_2,theta_2,t) (sqrt(r_1.^2+r_2.^2-2.*r_1.*r_2.*cos(theta_2))<=c*t)./(2*pi*c.*sqrt(abs(c^2*t.^2-r_1.^2-r_2.^2+2.*r_1.*r_2.*cos(theta_2))));
 +
k_2_eps=@(r_1,r_2,theta_2,t) (sqrt(r_1.^2+r_2.^2-2.*r_1.*r_2.*cos(theta_2))<=c*t)./(0.01*(c^2*t.^2-r_1.^2-r_2.^2+2.*r_1.*r_2.*cos(theta_2)<eps)+2*pi*c.*sqrt(abs(c^2*t.^2-r_1.^2-r_2.^2+2.*r_1.*r_2.*cos(theta_2))));
 +
 +
for j=1:length(tiempos)
 +
    tiempo=tiempos(j);
 +
    for i=1:length(radio)
 +
        u(i,:)=integral2(@(r_2,theta_2)k_2(radio(i),r_2,theta_2,tiempo), 0,1/2, 0,2*pi)*ones(1,length(theta_1));
 +
        u_eps(i,:)=integral2(@(r_2,theta_2)k_2_eps(radio(i),r_2,theta_2,tiempo), 0,1/2, 0,2*pi)*ones(1,length(theta_1));
 +
    end
 +
 +
    % Representación de la solución deshaciendo el cambio de variable
 +
    figura=figure(j)
 +
 +
    set(figura,'position',[2,2,1400,550])
 +
    subplot(1,2,1)
 +
    mesh(Radio.*cos(Theta),Radio.*sin(Theta),u)
 +
    colorbar
 +
    colormap('jet')
 +
    xlabel('x_1')
 +
    ylabel('x_2')
 +
    title("Usando la solución fundamental para n=2")
 +
 +
    subplot(1,2,2)
 +
    mesh(Radio.*cos(Theta),Radio.*sin(Theta),u_eps)
 +
    colorbar
 +
    colormap('jet')
 +
    xlabel('x_1')
 +
    ylabel('x_2')
 +
    title("Usando la regularización para n=2")
 +
 +
    sgtitle("Solución para t="+num2str(tiempo))
 +
 +
    foto=getframe(figura);
 +
    cadena="Ej2_Ap2_t_"+num2str(tiempo)+".png"
 +
    imwrite(foto.cdata,cadena)
 +
end
 +
 +
}}
 +
 +
Inicialmente, se definen los parámetros del problema, como la velocidad de propagación de la onda “c” y los límites del espacio y el ángulo. Luego, se discretiza el espacio y el tiempo, creando una malla.
 +
 +
La función calcula el término k-ésimo de la solución del problema para cada punto en la malla, mientras que “k_2_eps” realiza lo mismo pero con una regularización para evitar singularidades numéricas. Estos términos se integran numéricamente sobre la región de integración correspondiente a cada punto en el espacio y el ángulo, generando así la solución en cada tiempo dado.
 +
 +
La solución se visualiza en dos subgráficos, uno utilizando la solución fundamental y otro utilizando la regularización para comparar los resultados. Este proceso se repite para cada tiempo especificado en el vector "tiempos".
 +
 +
Cabe destacar que este código requiere bastante tiempo de ejecución, pues ha de calcular muchas integrales. Por ese mismo motivo, hemos representado en un primer lugar la solución tomando una discretización para los valores de <math> \theta_1 </math> de 100 puntos y para los valores de <math> r_1 </math> también de 100 puntos. Con ello, hemos obtenido los siguientes resultados.
 +
 +
[[Archivo: EQUIPOLUA Ej2 Ap2 t 0.png|600px|thumb|center| Representación de la solución en tiempo t=0]]
 +
[[Archivo: EQUIPOLUA Ej2 Ap2 t 0.5.png|600px|thumb|center| Representación de la solución en tiempo t=0.5]]
 +
[[Archivo: EQUIPOLUA Ej2 Ap2 t 1.png|600px|thumb|center| Representación de la solución en tiempo t=1]]
 +
[[Archivo: EQUIPOLUA Ej2 Ap2 t 2.png|600px|thumb|center| Representación de la solución en tiempo t=2]]
 +
 +
 +
En primer lugar, para entender estas gráficas es imprescindible destacar el fenómeno que sucede en la representación para <math> t=0.5 </math>.  Para ello, debemos recordar que la solución fundamental en dos dimensiones presentaba una singularidad. Ese es el motivo por el que en la primera gráfica hay una franja en la que aparentemente la solución no está definida. De hecho, lo que está sucediendo en concreto es que para dichos puntos Matlab devuelve el valor NaN pues el denominador de la expresión se anula.
 +
 +
Para solventar este inconveniente en la representación, hemos optado en apoyarnos en la regularización ya mencionada anteriormente en este trabajo. Sin embargo, no hemos aplicado dicha expresión exactamente, pues hemos optado por sumar <math> \epsilon=0.1 </math> cuando el denominador se anula. Esta modificación realizada es la que se puede apreciar en la segunda gráfica de cada una de las imágenes.
 +
 +
Con respecto a la primera de las imágenes, en la que se representa la solución para <math> t=0 </math>, esta nos permite confirmar que la solución cumple la condición inicial impuesta en el problema.
 +
 +
Finalmente, hemos representado la solución tomando una  discretización para los valores de <math> \theta_1 </math> de 1000 puntos y para los valores de <math> r_1 </math> también de 1000 puntos. Los resultados obtenidos han sido los siguientes:
 +
 +
[[Archivo: EQUIPOLUA Ej2 Ap2 t 0 1000.png|600px|thumb|center| Representación de la solución en tiempo t=0]]
 +
[[Archivo: EQUIPOLUA Ej2 Ap2 t 0.5 1000.png|600px|thumb|center| Representación de la solución en tiempo t=0.5]]
 +
[[Archivo: EQUIPOLUA Ej2 Ap2 t 1 1000.png|600px|thumb|center| Representación de la solución en tiempo t=1]]
 +
[[Archivo: EQUIPOLUA Ej2 Ap2 t 2 1000.png|600px|thumb|center| Representación de la solución en tiempo t=2]]
 +
 +
Como se puede apreciar, la singularidad que se producía en la representación anterior ya no es apreciable al tomar esta nueva discretización, pues la malla de puntos es más fina. Realmente, la franja de puntos que no se representan en la gráfica sigue existiendo, pero al ser la malla mucho más fina, dicha región pasa a ser inapreciable.
 +
 +
 +
==Referencias==
 +
*[Partial differential equations in action from modelling to theory. Sandro Salsa]
 +
*[Principio de Huygens: https://www.educa2.madrid.org/web/fmartinezgarcia/oscilaciones-y-ondas/-/book/oscilaciones-y-ondas?_book_viewer_WAR_cms_tools_chapterIndex=a1b08c93-73f1-4a0c-b702-44d2f349c9fc]
 +
 +
[[Categoría:EDP]]
 +
[[Categoría:EDP23/24]]

Revisión actual del 20:35 27 may 2024

Trabajo realizado por estudiantes
Título Ecuación de ondas. Grupo LUA
Asignatura EDP
Curso 2023-24
Autores Luis Carreras Hoyos, Lucía Gil Ruiz y Alejandra Hernández Sieber
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

La ecuación de ondas describe cómo se propaga una perturbación a lo largo del tiempo y el espacio en un medio dado. Su forma más común es una ecuación diferencial parcial de segundo orden, [math] u_{tt}-u_{xx} = 0[/math]. En este trabajo, representaremos diferentes soluciones de la ecuación de ondas en una dimensión y trataremos la solución fundamental de la ecuación de ondas en dimensiones [math]1, 2[/math] y [math] 3 [/math]. La cual nos servirá para interpretar el principio de Huygens.


2 Conceptos previos

Velocidad de propagación [math]c[/math]: La velocidad de propagación de una onda se refiere a la velocidad a la cual las perturbaciones o variaciones en la onda se desplazan a través del medio. Esta constante depende de las propiedades del medio en el que se propaga la onda.

Principio de Huygens: El principio de Huygens postula que todo punto alcanzado por una onda se comporta como un emisor de ondas. Basándose en este principio y utilizando un método geométrico Huygens explicó las propiedades de las ondas. (Reflexión, refracción, difracción e interferencias).

3 Ecuación de ondas en una dimensión

Consideramos una cuerda vibrante en el intervalo [math][0,1][/math] con densidad [math]d[/math] y tensión [math]\tau_0[/math] constante de manera que la velocidad de propagación [math]c = \frac{\tau_0} {d} = 1[/math]. Supondremos que la cuerda está fija en los extremos. Llamaremos [math]u_0(x)[/math] y [math]u_1(x)[/math] su posición e impulso iniciales respectivamente.

De esta manera el sistema de EDPs que modeliza el comportamiento de los desplazamientos transversales de la cuerda es:

[math] \begin{cases} u_{tt}-u_{xx} = 0~~~ \text{ con } x \in [0,1] \text{ y } t\gt0\\ u(0,t)=u(1,t)=0\\ u(x,0)=u_0(x) \\ u_t(x,0)=u_1(x) \\ \end{cases}[/math]

Para resolver este sistema aplicaremos separación de variables. De esta manera, suponemos que la solución [math]u(x,t)[/math] es de la forma [math]u(t,x)=X(x) \cdot T(t) [/math]. Por consiguiente, se obtienen dos problemas de valor inicial, uno asociado a la parte espacial de la solución [math]X(x)[/math] y el otro asociado a la parte temporal [math]T(t)[/math]. El primero de ellos será:

[math] \begin{cases} X^{’’}(x)+ \lambda \cdot X(x)~~~ \text{ en } x \in [0,1]\\ X(0)=X(1)=0\\ \end{cases}[/math]

Y el segundo problema será:

[math] \begin{cases} T^{’’}(t)+ \lambda \cdot T(t)~~~ \text{ con } t\gt0\\ \end{cases}[/math]

Teniendo en cuenta los posibles signos del valor de [math] \lambda [/math] se llega a que la solución general de [math]u(x,t)[/math] es de la forma:

[math] u(x,t)=\sum_{k=1}^{\infty} c_k \cdot sen(k \pi x) \cdot cos(k \pi t) + d_k \cdot sen(k \pi x) \cdot sen(k \pi t) [/math]

Donde [math] c_k, d_k [/math] son los coeficientes de la serie de Fourier que calcularemos imponiendo las condiciones iniciales.

Ahora, suponemos que tomamos como datos iniciales [math] u_{0}(x) = e^{ −100(x−1/2)^{2}} [/math] y [math] u_{1}(x) = 0 [/math]. Estos datos iniciales los incorporamos a nuestro programa de Matlab y dibujamos la solución en el intervalo de tiempo [math] t ∈ [0, 2] [/math] tomando los primeros [math] 50 [/math] términos de la serie. 

clear all
close all
 
% CONDICIONES INICIALES
u_0=@(x) exp(-100.*(x-1/2).^2);
u_1=@(x) 0;
 
% TÉRMINO K-ÉSIMO DE LA SOLUCIÓN DEL PROBLEMA HABIENDO RESTADO LA SOLUCIÓN ESTACIONARIA
u_k=@(x,t,k,coef_fourier_c,coef_fourier_d) coef_fourier_c.*sin(k.*pi.*x).*cos(k.*pi.*t)+coef_fourier_d.*sin(k.*pi.*x).*sin(k.*pi.*t); 
 
% DATOS PARA LA DISCRETIZACIÓN DEL ESPACIO
extremo_izquierdo=0;
extremo_derecho=1;
paso_espacio=0.001;
 
% DATOS PARA LA DISCRETIZACIÓN DEL TIEMPO
tiempo_inicial=0;
tiempo_final=2;
paso_tiempo=0.001;
 
% DISCRETIZACIÓN DEL ESPACIO Y DEL TIEMPO
[X,T]=meshgrid(extremo_izquierdo:paso_espacio:extremo_derecho,tiempo_inicial:paso_tiempo:tiempo_final); % matrices de mallado espacio-temporal
 
u=zeros(size(X)); % inicializamos la variable w (solución del problema habiendo restado la solución estacionaria)
 
% SOLUCIÓN EN EL INTERVALO DE TIEMPO [0,1] TOMANDO LOS k_end PRIMEROS TÉRMINOS DE LA SERIE 
k_end=50; % último valor de k del sumatorio, siendo k el índice del sumatorio de la solución
 
for ks=1:k_end % queremos resolver el problema para distintos valores de k
    c_k=trapz(X(1,:),sin(ks.*pi.*X(1,:)).*u_0(X(1,:)))./trapz(X(1,:),(sin(ks.*pi.*X(1,:))).^2);
    d_k=trapz(X(1,:),sin(ks.*pi.*X(1,:)).*u_1(X(1,:)))./trapz(X(1,:),(sin(ks.*pi.*X(1,:))).^2);
    u=u+u_k(X,T,ks,c_k,d_k);
end
grafica=figure(1)
mesh(X,T,u)
colorbar
colormap('jet')
xlabel('Espacio')
ylabel('Tiempo')
foto=getframe(grafica);
title('Representación de la solución')
imwrite(foto.cdata,"Ej1_Ap3.png")
 
 
% Crear un vídeo
pelicula=VideoWriter('Ej1_Ap3_Cuerda','MPEG-4'); %creo el video
pelicula.FrameRate=50; %controla velocidad
open(pelicula); %abrimos el vídeo
vector_tiempo=T(:,1)
for t=1:2:length(vector_tiempo)
    figura=figure(2)
    plot(X(1,:),u(t,:),'Linewidth',1.5)
    title(['Cuerda en t = ' num2str(vector_tiempo(t))]);
    axis([0 1 -1 1]);
 
    % Añadir el frame al vídeo
    imagen=getframe(figura);
    writeVideo(pelicula,imagen); %inserto la imagen
   
end
 
% Cerrar el vídeo
close(pelicula);


Representación de la solución obtenida para [math] u_{0}(x) = e^{ −100(x−1/2)^{2}} [/math] y [math] u_{1}(x) = 0 [/math] en el intervalo de tiempo [math] t \in [0,2] [/math] tomando los [math]50[/math] primeros términos de la serie.
Representación de la evolución de la solución obtenida para [math] u_{0}(x) = e^{ −100(x−1/2)^{2}} [/math] y [math] u_{1}(x) = 0 [/math] en el intervalo de tiempo [math] t \in [0,2] [/math] tomando los [math]50[/math] primeros términos de la serie.

Para el código, hemos comenzado definiendo las condiciones iniciales de la onda “u_0” y “u_1”. Luego, hemos definido una función “u_k”que representa el término k-ésimo de la solución del problema. A continuación, se discretiza el espacio y el tiempo para generar una malla espacio-temporal. Se inicializa la solución u y se calcula iterativamente la contribución de los primeros k términos a la solución. Después de obtener la solución en toda la malla, se grafica en forma de superficie tridimensional. Finalmente, se crea un video que muestra la evolución de la onda en el tiempo.

Podemos ver como la solución es periódica de tiempo [math] T=2 [/math] . Al inicio, se tiene una onda con una oscilación de valor máximo [math]1 [/math] en [math]x=0.5[/math].

Con el paso del tiempo, la onda se divide en dos y cada una se va trasladando a uno de los extremos de la cuerda. Cuando las ondas alcanzan los extremos se trasladan al eje [math] y [/math] negativo, produciéndose el mismo movimiento de manera simétrica e inversa con respecto a este eje.

Las ondas vuelven a la frontera, se trasladan otra vez al eje [math] y [/math] positivo original y se realiza el mismo movimiento transformándose ambas ondas en la onda original. Vuelve al estado inicial en tiempo [math] t=2 [/math].

Además, hemos representado el movimiento de la cuerda en [math] 3D [/math]. En el eje horizontal de la gráfica se representa tanto la posición espacial como temporal de la onda. Vemos como en tiempo [math]t=0[/math] la cuerda se encuentra en la posición inicial del video. A medida que pasa el tiempo, su posición máxima de [math]1[/math] unidad en el centro del espacio se va trasladando conservando siempre el valor [math]0[/math] en los extremos de la cuerda.


3.1 Cambio en los valores iniciales: Una onda que viaja en un sentido

Suponemos ahora que tomamos como datos iniciales los correspondientes a una onda que viaja en un solo sentido, es decir, la solución [math]u(x,t)[/math] es de la forma [math] u(x,t)=f(x-t)[/math].

Para ello, tomamos [math] u(x,0)=u_0(x)=f(x) [/math] y [math] u_t(x,0)=u_1(x)=-f’(x) [/math], siendo [math] u_0(x)[/math] como en el caso anterior, es decir, [math] u_0(x)=f(x)= e^{ −100(x−1/2)^{2}}[/math].

Implementando la solución en Matlab hemos obtenido los siguientes resultados: 

clear all
close all
 
% CONDICIONES INICIALES
u_0=@(x) exp(-100.*(x-1/2).^2);
u_1=@(x) 200.*(x-1/2).*exp(-100.*(x-1/2).^2);
 
% TÉRMINO K-ÉSIMO DE LA SOLUCIÓN DEL PROBLEMA HABIENDO RESTADO LA SOLUCIÓN ESTACIONARIA
u_k=@(x,t,k,coef_fourier_c,coef_fourier_d) coef_fourier_c.*sin(k.*pi.*x).*cos(k.*pi.*t)+coef_fourier_d/(k.*pi).*sin(k.*pi.*x).*sin(k.*pi.*t); 
 
% DATOS PARA LA DISCRETIZACIÓN DEL ESPACIO
extremo_izquierdo=0;
extremo_derecho=1;
paso_espacio=0.001;
 
% DATOS PARA LA DISCRETIZACIÓN DEL TIEMPO
tiempo_inicial=0;
tiempo_final=2;
paso_tiempo=0.001;
 
% DISCRETIZACIÓN DEL ESPACIO Y DEL TIEMPO
[X,T]=meshgrid(extremo_izquierdo:paso_espacio:extremo_derecho,tiempo_inicial:paso_tiempo:tiempo_final); % matrices de mallado espacio-temporal
 
u=zeros(size(X)); % inicializamos la variable w (solución del problema habiendo restado la solución estacionaria)
 
% SOLUCIÓN EN EL INTERVALO DE TIEMPO [0,1] TOMANDO LOS k_end PRIMEROS TÉRMINOS DE LA SERIE 
k_end=50; % último valor de k del sumatorio, siendo k el índice del sumatorio de la solución
 
for ks=1:k_end % queremos resolver el problema para distintos valores de k
    c_k=trapz(X(1,:),sin(ks.*pi.*X(1,:)).*u_0(X(1,:)))./trapz(X(1,:),(sin(ks.*pi.*X(1,:))).^2);
    d_k=trapz(X(1,:),sin(ks.*pi.*X(1,:)).*u_1(X(1,:)))./trapz(X(1,:),(sin(ks.*pi.*X(1,:))).^2);
    u=u+u_k(X,T,ks,c_k,d_k);
end
grafica=figure(1)
mesh(X,T,u)
colorbar
colormap('jet')
xlabel('Espacio')
ylabel('Tiempo')
foto=getframe(grafica);
title('Representación de la solución')
imwrite(foto.cdata,"Ej1_Ap4.png")
 
 
% Crear un vídeo
pelicula=VideoWriter('Ej1_Ap4_Cuerda','MPEG-4'); %creo el video
pelicula.FrameRate=50; %controla velocidad
open(pelicula); %abrimos el vídeo
vector_tiempo=T(:,1)
for t=1:2:length(vector_tiempo)
    figura=figure(2)
    plot(X(1,:),u(t,:),'Linewidth',1.5)
    title(['Cuerda en t = ' num2str(vector_tiempo(t))]);
    axis([0 1 -8 8]);
 
    % Añadir el frame al vídeo
    imagen=getframe(figura);
    writeVideo(pelicula,imagen); %inserto la imagen
   
end
 
% Cerrar el vídeo
close(pelicula);


Representación de la solución obtenida para [math] u_{0}(x) = e^{ −100(x−1/2)^{2}} [/math] suponiendo que la onda viaja en un solo sentido.
Representación de la evolución de la solución obtenida para [math] u_{0}(x) = e^{ −100(x−1/2)^{2}} [/math] suponiendo que la onda viaja en un solo sentido

En la gráfica, la solución muestra un comportamiento oscilatorio, típico de las ondas, con los puntos máximos y mínimos moviéndose a través del dominio espacial a medida que el tiempo avanza. Además, se puede apreciar como la solución viaja a velocidad constante con velocidad [math] c=1 ~ m/s[/math]. Por otro lado, se puede observar como cuando la onda llega a la frontera rebota y se comienza a propagar con sentido y amplitud contrarios. Esto ocurre porque para mantener [math]u(0,t)=u(1,t)=0[/math] , la única forma en que la onda puede seguir cumpliendo estas condiciones después de la reflexión es si la parte reflejada de la onda tiene la amplitud opuesta. Se produce ese comportamiento de manera constante en el tiempo. Hay que destacar que, en tiempo [math] t=2 [/math] la onda ha regresado a su posición original.

3.2 Cambio en condiciones frontera: Incorporación de condiciones de frontera tipo Neumann

En este último caso, vamos a realizar un cambio en las condiciones frontera del sistema. Incorporamos condiciones fronteras tipo Neumann para sustituir a las condiciones tipo Dirichlet presentes en el problema original. De esta manera, tendremos el siguiente problema:

[math] \begin{cases} u_{tt}-u_{xx} = 0~~~ \text{ con } x \in [0,1] \text{ y } t\gt0\\ u_x(0,t)=u_x(1,t)=0 \\ u(x,0)=u_0(x) \\ u_t(x,0)=u_1(x) \\ \end{cases}[/math]

Implementamos la solución en Matlab y obtenemos el siguiente resultado: 

clear all
close all
 
% CONDICIONES NEUMANN
u_0=@(x) exp(-100.*(x-1/2).^2);
u_1=@(x) 200.*(x-1/2).*exp(-100.*(x-1/2).^2);
 
% TÉRMINO K-ÉSIMO DE LA SOLUCIÓN DEL PROBLEMA
u_k=@(x,t,k,coef_fourier_c,coef_fourier_d) coef_fourier_c.*cos(k.*pi.*x).*cos(k.*pi.*t)+coef_fourier_d/(k.*pi).*cos(k.*pi.*x).*sin(k.*pi.*t); 
 
% DATOS PARA LA DISCRETIZACIÓN DEL ESPACIO
extremo_izquierdo=0;
extremo_derecho=1;
paso_espacio=0.001;
 
% DATOS PARA LA DISCRETIZACIÓN DEL TIEMPO
tiempo_inicial=0;
tiempo_final=2;
paso_tiempo=0.001;
 
% DISCRETIZACIÓN DEL ESPACIO Y DEL TIEMPO
[X,T]=meshgrid(extremo_izquierdo:paso_espacio:extremo_derecho,tiempo_inicial:paso_tiempo:tiempo_final); % matrices de mallado espacio-temporal
 
u=zeros(size(X)); % inicializamos la variable w (solución del problema habiendo restado la solución estacionaria)
 
% SOLUCIÓN EN EL INTERVALO DE TIEMPO [0,1] TOMANDO LOS k_end PRIMEROS TÉRMINOS DE LA SERIE 
k_end=50; % último valor de k del sumatorio, siendo k el índice del sumatorio de la solución
 
% Añadimos a la solución el término cuando k=0
c_0=trapz(X(1,:),cos(0.*pi.*X(1,:)).*u_0(X(1,:)))./trapz(X(1,:),(cos(0.*pi.*X(1,:))).^2);
u=u+c_0.*cos(0.*pi.*X).*cos(0.*pi.*T);
 
for ks=1:k_end % queremos resolver el problema para distintos valores de k
    c_k=trapz(X(1,:),cos(ks.*pi.*X(1,:)).*u_0(X(1,:)))./trapz(X(1,:),(cos(ks.*pi.*X(1,:))).^2);
    d_k=trapz(X(1,:),cos(ks.*pi.*X(1,:)).*u_1(X(1,:)))./trapz(X(1,:),(cos(ks.*pi.*X(1,:))).^2);
    u=u+u_k(X,T,ks,c_k,d_k);
end
grafica=figure(1)
mesh(X,T,u)
colorbar
colormap('jet')
xlabel('Espacio')
ylabel('Tiempo')
title('Representación de la solución')
foto=getframe(grafica);
imwrite(foto.cdata,"Ej1_Ap5.png")
 
 
% Crear un vídeo
pelicula=VideoWriter('Ej1_Ap5_Cuerda','MPEG-4'); %creo el video
pelicula.FrameRate=50; %controla velocidad
open(pelicula); %abrimos el vídeo
vector_tiempo=T(:,1)
for t=1:2:length(vector_tiempo)
    figura=figure(2)
    plot(X(1,:),u(t,:),'Linewidth',1.5)
    title(['Cuerda en t = ' num2str(vector_tiempo(t))]);
    axis([0 1 -8 8]);
 
    % Añadir el frame al vídeo
    imagen=getframe(figura);
    writeVideo(pelicula,imagen); %inserto la imagen
   
end
 
% Cerrar el vídeo
close(pelicula);
Representación de la solución obtenida
Representación de la evolución de la solución obtenida


En este código, se sigue el mismo algoritmo que el anterior, pero cambiando las condiciones de contorno de Neumann en lugar de condiciones de contorno de Dirichlet.

El término k-ésimo de la solución “u_k” se redefine para adaptarse a estas nuevas condiciones de contorno. Nuestra base de Fourier para el problema anterior estaba formada por el conjunto [math] \left\{sin(k \pi x)\right\}[/math]. Sin embargo, imponiendo las condiciones de Neuman nuestra nueva base es [math]\left\{1,cos(k \pi x) \right\} [/math]. Se muestra gráficamente y se genera un video de la evolución temporal de la onda.

En la gráfica, la solución muestra un comportamiento oscilatorio, esta vez solo de valores máximos moviéndose a través del dominio espacial a medida que el tiempo avanza. La condición nueva [math]u_x(0,t) = u_x(1,t) = 0 [/math] indica que la derivada espacial de u es cero en los extremos.

Cabe destacar que no se produce el fenómeno de reflexión con respecto al eje x una vez la onda alcanza el extremo izquierdo de la cuerda, tal y como sucedía en el caso anterior. Esto se debe a que ya no tenemos las condiciones de Dirichlet que fijaban la condición frontera para el valor 0, en este nuevo caso la condición frontera puedo tomar valores distintos de 0.

Este comportamiento, lo podemos apreciar también en el vídeo, donde vemos que la onda se deplaza y alcanzan su máximo valor cuando dicho desplazamiento alcanza los extremos de la cuerda. En conreto, la solución toma un valor de [math]2[/math].

4 Solución fundamental de ondas

La solución fundamental de ondas es la solución que se obtiene al dar un impulso inicial muy localizado en [math] x=0 [/math]. A lo largo de esta subsección analizaremos la solución fundamental de la ecuación de ondas en dimensiones [math]1, 2[/math] y [math]3[/math].

Matemáticamente, la solución fundamental es la solución del sistema relativo a la ecuación de ondas:

[math]\left \{ \begin{array}{ll} u_{tt}-c^2\Delta u=0, \quad x \in \mathbb{R}^n, t\gt0, \\ u(x,0)=0, u_t(x,0)=\delta (x), \quad x \in \mathbb{R}^n, \end{array} \right. [/math]


donde [math] \delta(x) [/math] es la delta de Dirac. Formalmente, se define como el siguiente límite:

[math] \delta (x) = \lim \limits_{r \rightarrow 0} \frac{1}{|B(0,r)|} \chi_{B(0,r)}(x) \sim \left \{ \begin{array}{ll} \infty \quad x=0, \\ 0 \quad x\neq 0, \end{array} \right. [/math]

donde [math]\chi_{B(0,r)}(x) [/math] es la función característica de la bola centrada en [math] 0 [/math] de radio [math] r [/math] y [math]|B(0,r)|[/math] es el volumen de la bola.

Nos centraremos ahora en la solución fundamental en dimensiones [math] 1,2 [/math] y [math] 3 [/math], dado que son las que Podemos representar con facilidad. Sus expresiones asociadas son:

En dimensión [math]n=1[/math]:

[math] K_1(x,t) = \frac{1}{2c}[H(x+ct)-H(x-ct)], [/math]

Donde [math]H(s)[/math] es la función de Heaviside,

[math] H(s) = \left \{ \begin{array}{ll} 0, \quad x\lt0 \\ 1, \quad x \geq 0 \end{array} \right. [/math]

En dimensión [math] n=2 [/math], la expresión de la solución fundamental es:

[math] K_2(x,t)=\frac{1}{2\pi c \sqrt{c^2t^2-|x|^2}} \chi_{B(0,ct)}(x), [/math]

donde [math] \chi_{B(0,ct)}(x) [/math] es la función característica de la bola de centro [math] 0 [/math] y radio [math] ct [/math].

Por último, en dimensión [math] 3 [/math], la expresión es:

[math] K_3(x,t)=\frac{\delta (|x|-ct)}{4\pi c|x|}, t\gt0. [/math]

4.1 Solución fundamental de ondas en una dimensión

En primer lugar, estudiamos la solución fundamental de ondas en una dimensión. Recordemos que dicha expresión era:

[math] K_1(x,t) = \frac{1}{2c}[H(x+ct)-H(x-ct)], [/math]

Donde [math]H(s)[/math] es la función de Heaviside,

[math] H(s) = \left \{ \begin{array}{ll} 0, \quad x\lt0 \\ 1, \quad x \geq 0 \end{array} \right. [/math]

Gráficamente la solución es la siguiente:

Representación de la solución fundamental de ondas en dimensión [math] n=1 [/math]
clear all
close all
 
c=1;
 
% Definimos la solución fundamental para n=1
k_1=@(x,t) 1/(2.*c).*((x+c.*t>=0)-(x-c.*t>=0));
 
% Datos inicial y final del espacio
espacio_inicial=-100;
espacio_final=100;
 
% Datos inicial y final del tiempo
tiempo_inicial=0;
tiempo_final=100;
 
% Discretización del espacio y del tiempo
espacio=linspace(espacio_inicial,espacio_final,1000);
tiempo=linspace(tiempo_inicial,tiempo_final,1000);
 
% Mallado
[Espacio,Tiempo]=meshgrid(espacio,tiempo);
 
% Representamos la solución
grafica=figure(1)
mesh(Espacio,Tiempo,k_1(Espacio,Tiempo))
colorbar
colormap('jet')
xlabel('Espacio')
ylabel('Tiempo')
title('Representación de la solución fundamental para n=1')
foto=getframe(grafica);
imwrite(foto.cdata,"Ej2_Ap1_n1.png")


Podemos observar como la solución obtenida es radial, es decir, no depende de [math] x [/math], únicamente del valor absoluto de x, [math] |x| [/math]. Esto puede verse fácilmente observando la paridad de la función en la representación. Analíticamente vemos como:

[math] K_1(-x,t) = \frac{1}{2c}[H(-x+ct)-H(-x-ct)] = \frac{1}{2c}[H(-(x-ct))-H(-(x+ct))] = \frac{1}{2c}[-H(x-ct)+H(x+ct)] = \frac{1}{2c}[H(x+ct)-H(x-ct)] = K_1(x,t)[/math]

Donde en la tercera igualdad se ha usado que, según la definición de la función de Heaviside, [math] H(-s)= 1-H(s) [/math].

Con respecto al código, primero se define la solución fundamental en una dimensión además de los datos necesarios: puntos del espacio y tiempo iniciales y finales. Con ello, se define una discretización espacio-temporal y su correspondiente mallado. Por último, representamos la solución gráficamente e introducimos un comando para que nos guarde directamente la foto en el ordenador.

4.2 Solución fundamental de ondas en dos dimensiones

En segundo lugar, estudiamos la solución fundamental de ondas en dos dimensiones. Su expresión es:

[math] K_2(x,t)=\frac{1}{2\pi c \sqrt{c^2t^2-|x|^2}} \chi_{B(0,ct)}(x), [/math]

Podemos ver como la solución es radial ya que solo depende de la norma de [math]x[/math]. Esta solución vista como solución radial simplifica considerablemente el problema, en particular permite la representación de la solución con el transcurso del tiempo.

Si implementamos nuestra solución en Matlab, tenemos que añadir la condición de [math]c=\epsilon=0.01[/math] para evitar la singularidad.

Representación de la solución fundamental de ondas en dimensión [math] n=2 [/math]
clear all
close all
 
c=1;
 
% Definimos la solución fundamental para n=2
k_2=@(r,t) (r<=c*t)./(0.01+2*pi*c.*sqrt(c^2*t.^2-r.^2));
 
% Datos inicial y final del espacio
radio_inicial=0;
radio_final=1;
 
% Datos inicial y final del tiempo
tiempo_inicial=0;
tiempo_final=1;
 
% Discretización del espacio y del tiempo
radio=linspace(radio_inicial,radio_final,1000);
tiempo=linspace(tiempo_inicial,tiempo_final,1000);
 
% Mallado
[Radio,Tiempo]=meshgrid(radio,tiempo);
 
% Representamos la solución
grafica=figure(1)
mesh(Radio,Tiempo,k_2(Radio,Tiempo))
colorbar
colormap('jet')
xlabel('Radio')
ylabel('Tiempo')
title('Representación de la solución fundamental para n=2')
foto=getframe(grafica);
imwrite(foto.cdata,"Ej2_Ap1_n2.png")


En esta imagen, apreciamos el comportamiento de la onda, el cual se inicia en el extremo izquierdo con un valor de 100 en [math](x,t)=(0,0)[/math]. Sin embargo, en cuanto aumenta el radio, el valor de la función cambia abruptamente a cero. Con el paso del tiempo, este valor máximo se desplaza hacia el lado derecho hasta el tiempo [math]t=1[/math] donde finalmente alcanza el extremo derecho.

En relación con el código empleado, es muy similar al de para una dimensión. La única diferencia radica en la definición de la solución fundamental. En este caso, se define en ella la función característica en la bola como [math]r[/math]=[math]c[/math]*[math]t[/math] y se define además la regularización de la solución fundamental para evitar la singularidad en su representación.


4.3 Solución fundamental de ondas en tres dimensiones

En este último caso, representamos la solución fundamental de ondas en tres dimensiones. Nuevamente, recordamos que su expresión era:

[math] K_3(x,t)=\frac{\delta (|x|-ct)}{4\pi c|x|}, t\gt0 [/math]

Para representar gráficamente la solución, hemos tenido que optar por emplear la aproximación de la delta de Dirac. Esto es dado que la delta de Dirac presenta una singularidad. Tomamos entonces:

[math] \delta (s) \sim \phi _k (s) = \sqrt{\frac{k}{\pi}}e^{-ks^2}, ~~~~~ k\gt\gt1, [/math]

Concretamente, hemos tomado [math]k=1000[/math]. Además, podemos observar como se verifica que [math] \int_{-\infty}^{+\infty} \phi _k (s) ~ \text{ds}=1[/math], criterio que debe verificar la delta de Dirac.

La solución obtenida en este caso es la siguiente:

Representación de la solución fundamental de ondas en dimensión [math] n=3 [/math]
clear all
close all
 
c=1;
 
% Definimos la aproximación de la delta de Dirac
aprox_dirac=@(k,s) sqrt(k/pi)*exp(-k*s.^2);
 
% Definimos la solución fundamental para n=3
k_3=@(r,t) aprox_dirac(1000,r-c.*t)./(4*pi*c.*r);
 
% Datos inicial y final del espacio
radio_inicial=0;
radio_final=1;
 
% Datos inicial y final del tiempo
tiempo_inicial=0;
tiempo_final=1;
 
% Discretización del espacio y del tiempo
radio=linspace(0,1,1000);
tiempo=linspace(0,1,1000);
 
% Mallado
[Radio,Tiempo]=meshgrid(radio,tiempo);
 
% Representamos la solución
grafica=figure(1)
mesh(Radio,Tiempo,k_3(Radio,Tiempo))
colorbar
colormap('jet')
xlabel('Radio')
ylabel('Tiempo')
title('Representación de la solución fundamental para n=3')
foto=getframe(grafica);
imwrite(foto.cdata,"Ej2_Ap1_n3.png")


Análogamente, la solución fundamental es radial, dependiendo únicamente de la norma de [math] \textbf{x} [/math].

Por otro lado, podemos observar como la solución fundamental tiene una singularidad en el punto [math](\textbf{x},t)=(\textbf{0},0) [/math] debido a la definición de la delta de Dirac. Al haber tomado la aproximación de la delta de Dirac, esta singularidad no se produce gráficamente, aunque sí que toma valores bastante grandes con respecto a los demás valores espacio-temporales del dominio representado.

Con respecto al código, de nuevo, el funcionamiento es análogo a los dos anteriores. En primer lugar, se define la aproximación de la delta de Dirac y la solución fundamental en dimensión 3, además de los datos necesarios: puntos del espacio y tiempo iniciales y finales. Con ello, se define una discretización espacio-temporal y el mallado a emplear. Por último, representamos la solución gráficamente.

4.4 Solución del sistema usando la solución fundamental para n=2

En esta última sección, se busca calcular y representar la solución del problema en dimensión 2

[math]\left \{ \begin{array}{ll} u_{tt}-c^2\Delta u=0, \quad x \in \mathbb{R}^n, t\gt0, \\ u(x,0)=0, u_t(x,0)=h (x)=\chi_{B(0,1/2)}(x) , \quad x \in \mathbb{R}^n, \end{array} \right. [/math]

que viene dada por la convolución [math] u(x,t)=\int_{\mathbb{R^2}} K_2(x-y,t) h(y) dy [/math] .

En primer lugar, es importante destacar que dicha solución es radial en x pues tanto la solución fundamental en dos dimensiones como la función [math] h(x)=\chi_{B(0,1/2)}(x) [/math] son funciones radiales y, con ello, la convolución también será una función radial.

Por lo tanto, si suponemos que [math] x=(r_1 \cdot cos(\theta_1), r_1 \cdot sin(\theta_1)) [/math] al usar coordenadas polares, esta propiedad nos permite que a la hora de calcular la solución del problema planteado, basta con calcular la solución cuando [math] \theta_1=0 [/math] y posteriormente tomar el mismo valor de la solución para cualquier valor [math] \theta_1 \in [0,2\pi] [/math].

Retomando la integral con la que definimos la solución, vamos a calcularla usando coordenadas polares, de manera que sean [math] x=(r_1 cos(0), r_1 sin(0)) [/math] e [math] y= (r_2 cos(\theta_2),r_2 sin(\theta_2)) [/math], podemos expresar [math] |x-y|^2=r_1^2+r_2^2-2 r_1 r_2 cos(\theta_2) [/math]. Con ello, la solución del problema viene dada en coordenadas polares por la siguiente expresión:

[math] u(r_1,t)=\int_{0}^{2 \pi} \int_{0}^{1/2} \frac{1_{B(0,ct)} (r_1 cos(\theta_1)-r_2 cos(\theta_2),-r_2 sin(\theta_2))}{2 \pi c \sqrt{c^2 t^2 – (r_1^2+r_2^2-2 r_1 r_2 cos(\theta_2))}} dr_2 d\theta_2 [/math]

Teniendo esto en cuenta, hemos elaborado el siguiente código de Matlab que representa la solución para los tiempos [math] t=0,0.5,1,2 [/math] y velocidad de propagación [math] c=1 [/math], sin olvidar que hay que deshacer los cambios de variables realizados. 

clear all
close all
 
warning('off','all')
 
c=1;
 
% Datos inicial y final del espacio
radio_inicial=0;
radio_final=1;
 
% Datos inicial y final de theta
theta_inicial=0;
theta_final=2*pi;
 
% Tiempos a estudiar
tiempos=[0,0.5,1,2];
 
% Discretización del espacio y del tiempo
theta_1=linspace(theta_inicial,theta_final,1000);
radio=linspace(radio_inicial,radio_final,1000);
 
% Mallado
[Theta,Radio]=meshgrid(theta_1,radio);
 
u=zeros(size(Theta));
u_eps=zeros(size(Theta));
 
k_2=@(r_1,r_2,theta_2,t) (sqrt(r_1.^2+r_2.^2-2.*r_1.*r_2.*cos(theta_2))<=c*t)./(2*pi*c.*sqrt(abs(c^2*t.^2-r_1.^2-r_2.^2+2.*r_1.*r_2.*cos(theta_2))));
k_2_eps=@(r_1,r_2,theta_2,t) (sqrt(r_1.^2+r_2.^2-2.*r_1.*r_2.*cos(theta_2))<=c*t)./(0.01*(c^2*t.^2-r_1.^2-r_2.^2+2.*r_1.*r_2.*cos(theta_2)<eps)+2*pi*c.*sqrt(abs(c^2*t.^2-r_1.^2-r_2.^2+2.*r_1.*r_2.*cos(theta_2))));
 
for j=1:length(tiempos)
    tiempo=tiempos(j);
    for i=1:length(radio)
        u(i,:)=integral2(@(r_2,theta_2)k_2(radio(i),r_2,theta_2,tiempo), 0,1/2, 0,2*pi)*ones(1,length(theta_1)); 
        u_eps(i,:)=integral2(@(r_2,theta_2)k_2_eps(radio(i),r_2,theta_2,tiempo), 0,1/2, 0,2*pi)*ones(1,length(theta_1)); 
    end
 
    % Representación de la solución deshaciendo el cambio de variable
    figura=figure(j)
 
    set(figura,'position',[2,2,1400,550])
    subplot(1,2,1)
    mesh(Radio.*cos(Theta),Radio.*sin(Theta),u)
    colorbar
    colormap('jet')
    xlabel('x_1')
    ylabel('x_2')
    title("Usando la solución fundamental para n=2")
 
    subplot(1,2,2)
    mesh(Radio.*cos(Theta),Radio.*sin(Theta),u_eps)
    colorbar
    colormap('jet')
    xlabel('x_1')
    ylabel('x_2')
    title("Usando la regularización para n=2")
 
    sgtitle("Solución para t="+num2str(tiempo))
 
    foto=getframe(figura);
    cadena="Ej2_Ap2_t_"+num2str(tiempo)+".png"
    imwrite(foto.cdata,cadena)
end


Inicialmente, se definen los parámetros del problema, como la velocidad de propagación de la onda “c” y los límites del espacio y el ángulo. Luego, se discretiza el espacio y el tiempo, creando una malla.

La función calcula el término k-ésimo de la solución del problema para cada punto en la malla, mientras que “k_2_eps” realiza lo mismo pero con una regularización para evitar singularidades numéricas. Estos términos se integran numéricamente sobre la región de integración correspondiente a cada punto en el espacio y el ángulo, generando así la solución en cada tiempo dado.

La solución se visualiza en dos subgráficos, uno utilizando la solución fundamental y otro utilizando la regularización para comparar los resultados. Este proceso se repite para cada tiempo especificado en el vector "tiempos".

Cabe destacar que este código requiere bastante tiempo de ejecución, pues ha de calcular muchas integrales. Por ese mismo motivo, hemos representado en un primer lugar la solución tomando una discretización para los valores de [math] \theta_1 [/math] de 100 puntos y para los valores de [math] r_1 [/math] también de 100 puntos. Con ello, hemos obtenido los siguientes resultados.

Representación de la solución en tiempo t=0
Representación de la solución en tiempo t=0.5
Representación de la solución en tiempo t=1
Representación de la solución en tiempo t=2


En primer lugar, para entender estas gráficas es imprescindible destacar el fenómeno que sucede en la representación para [math] t=0.5 [/math]. Para ello, debemos recordar que la solución fundamental en dos dimensiones presentaba una singularidad. Ese es el motivo por el que en la primera gráfica hay una franja en la que aparentemente la solución no está definida. De hecho, lo que está sucediendo en concreto es que para dichos puntos Matlab devuelve el valor NaN pues el denominador de la expresión se anula.

Para solventar este inconveniente en la representación, hemos optado en apoyarnos en la regularización ya mencionada anteriormente en este trabajo. Sin embargo, no hemos aplicado dicha expresión exactamente, pues hemos optado por sumar [math] \epsilon=0.1 [/math] cuando el denominador se anula. Esta modificación realizada es la que se puede apreciar en la segunda gráfica de cada una de las imágenes.

Con respecto a la primera de las imágenes, en la que se representa la solución para [math] t=0 [/math], esta nos permite confirmar que la solución cumple la condición inicial impuesta en el problema.

Finalmente, hemos representado la solución tomando una discretización para los valores de [math] \theta_1 [/math] de 1000 puntos y para los valores de [math] r_1 [/math] también de 1000 puntos. Los resultados obtenidos han sido los siguientes:

Representación de la solución en tiempo t=0
Representación de la solución en tiempo t=0.5
Representación de la solución en tiempo t=1
Representación de la solución en tiempo t=2

Como se puede apreciar, la singularidad que se producía en la representación anterior ya no es apreciable al tomar esta nueva discretización, pues la malla de puntos es más fina. Realmente, la franja de puntos que no se representan en la gráfica sigue existiendo, pero al ser la malla mucho más fina, dicha región pasa a ser inapreciable.


5 Referencias

Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Ecuación_de_Ondas,_(equipo_LUA)&oldid=73110»