Diferencia entre revisiones de «Ecuación de ondas (GRwM)»
(→Cambio a condición frontera de tipo Neumann) |
(→Resolución del sistema) |
||
| (No se muestran 45 ediciones intermedias de 3 usuarios) | |||
| Línea 15: | Línea 15: | ||
= Ecuación de ondas I= | = Ecuación de ondas I= | ||
En esta sección, vamos a representar diferentes soluciones de la ecuación de ondas en una dimensión. | En esta sección, vamos a representar diferentes soluciones de la ecuación de ondas en una dimensión. | ||
| − | Para ello, vamos a considerar una cuerda vibrante en el intervalo <math> [0,1]</math>, con densidad <math> d</math> y tensión <math> \tau_0 </math> constante. De modo que la velocidad de propagación <math> c=1 m/s</math> . Además, vamos a suponer que la cuerda está fija en los extremos y vamos a denotar <math> u_0(x)</math> y <math> u_1(x)</math> su posición e impulso iniciales respectivamente. | + | Para ello, vamos a considerar una cuerda vibrante en el intervalo <math> [0,1]</math>, con densidad <math> d</math> y tensión <math> \tau_0 </math> constante. De modo que la velocidad de propagación es <math> c=1 ~~m/s</math> . Además, vamos a suponer que la cuerda está fija en los extremos y vamos a denotar por <math> u_0(x)</math> y <math> u_1(x)</math> su posición e impulso iniciales respectivamente. |
==Planteamiento del problema == | ==Planteamiento del problema == | ||
| − | Comenzamos | + | Comenzamos planteando el sistema de EDP que modeliza el comportamiento de los desplazamientos transversales de la cuerda. Para ello, basta considerar la ecuación de ondas <math> u_{tt} – c^2u_{xx}=0</math> y añadir las condiciones iniciales y de frontera descritas en el apartado anterior. De esta manera se obtiene que el sistema EDP a resolver es: |
<center><math>\left \{ \begin{array}{ll} | <center><math>\left \{ \begin{array}{ll} | ||
| Línea 47: | Línea 47: | ||
| − | <math> d_k= \frac{\int_0^1 sin (k \pi x) u_1 (x) dx}{\int_0 ^1 sin^2(k \pi x) dx}</math> | + | <math> d_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> |
Nótese que al considerar t como un múltiplo de dos, por la estructura de la solución, se obtiene de nuevo la condición inicial. Es decir, la solución es periódica con periodo dos. Este resultado lo analizaremos con más detalle en los siguientes apartados. | Nótese que al considerar t como un múltiplo de dos, por la estructura de la solución, se obtiene de nuevo la condición inicial. Es decir, la solución es periódica con periodo dos. Este resultado lo analizaremos con más detalle en los siguientes apartados. | ||
| Línea 64: | Línea 64: | ||
| − | Como podemos observar | + | Como podemos observar, la amplitud máxima de la onda es <math> 1 </math> y se alcanza para valores enteros del tiempo. Además, para <math>x=0 </math> y <math>x=1</math> se tiene que <math>u(x,t)=0</math>, pues por hipótesis la cuerda está fija en los extremos. |
| − | Además, la solución es periódica en tiempo. Para poder apreciar esto, | + | Además, la solución es periódica en tiempo. Para poder apreciar esto, se presenta el siguiente vídeo. |
| Línea 153: | Línea 153: | ||
| − | [[Archivo: | + | [[Archivo:Ejercicio1apartado4bueno.png|400px|thumb|center| Representación de la solución de la ecuación de ondas con datos iniciales <math> u_0(x)= f(x)</math> y <math> u_1(x)= - f’(x)</math>, con <math>c=1</math>, <math> t \in [0,2] </math> y considerando <math>1000</math> términos de la serie.]] |
| + | |||
| − | |||
| − | |||
| − | Además, se tiene que el perfil de la solución avanza con una velocidad de propagación <math>c= 1 m/s </math>. Para poder apreciar esto, veamos la gráfica anterior desde otra perspectiva: | + | Además, se tiene que el perfil de la solución avanza con una velocidad de propagación <math>c= 1~~ m/s </math>. Para poder apreciar esto, veamos la gráfica anterior desde otra perspectiva: |
| Línea 172: | Línea 171: | ||
| − | La solución es periódica de periodo dos | + | La solución es periódica de periodo dos. Se cumple que cuando la onda llega al extremo, rebota simétricamente, sufriendo un cambio en el signo de la tensión. Además, justo cuando la onda llega a la frontera, se puede observar cómo se anula la solución en todo el intervalo del espacio. |
| Línea 233: | Línea 232: | ||
==Cambio a condición frontera de tipo Neumann== | ==Cambio a condición frontera de tipo Neumann== | ||
| − | Por último, vamos a estudiar cómo varía la solución al considerar condiciones de frontera de tipo Neumann en lugar de tipo Dirichlet. Para ello, en este caso vamos a tomar <math>u_x(0,t) = u_x(1,t) =0</math>. | + | Por último, vamos a estudiar cómo varía la solución al considerar condiciones de frontera de tipo Neumann en lugar de tipo Dirichlet. Para ello, en este caso vamos a tomar <math>u_x(0,t) = u_x(1,t) =0</math>. Entonces, tenemos el siguiente sistema: |
| + | |||
| + | <center><math>\left \{ \begin{array}{ll} | ||
| + | |||
| + | u_{tt} – u_{xx} =0 | ||
| + | |||
| + | \\ | ||
| + | |||
| + | u_x(0,t)=u_x(1,t)=0 | ||
| + | \\ | ||
| + | |||
| + | u(x,0)= e^{-100 (x- \frac{1}{2})^2} | ||
| + | \\ | ||
| + | |||
| + | u_t(x,0)= 200(x- \frac{1}{2})e^{-100 (x- \frac{1}{2})^2} | ||
| + | |||
| + | \end{array} \right. | ||
| + | |||
| + | </math></center> | ||
En este caso, aplicando separación de variables obtenemos que la solución es: | En este caso, aplicando separación de variables obtenemos que la solución es: | ||
| + | <center><math> u(x,t)= \sum _{k=1}^{\infty} c_k cos(k \pi x) cos (k \pi t) + d_k cos (k \pi x) sin(k \pi t) </math></center> | ||
| + | siendo | ||
| + | <math> c_k= \frac{\int_0^1 cos (k \pi x) u_0 (x) dx}{\int_0 ^1 cos^2(k \pi x) dx}</math> | ||
| + | <math> d_k= \frac{\int_0^1 cos (k \pi x) u_1 (x) dx}{\int_0 ^1 cos^2(k \pi x) dx}</math> | ||
| + | Nótese que, al igual que para el problema con condiciones Dirichlet, en este caso la solución también es periódica de periodo dos. | ||
La gráfica de la solución en este caso es: | La gráfica de la solución en este caso es: | ||
| − | [[Archivo: |400px|thumb|center| Representación de | + | |
| + | [[Archivo: Ejercicio1apartado5.png|400px|thumb|center| Representación de la solución de la ecuación de ondas considerando <math>u_x(0,t) = u_x(1,t) =0</math>, con <math>c=1</math> | ||
| + | , <math>t \in [0,2] </math> y considerando <math>1000</math> términos de la serie. ]] | ||
| + | |||
Para poder apreciar mejor el comportamiento de la onda, hemos realizado el siguiente vídeo: | Para poder apreciar mejor el comportamiento de la onda, hemos realizado el siguiente vídeo: | ||
| Línea 253: | Línea 278: | ||
[[Archivo:Vídeo-Ej-1.5.gif|400px|thumb|center| Representación de la solución para <math> t \in [0,4]</math>.]] | [[Archivo:Vídeo-Ej-1.5.gif|400px|thumb|center| Representación de la solución para <math> t \in [0,4]</math>.]] | ||
| − | Como podemos observar, al igual que en los casos anteriores la solución es de periodo <math> 2 | + | |
| − | Además, este caso, como los extremos no están fijos, cuando la onda llega a uno de los extremos, este varía. La razón de esto es que se debe cumplir que <math> u_x=0 </math>. | + | Como podemos observar, al igual que en los casos anteriores, la solución es de periodo <math> 2 </math>. |
| + | Además, en este caso, como los extremos no están fijos, cuando la onda llega a uno de los extremos, este varía. La razón de esto es que se debe cumplir que <math> u_x=0 </math>. | ||
| Línea 314: | Línea 340: | ||
=Ecuación de ondas II= | =Ecuación de ondas II= | ||
| − | |||
| − | |||
| − | + | En esta sección, vamos a estudiar la solución fundamental de la ecuación de ondas en dimensiones 1, 2 y 3 y una aplicación de esta en dimensión 2. Hay que tener en cuenta que la solución fundamental es la solución que se obtiene al dar un impulso inicial muy localizado en <math> x=0 </math>. | |
| + | |||
| + | ==Solución fundamental== | ||
| + | En primer lugar, consideramos el siguiente sistema: | ||
<center><math>\left \{ \begin{array}{ll} | <center><math>\left \{ \begin{array}{ll} | ||
| Línea 326: | Línea 353: | ||
</math></center> | </math></center> | ||
| − | donde <math> \delta(x) </math> es la delta de Dirac | + | donde <math> \delta(x) </math> es la delta de Dirac. |
| − | + | La expresión de la solución fundamental es: | |
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | La expresión de | + | |
<math>1)</math> En dimensión <math>n=1</math>: | <math>1)</math> En dimensión <math>n=1</math>: | ||
| − | <center><math> K_1(x,t) = \frac{1}{2c}[ | + | <center><math> K_1(x,t) =\frac{1}{2c}\textbf{1}_{[-ct,ct]}(x)=\frac{1}{2c}\textbf{1}_{[0,ct]}(|x|), </math></center> |
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | 0, | + | |
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | </math></center> | + | |
| + | Que se puede observar que tiene una expresión analítica radial, luego la solución fundamental para dimensión 1 es radial. | ||
Gráficamente la solución fundamental de la ecuación de ondas en <math> dim =1 </math> es: | Gráficamente la solución fundamental de la ecuación de ondas en <math> dim =1 </math> es: | ||
| Línea 362: | Línea 372: | ||
<center><math> | <center><math> | ||
| − | K_2(x,t)=\frac{1}{2\pi c \sqrt{c^2t^2-|x|^2}} \chi_{B(0,ct)}(x), | + | K_2(x,t)=\frac{1}{2\pi c \sqrt{c^2t^2-|x|^2}} \chi_{B(0,ct)}(x)=\frac{1}{2\pi c \sqrt{c^2t^2-|x|^2}} \textbf{1}_{[0,ct]}(|x|), |
</math></center> | </math></center> | ||
| − | donde <math> \chi_{B(0,ct)}(x) </math> es la función característica de la bola de centro 0 y radio <math> ct </math>. | + | donde <math> \chi_{B(0,ct)}(x) </math> es la función característica de la bola de centro 0 y radio <math> ct </math>. Como se puede observar, también tiene una expresión analítica radial y, por ello, es una solución radial. |
| Línea 385: | Línea 395: | ||
</math></center> | </math></center> | ||
| + | |||
| + | Que ya de por sí tiene también una expresión analítica radial. | ||
Por tanto, la solución fundamental en <math> dim=3</math> se representa de la siguiente forma: | Por tanto, la solución fundamental en <math> dim=3</math> se representa de la siguiente forma: | ||
| Línea 398: | Línea 410: | ||
Con <math>k=1000</math>. | Con <math>k=1000</math>. | ||
| − | == | + | ===Programas=== |
| + | {{matlab|codigo= | ||
| + | clear all | ||
| + | close all | ||
| + | |||
| + | Código de la solución fundamental para dimensión 1. | ||
| + | |||
| + | %%% EJERCICIO 2.1 dim1 %%%% | ||
| + | |||
| + | % Intervalos de definición | ||
| + | x0=-1; xf=1; % Valor inicial y final de x | ||
| + | t0=0; tf=1; % Valor inicial y final de t | ||
| + | difx=10^-2; dift=10^-3; % División del intervalo de x y t respectivamente | ||
| + | X=x0:difx:xf; % Discretización del intervalo de x | ||
| + | T=t0:dift:tf; % Discretización del intervalo de t | ||
| + | |||
| + | |||
| + | [XX,TT]=meshgrid(X,T); | ||
| + | |||
| + | % Velocidad de propagación | ||
| + | c=1; | ||
| + | |||
| + | % Cálculo de la solución fundamental | ||
| + | K=(XX<=c*TT).*(XX>=-c*TT)./(2*c); | ||
| + | |||
| + | surf(XX,TT,K) | ||
| + | shading interp | ||
| + | title("Solución fundamental dim 1") | ||
| + | xlabel("x") | ||
| + | ylabel("t") | ||
| + | |||
| + | }} | ||
| + | |||
| + | Código para representar la solución fundamental para dimensión 2. | ||
| + | |||
| + | {{matlab|codigo= | ||
| + | clear all | ||
| + | close all | ||
| + | |||
| + | %%% EJERCICIO 2.1 dim2 %%%% | ||
| + | |||
| + | |||
| + | % Intervalos de definición | ||
| + | |||
| + | r0=0; rf=1; % Valor inicial y final de r | ||
| + | t0=0; tf=1; % Valor inicial y final de t | ||
| + | divr=10^-2;divt=10^-3; % División del intervalo de r y t respectivamente | ||
| + | R=r0:divr:rf; % Discretización del intervalo de r | ||
| + | T=t0:divt:tf; % Discretización del intervalo de t | ||
| + | |||
| + | [RR,TT]=meshgrid(R,T); | ||
| + | |||
| + | % Velocidad de propagación | ||
| + | c=1; | ||
| + | |||
| + | % Cálculo de la solución fundamental | ||
| + | K=(RR<=c*TT)./(eps + 2*pi*c.*sqrt(c^2*TT.^2-RR.^2)); | ||
| + | |||
| + | surf(RR,TT,K) | ||
| + | shading interp | ||
| + | title("Solución fundamental dim 2 en función del radio") | ||
| + | xlabel("r") | ||
| + | ylabel("t") | ||
| + | |||
| + | |||
| + | }} | ||
| + | |||
| + | |||
| + | Código de la solución fundamental para dimensión 3. | ||
| + | |||
| + | {{matlab|codigo= | ||
| + | clear all | ||
| + | close all | ||
| + | |||
| + | %%% EJERCICIO 2.1 dim3 %%%% | ||
| + | |||
| + | % Intervalos de definición | ||
| + | |||
| + | r0=0; rf=1; % Valor inicial y final de r | ||
| + | t0=0; tf=1; % Valor inicial y final de t | ||
| + | divr=10^-2;divt=10^-3; % División del intervalo de r y t respectivamente | ||
| + | R=r0:divr:rf; % Discretización del intervalo de r | ||
| + | T=t0:divt:tf; % Discretización del intervalo de t | ||
| + | |||
| + | [RR,TT]=meshgrid(R,T); | ||
| + | |||
| + | c=1; % Velocidad de propagación | ||
| + | k=1000; % Parámetro empleado en la aproximación de la delta de Dirac | ||
| + | |||
| + | % Cálculo de la solución fundamental | ||
| + | K=sqrt(k/pi).*exp(-k*(RR-c*TT).^2)./(4*pi*c*RR); | ||
| + | |||
| + | surf(RR,TT,K) | ||
| + | shading interp | ||
| + | title("Solución fundamental dim 3 en función del radio") | ||
| + | xlabel("r") | ||
| + | ylabel("t") | ||
| + | |||
| + | |||
| + | }} | ||
| + | |||
| + | ==Aplicación de la solución fundamental== | ||
Ahora vamos a considerar el siguiente problema en dimensión 2: | Ahora vamos a considerar el siguiente problema en dimensión 2: | ||
| Línea 419: | Línea 532: | ||
| − | [[Archivo:Solfund0ej2.png|400px|thumb|center| Representación de la solución | + | [[Archivo:Solfund0ej2.png|400px|thumb|center| Representación de la solución para <math> t=0 </math>.]] |
| − | [[Archivo:Solfund0.5ej2.png|400px|thumb|center| Representación de la solución | + | [[Archivo:Solfund0.5ej2.png|400px|thumb|center| Representación de la solución para <math> t=0.5 </math>.]] |
| − | [[Archivo:Solfund1ej2.png|400px|thumb|center| Representación de la solución | + | [[Archivo:Solfund1ej2.png|400px|thumb|center| Representación de la solución para <math> t=1 </math>.]] |
| − | [[Archivo:Solfund2ej2.png|400px|thumb|center| Representación de la solución | + | [[Archivo:Solfund2ej2.png|400px|thumb|center| Representación de la solución para <math> t=2 </math>.]] |
| + | |||
| + | |||
A continuación, se deja un vídeo de la solución a lo largo del tiempo. | A continuación, se deja un vídeo de la solución a lo largo del tiempo. | ||
| − | |||
| − | Se puede observar | + | [[Archivo:Vídeo-Ej-2.2.gif|400px|thumb|center| Representación de la solución con <math> r \in [0,2] </math> y <math> t \in [0,4] </math>.]] |
| + | |||
| + | |||
| + | Se puede observar cómo, al contrario que en las primeras secciones, esta solución no es periódica. Esto se debe exclusivamente al dominio, ya que al ser infinito, no se produce un 'rebote' de la onda, y por ello no vuelve al punto de partida. | ||
La solución de esta ecuación es radial, no depende del ángulo, esto se debe a que las condiciones del problema se pueden poner en función del radio de la siguiente forma: | La solución de esta ecuación es radial, no depende del ángulo, esto se debe a que las condiciones del problema se pueden poner en función del radio de la siguiente forma: | ||
| Línea 448: | Línea 565: | ||
</math></center> | </math></center> | ||
| − | Por tanto la solución no depende del radio. | + | Por tanto la solución no depende del radio. Pero esto también se obtiene de que la solución es la convolución de dos funciones radiales, y por ello, se obtiene un resultado radial. |
| − | Además los programas se han hecho conociendo que la solución es radial, ya que reduce gran parte de la complejidad del código permitiendo que tarde menos. | + | Además los programas se han hecho conociendo que la solución es radial, ya que reduce gran parte de la complejidad del código, permitiendo que tarde menos. |
| + | |||
| + | ===Programas=== | ||
| + | |||
| + | Código para representar la solución: | ||
| + | |||
| + | {{matlab|codigo= | ||
| + | clear all | ||
| + | close all | ||
| + | |||
| + | %%% EJERCICIO 2.2 %%%% | ||
| + | |||
| + | c=1; % Velocidad de propagación. | ||
| + | |||
| + | % Definimos los intervalos de integración y sus discretizaciones | ||
| + | |||
| + | r0=0; rf=1; % Valor inicial y final de r | ||
| + | t0=0; tf=1; % Valor inicial y final de t | ||
| + | th0=0; thf=2*pi; % Valor inicial y final de theta | ||
| + | T=[0,0.5,1,2]; % Discretización del intervalo de theta | ||
| + | Th=linspace(th0,thf); % Discretización del intervalo de theta | ||
| + | R=linspace(r0,rf,50); % Discretización del intervalo del radio | ||
| + | |||
| + | [TH,RR]=meshgrid(Th,R); | ||
| + | |||
| + | s=size(TH); | ||
| + | |||
| + | U=zeros(s); % Inicializamos los valores de u | ||
| + | |||
| + | XX=RR.*cos(TH); | ||
| + | YY=RR.*sin(TH); | ||
| + | |||
| + | % Solución fundamental en polares preparada para la convolución | ||
| + | K=@(a,l,r,t)(sqrt(r.^2+l.^2-2*r.*l.*cos(a))<c*t)./(0.01*(c^2*t.^2-r.^2- l.^2 + 2.*r.*l.*cos(a)<10^-16)+2*pi*c*sqrt(abs(c^2*t.^2- r.^2- l.^2 + 2.*r.*l.*cos(a)))); | ||
| + | |||
| + | % Solución para cada valor de t | ||
| + | for j=1:length(T) | ||
| + | for i=1:length(R) | ||
| + | U(i,:)=integral2(@(a,l)K(a,l,R(i),T(j)), 0,2*pi, 0,1/2)*ones(1,length(Th)); % Evaluamos la integral | ||
| + | j,i | ||
| + | end | ||
| + | |||
| + | figure(j) | ||
| + | surf(XX,YY,U) | ||
| + | shading interp | ||
| + | title("Solución para dim 2") | ||
| + | subtitle("t="+num2str(T(j))) | ||
| + | |||
| + | end | ||
| + | }} | ||
| + | |||
| + | ==Conclusión== | ||
| + | La solución fundamental es una herramienta muy importante en la resolución de la ecuación de ondas, que además permite la interpretación del principio de Huygens. La relación entre ambos conceptos es que la propagación de un frente de onda puede ser vista como la superposición de las soluciones fundamentales emitidas desde cada punto del frente de onda anterior, conforme al principio de Huygens. | ||
=Referencias= | =Referencias= | ||
| Línea 455: | Línea 624: | ||
* [https://www.lifeder.com/velocidad-propagacion-onda/ Velocidad de propagación de una onda] | * [https://www.lifeder.com/velocidad-propagacion-onda/ Velocidad de propagación de una onda] | ||
* [https://mat.caminos.upm.es/wiki/Series_de_Fourier_(GRwM) Series de Fourier] | * [https://mat.caminos.upm.es/wiki/Series_de_Fourier_(GRwM) Series de Fourier] | ||
| − | + | * [https://www.fisicalab.com/apartado/principio-huygens Principio de Huygens] | |
[[Categoría:EDP]] | [[Categoría:EDP]] | ||
[[Categoría:EDP23/24]] | [[Categoría:EDP23/24]] | ||
Revisión actual del 10:50 26 may 2024
Contenido
1 Introducción
A lo largo de este trabajo vamos a representar diferentes soluciones de la ecuación de ondas en una dimensión, interpretaremos los resultados obtenidos. También vamos a estudiar la solución fundamental de dicha ecuación en 1, 2 y 3 dimensiones, gracias a la cual podremos interpretar el principio de Huygens.
Además, es importante destacar que todos los códigos y gráficas que se verán durante el trabajo han sido realizados con MatLab.
2 Conceptos previos
Antes de introducirnos de lleno en el trabajo, vamos a ver algunos conceptos que son necesarios para la comprensión de este:
Velocidad de propagación: La velocidad de propagación de una onda es la magnitud que mide la velocidad a la que se propaga la perturbación de la onda a lo largo de su desplazamiento. La velocidad a la que se propaga la onda depende tanto del tipo de onda como del medio por el que esta se propaga.
Principio de Huygens: El principio de Huygens es un concepto fundamental en la óptica y la física de ondas, formulado por el físico holandés Christiaan Huygens en 1678. Este principio establece que cada punto de un frente de onda actúa como una fuente secundaria de ondas esféricas, y que el frente de onda en un momento posterior se puede considerar como la envolvente de todas estas ondas secundarias.
Además, en este trabajo vamos a emplear series de Fourier.
3 Ecuación de ondas I
En esta sección, vamos a representar diferentes soluciones de la ecuación de ondas en una dimensión. Para ello, vamos a considerar una cuerda vibrante en el intervalo [math] [0,1][/math], con densidad [math] d[/math] y tensión [math] \tau_0 [/math] constante. De modo que la velocidad de propagación es [math] c=1 ~~m/s[/math] . Además, vamos a suponer que la cuerda está fija en los extremos y vamos a denotar por [math] u_0(x)[/math] y [math] u_1(x)[/math] su posición e impulso iniciales respectivamente.
3.1 Planteamiento del problema
Comenzamos planteando el sistema de EDP que modeliza el comportamiento de los desplazamientos transversales de la cuerda. Para ello, basta considerar la ecuación de ondas [math] u_{tt} – c^2u_{xx}=0[/math] y añadir las condiciones iniciales y de frontera descritas en el apartado anterior. De esta manera se obtiene que el sistema EDP a resolver es:
3.2 Resolución del sistema
Para resolver el sistema, vamos a aplicar el método de separación. Para ello, vamos a considerar que la solución es de la forma [math] u(x,t) = T(t) X(x)[/math] .
Haciendo las cuentas pertinentes y aplicando el principio de superposición obtenemos que la solución del sistema es:
siendo
[math] c_k= \frac{\int_0^1 sin (k \pi x) u_0 (x) dx}{\int_0 ^1 sin^2(k \pi x) dx}[/math]
[math] d_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]
Nótese que al considerar t como un múltiplo de dos, por la estructura de la solución, se obtiene de nuevo la condición inicial. Es decir, la solución es periódica con periodo dos. Este resultado lo analizaremos con más detalle en los siguientes apartados.
3.3 Particularización del problema
Ahora vamos a considerar como datos iniciales [math] u_0(x)= e^{-100 (x- \frac{1}{2})^2}[/math] y [math] u_1(x)= 0[/math] y vamos a representar la solución en el intervalo de tiempo [math] [0,2][/math]. Teniendo en cuenta la solución general obtenida en el caso anterior, se tiene que la solución de este problema es:
ya que en la expresión de [math] d_k[/math] aparece la función [math] u_1(x)[/math], cuyo valor es cero por hipótesis. Para obtener el valor de [math] c_k [/math] vamos a emplear el método del trapecio, que al ser un método numérico, lleva asociado un error. En este trabajo no estudiaremos este error, para un análisis detallado del mismo consulte el trabajo Ecuación de Laplace (GRwM).
Teniendo esto en cuenta, hemos representado la solución en el intervalo de tiempo mencionado:
Como podemos observar, la amplitud máxima de la onda es [math] 1 [/math] y se alcanza para valores enteros del tiempo. Además, para [math]x=0 [/math] y [math]x=1[/math] se tiene que [math]u(x,t)=0[/math], pues por hipótesis la cuerda está fija en los extremos.
Además, la solución es periódica en tiempo. Para poder apreciar esto, se presenta el siguiente vídeo.
Como podemos observar, al considerar la solución al variar el tiempo entre [math]0[/math] y [math]4[/math], partiendo de una posición, la onda vuelve a esa misma posición dos veces, es decir, su periodo es dos, como ya habíamos deducido en el apartado anterior.
Además, cuando la onda llega a cada uno de los extremos, como estos se mantienen fijos por hipótesis, la onda rebota simétricamente , sufriendo un cambio en el signo de la tensión.
3.3.1 Programa
Nota: En este programa y en los siguientes la solución se obtiene de forma numérica. Discretizaremos el espacio mediante una malla en la cual obtenemos un valor aproximado de la solución. En los siguientes programas los comentarios indicarán los valores que se pueden modificar y qué cambian.
Código para representar la solución del sistema de EDP's:
% EJERCICIO 1 ONDAS Ap3
clear all
close all
% Intervalos de definición
x0=0; xf=1; % Valor inicial y final de x
t0=0; tf=2; % Valor inicial y final de t
difx=10^-3; dift=10^-3; % División del intervalo de x y t respectivamente
X=x0:difx:xf; % Discretización del intervalo de x
T=t0:dift:tf; % Discretización del intervalo de x
[XX,TT]=meshgrid(X,T);
% Intervalo de integración
difint=10^-3;
Int=x0:difint:xf;
% Condiciones iniciales
u0=@(x)exp(-100*(x-1/2).^2);
U0=u0(Int);
U=zeros(size(XX));
N=1000; % Número de términos de la serie que sumamos
% Obtención de la solución
for k=1:N
S=sin(k*pi*XX);
ck=2*trapz(Int,U0.*S(1,:)); % Coeficientes de Fourier ck con la fórmula del trapecio
U=U+ck*cos(k*pi*TT).*S;
end
% Dibujamos
surf(XX,TT,U)
shading interp
title("Solución para t \in [0,2]")
xlabel("Espacio")
ylabel("Tiempo")
zlabel("u(x,t)")
3.4 Otra particularización del problema
Ahora vamos a considerar que la onda viaja en un solo sentido, es decir, [math] u(x,t)=f(x-t)[/math]. Para hacerlo, vamos a tomar como datos iniciales [math] u_0(x)= f(x)[/math] y [math] u_1(x)= - f’(x)[/math], con [math] f(x)= e^{-100 (x- \frac{1}{2})^2}[/math].
En este caso el sistema es:
Resolviendo al igual que en los casos anteriores, obtenemos una solución cuya representación es la siguiente:
Además, se tiene que el perfil de la solución avanza con una velocidad de propagación [math]c= 1~~ m/s [/math]. Para poder apreciar esto, veamos la gráfica anterior desde otra perspectiva:
Observamos que el tiempo que tarda en llegar la onda de un extremo a otro es una unidad de tiempo. Como además estamos considerando que la cuerda mide una unidad de espacio, la velocidad de propagación es uno.
Además, al igual que en los casos anteriores, vamos a representar un vídeo de cómo varía la solución en función del tiempo:
La solución es periódica de periodo dos. Se cumple que cuando la onda llega al extremo, rebota simétricamente, sufriendo un cambio en el signo de la tensión. Además, justo cuando la onda llega a la frontera, se puede observar cómo se anula la solución en todo el intervalo del espacio.
3.4.1 Programa
Código para representar la solución del sistema de EDP's anterior.
% EJERCICIO 1 ONDAS Ap4
clear all
close all
% Intervalos de definición
x0=0; xf=1; % Valor inicial y final de x
t0=0; tf=2; % Valor inicial y final de t
difx=10^-3; dift=10^-3; % División del intervalo de x y t respectivamente
X=x0:difx:xf; % Discretización del intervalo de x
T=t0:dift:tf; % Discretización del intervalo de x
[XX,TT]=meshgrid(X,T);
% Intervalo de integración
difint=10^-3;
Int=x0:difint:xf;
% Condiciones iniciales
u0=@(x)exp(-100*(x-1/2).^2);
U0=u0(Int);
u1=@(x)100*2*(x-1/2).*exp(-100*(x-1/2).^2);
U1=u1(Int);
U=zeros(size(XX));
N=1000; % Número de términos de la serie que sumamos
% Obtención de la solución
for k=1:N
S=sin(k*pi*XX);
s=S(1,:);
ck=2*trapz(Int,U0.*s); % Coeficientes de Fourier ck con la fórmula del trapecio
dk=2/(k*pi)*trapz(Int,U1.*s); % Coeficientes de Fourier dk con la fórmula del trapecio
U=U+(ck*cos(k*pi*TT)+dk*sin(k*pi*TT)).*S;
end
% Dibujamos
surf(XX,TT,U)
shading interp
title("Solución para t \in [0,2]")
xlabel("Espacio")
ylabel("Tiempo")
zlabel("u(x,t)")
3.5 Cambio a condición frontera de tipo Neumann
Por último, vamos a estudiar cómo varía la solución al considerar condiciones de frontera de tipo Neumann en lugar de tipo Dirichlet. Para ello, en este caso vamos a tomar [math]u_x(0,t) = u_x(1,t) =0[/math]. Entonces, tenemos el siguiente sistema:
En este caso, aplicando separación de variables obtenemos que la solución es:
siendo
[math] c_k= \frac{\int_0^1 cos (k \pi x) u_0 (x) dx}{\int_0 ^1 cos^2(k \pi x) dx}[/math]
[math] d_k= \frac{\int_0^1 cos (k \pi x) u_1 (x) dx}{\int_0 ^1 cos^2(k \pi x) dx}[/math]
Nótese que, al igual que para el problema con condiciones Dirichlet, en este caso la solución también es periódica de periodo dos.
La gráfica de la solución en este caso es:
Para poder apreciar mejor el comportamiento de la onda, hemos realizado el siguiente vídeo:
Como podemos observar, al igual que en los casos anteriores, la solución es de periodo [math] 2 [/math].
Además, en este caso, como los extremos no están fijos, cuando la onda llega a uno de los extremos, este varía. La razón de esto es que se debe cumplir que [math] u_x=0 [/math].
3.5.1 Programa
Representación de la solución del sistema con las condiciones de tipo Neumann:
% EJERCICIO 1 ONDAS Ap4
clear all
close all
% Intervalos de definición
x0=0; xf=1; % Valor inicial y final de x
t0=0; tf=2; % Valor inicial y final de t
difx=10^-3; dift=10^-3; % División del intervalo de x y t respectivamente
X=x0:difx:xf; % Discretización del intervalo de x
T=t0:dift:tf; % Discretización del intervalo de x
[XX,TT]=meshgrid(X,T);
% Intervalo de integración
difint=10^-3;
Int=x0:difint:xf;
% Condiciones iniciales
u0=@(x)exp(-100*(x-1/2).^2);
U0=u0(Int);
u1=@(x)100*2*(x-1/2).*exp(-100*(x-1/2).^2);
U1=u1(Int);
N=1000; % Número de términos de la serie que sumamos
% Obtención de la solución
U=trapz(U0,Int)*ones(size(XX));
for k=1:N
C=cos(k*pi*XX);
c=C(1,:);
ck=2*trapz(Int,U0.*c); % Coeficientes de Fourier ck con la fórmula del trapecio
dk=2/(k*pi)*trapz(Int,U1.*c); % Coeficientes de Fourier dk con la fórmula del trapecio
U=U+(ck*cos(k*pi*TT)+dk*sin(k*pi*TT)).*C;
end
% Dibujamos
surf(XX,TT,U)
shading interp
title("Solución para t \in [0,2]")
xlabel("Espacio")
ylabel("Tiempo")
zlabel("u(x,t)")
4 Ecuación de ondas II
En esta sección, vamos a estudiar la solución fundamental de la ecuación de ondas en dimensiones 1, 2 y 3 y una aplicación de esta en dimensión 2. Hay que tener en cuenta que la solución fundamental es la solución que se obtiene al dar un impulso inicial muy localizado en [math] x=0 [/math].
4.1 Solución fundamental
En primer lugar, consideramos el siguiente sistema:
donde [math] \delta(x) [/math] es la delta de Dirac. La expresión de la solución fundamental es:
[math]1)[/math] En dimensión [math]n=1[/math]:
Que se puede observar que tiene una expresión analítica radial, luego la solución fundamental para dimensión 1 es radial.
Gráficamente la solución fundamental de la ecuación de ondas en [math] dim =1 [/math] es:
[math]2)[/math] En dimensión [math]n=2[/math]:
donde [math] \chi_{B(0,ct)}(x) [/math] es la función característica de la bola de centro 0 y radio [math] ct [/math]. Como se puede observar, también tiene una expresión analítica radial y, por ello, es una solución radial.
Luego la representación de la solución fundamental en [math] dim=2[/math] es:
donde hemos considerado la siguiente regularización para evitar la singularidad:
[math]3)[/math] En dimensión [math]n=3[/math]:
Que ya de por sí tiene también una expresión analítica radial.
Por tanto, la solución fundamental en [math] dim=3[/math] se representa de la siguiente forma:
donde hemos sustituido la delta de Dirac por su aproximación
Con [math]k=1000[/math].
4.1.1 Programas
clear all
close all
Código de la solución fundamental para dimensión 1.
%%% EJERCICIO 2.1 dim1 %%%%
% Intervalos de definición
x0=-1; xf=1; % Valor inicial y final de x
t0=0; tf=1; % Valor inicial y final de t
difx=10^-2; dift=10^-3; % División del intervalo de x y t respectivamente
X=x0:difx:xf; % Discretización del intervalo de x
T=t0:dift:tf; % Discretización del intervalo de t
[XX,TT]=meshgrid(X,T);
% Velocidad de propagación
c=1;
% Cálculo de la solución fundamental
K=(XX<=c*TT).*(XX>=-c*TT)./(2*c);
surf(XX,TT,K)
shading interp
title("Solución fundamental dim 1")
xlabel("x")
ylabel("t")
Código para representar la solución fundamental para dimensión 2.
clear all
close all
%%% EJERCICIO 2.1 dim2 %%%%
% Intervalos de definición
r0=0; rf=1; % Valor inicial y final de r
t0=0; tf=1; % Valor inicial y final de t
divr=10^-2;divt=10^-3; % División del intervalo de r y t respectivamente
R=r0:divr:rf; % Discretización del intervalo de r
T=t0:divt:tf; % Discretización del intervalo de t
[RR,TT]=meshgrid(R,T);
% Velocidad de propagación
c=1;
% Cálculo de la solución fundamental
K=(RR<=c*TT)./(eps + 2*pi*c.*sqrt(c^2*TT.^2-RR.^2));
surf(RR,TT,K)
shading interp
title("Solución fundamental dim 2 en función del radio")
xlabel("r")
ylabel("t")
Código de la solución fundamental para dimensión 3.
clear all
close all
%%% EJERCICIO 2.1 dim3 %%%%
% Intervalos de definición
r0=0; rf=1; % Valor inicial y final de r
t0=0; tf=1; % Valor inicial y final de t
divr=10^-2;divt=10^-3; % División del intervalo de r y t respectivamente
R=r0:divr:rf; % Discretización del intervalo de r
T=t0:divt:tf; % Discretización del intervalo de t
[RR,TT]=meshgrid(R,T);
c=1; % Velocidad de propagación
k=1000; % Parámetro empleado en la aproximación de la delta de Dirac
% Cálculo de la solución fundamental
K=sqrt(k/pi).*exp(-k*(RR-c*TT).^2)./(4*pi*c*RR);
surf(RR,TT,K)
shading interp
title("Solución fundamental dim 3 en función del radio")
xlabel("r")
ylabel("t")
4.2 Aplicación de la solución fundamental
Ahora vamos a considerar el siguiente problema en dimensión 2:
cuya solución viene dada por la convolución
A continuación se muestra la solución obtenida para los tiempos [math] t=0,0.5,1,2 [/math]:
A continuación, se deja un vídeo de la solución a lo largo del tiempo.
Se puede observar cómo, al contrario que en las primeras secciones, esta solución no es periódica. Esto se debe exclusivamente al dominio, ya que al ser infinito, no se produce un 'rebote' de la onda, y por ello no vuelve al punto de partida.
La solución de esta ecuación es radial, no depende del ángulo, esto se debe a que las condiciones del problema se pueden poner en función del radio de la siguiente forma:
Por tanto la solución no depende del radio. Pero esto también se obtiene de que la solución es la convolución de dos funciones radiales, y por ello, se obtiene un resultado radial. Además los programas se han hecho conociendo que la solución es radial, ya que reduce gran parte de la complejidad del código, permitiendo que tarde menos.
4.2.1 Programas
Código para representar la solución:
clear all
close all
%%% EJERCICIO 2.2 %%%%
c=1; % Velocidad de propagación.
% Definimos los intervalos de integración y sus discretizaciones
r0=0; rf=1; % Valor inicial y final de r
t0=0; tf=1; % Valor inicial y final de t
th0=0; thf=2*pi; % Valor inicial y final de theta
T=[0,0.5,1,2]; % Discretización del intervalo de theta
Th=linspace(th0,thf); % Discretización del intervalo de theta
R=linspace(r0,rf,50); % Discretización del intervalo del radio
[TH,RR]=meshgrid(Th,R);
s=size(TH);
U=zeros(s); % Inicializamos los valores de u
XX=RR.*cos(TH);
YY=RR.*sin(TH);
% Solución fundamental en polares preparada para la convolución
K=@(a,l,r,t)(sqrt(r.^2+l.^2-2*r.*l.*cos(a))<c*t)./(0.01*(c^2*t.^2-r.^2- l.^2 + 2.*r.*l.*cos(a)<10^-16)+2*pi*c*sqrt(abs(c^2*t.^2- r.^2- l.^2 + 2.*r.*l.*cos(a))));
% Solución para cada valor de t
for j=1:length(T)
for i=1:length(R)
U(i,:)=integral2(@(a,l)K(a,l,R(i),T(j)), 0,2*pi, 0,1/2)*ones(1,length(Th)); % Evaluamos la integral
j,i
end
figure(j)
surf(XX,YY,U)
shading interp
title("Solución para dim 2")
subtitle("t="+num2str(T(j)))
end
4.3 Conclusión
La solución fundamental es una herramienta muy importante en la resolución de la ecuación de ondas, que además permite la interpretación del principio de Huygens. La relación entre ambos conceptos es que la propagación de un frente de onda puede ser vista como la superposición de las soluciones fundamentales emitidas desde cada punto del frente de onda anterior, conforme al principio de Huygens.