Ecuación de ondas. Otelo, Yan y Mika

1 Ecuación de ondas I

La idea principal de este artículo será dibujar diferentes soluciones de la ecuación de ondas en una dimensión. Para ello, nos basaremos en el enunciado descrito a continuación.

1.1 Planteamiento del problema

Se considera una cuerda vibrante que ocupa 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 es [math] c = \tau_0/d = 1 [/math] . Supondremos además 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.

Lo primero que haremos, será escribir el sistema de EDPS que modeliza el comportamiento de los desplazamiento transversales de la cuerda. Esto es:

[math] \left\{ \begin{array}{ll} u_{tt}-u_{xx}=0, \hspace{6mm} t\gt 0, x\in [0,1], \\ u(0,t)=u(1,t)=0, \hspace{6mm} t\gt0, \\ u(x,0)=u_0(x), \hspace{3mm} u_t(x,0)=u_1(x) \hspace{6mm} x\in [0,1]. \\ \end{array} \right. [/math]

Cuya solución encontrada gracias al método de separación de las variables, el principio de superposición y las condiciones iniciales es

[math] u(x,t)=\sum_{k=1}^\infty \left[\frac{1}{k\pi}u_{1,k}\sin(k\pi t)+u_{0,k}\cos(k\pi t)\right]\sin(k\pi x) [/math]

donde:

[math] u_{1,k} = 2 \int_0^1 u_1(x) \sin(k\pi x) dx [/math]

[math] u_{0,k} = 2 \int_0^1 u_0(x) \sin(k\pi x) dx[/math]

1.2 Dibujo de la solución dadas las condiciones iniciales

Considerando [math] u_0(x,t)= e^{-100(x-1/2)^2}, \hspace{3mm} u_1(x) = 0[/math], entonces la solución tomando los primeros [math]50 [/math] términos de la serie es

[math] u(x,t)=\sum_{k=1}^{50} 2 \left[\int_0^1 u_0(x) \sin(k\pi x) dx\right]cos(k\pi t) \sin(k\pi x) [/math]

Y el dibujo de la solución en el intervalo de tiempo [math] [0,2] [/math] queda así:

Solución ecuación de ondas

Podemos observar que solución para estas condiciones iniciales es periódica en el tiempo con periodo [math] T = 2 [/math] ya que en la gráfica se repite el patrón inicial.

Esto también se puede apreciar si repetimos el dibujo en el intervalo temporal [math] [0,4] [/math]:

Solución ecuación de ondas

1.3 Otros datos iniciales más particulares para el sistema

Ahora, consideraremos los datos iniciales correspondientes a una onda que viaja en un sólo sentido [math]\left(u(x,t) = f(x-t)\right) [/math]. Para ello, tomamos:

[math] u_0(x) = f(x), \hspace{3mm} u_1(x) = -f'(x) \hspace{2mm} con \hspace{2mm} f(x) = e^{-100(x-1/2)^2} [/math]

Es decir, la misma [math]u_0 [/math] que antes pero esta vez tomamos como [math]u_1 [/math] su derivada respecto el espacio cambiada de signo:

[math] u_0(x) = e^{-100(x-1/2)^2}, \hspace{3mm} u_1(x) = (200x-100)e^{-100(x-1/2)^2} \hspace{2mm} [/math] [math] u(x,t)=\sum_{k=1}^\infty \left[\frac{1}{k\pi}u_{1,k}\sin(k\pi t)+u_{0,k}cos(k\pi t)\right]\sin(k\pi x) [/math]

Onda que viaja en un solo sentido

Podemos observar que claramente en el dibujo se muestra una onda que viaja en un solo sentido con longitud de onda [math]\lambda = 2 [/math] metros. Por tanto la velocidad es [math]c = \frac{\lambda}{T} = \frac{2}{2} = 1 m/s [/math]

La onda parte de [math]x = 0.5 [/math]m y se dirige hacia [math]x=1 [/math]. En [math]t = 0.5 [/math]s llega a esa frontera y allí cambia de signo y de dirección, para dirigirse hacia el extremo [math]x=0[/math], que tarda en alcanzar [math]1 [/math]s. Por último, cuando llega a este extremo, pasa de nuevo a tomar valores positivos y toma la dirección inicial para volver a la posición de partida [math]x = 0.5 [/math]. De nuevo probando. que el periodo es de [math]2 [/math] segundos.

Todo esto se ve mejor observando la onda desde otro ángulo:

Onda que viaja en un solo sentido

Archivo:Ejercicio44.gif

En el video anterior se puede observar como tarda una unidad de tiempo en recorrer una unidad de espacio ya que el pico de la onda en t = 0 empieza en x = 0.5 y vuelve a x = 0.5 en t = 1. Sin embargo, una vez alcanzada la frontera de x = 1, la onda cambia de signo y de sentido. Este hecho vuelve a ocurrir cuando alcanza la próxima frontera x = 0, donde vuelve a tomar el signo positivo y el sentido inicial. Esto se repite cada vez que alcanza una frontera. Además, se puede observar que en t = 2, la onda vuelve a ser igual a la inicial, comprobando así que el periodo T es igual a 2.

1.4 Diferencia con la condición Neumann

Por último, veremos qué ocurre si repetimos lo anterior cambiando las condiciones Dirichlet por Neumann:

[math] u_x(0,t) = u_x(1,t) = 0 [/math]

Con ello, el sistema pasa a ser:

[math] \left\{ \begin{array}{ll} u_{tt}-u_{xx}=0, \hspace{6mm} t\gt 0, x\in [0,1], \\ u_x(0,t)=u_x(1,t)=0, \hspace{6mm} t\gt0, \\ u(x,0)=e^{-100(x-1/2)^2}, \hspace{3mm} u_t(x,0)=(200x-100)e^{-100(x-1/2)^2} \hspace{6mm} x\in [0,1]. \\ \end{array} \right. [/math]

Aplicamos de nuevo separación de las variables y el principio de superposición peor debido al cambio en las condiciones obtenemos esta solución:

[math] u(x,t)=\sum_{k=1}^\infty \left[\frac{1}{k\pi}u_{1,k}\sin(k\pi t)+u_{0,k}cos(k\pi t)\right]cos(k\pi x) [/math]

donde, ahora:

[math] u_{1,k} = 2 \int_0^1 u_1(x) \cos(k\pi x) dx [/math]

[math] u_{0,k} = 2 \int_0^1 u_0(x) \cos(k\pi x) dx[/math]

En este caso, la gráfica de la solución queda así:

Solución ecuación de ondas

En este caso, la velocidad de la onda es la misma y el desplazamiento ocurre en el mismo sentido y dirección. Lo único que cambia es que cuando la onda llega a los extremos del intervalo espacial, en vez de cambiar de signo, aumenta hasta el máximo y regresa al mismo valor con el cambio de dirección. Esto hace que en los extremos se cumpla la condición Neumann que hemos impuesto.

2 Ecuación de ondas II

En este apartado dibujaremos la solución fundamental de la ecuación de ondas en dimensiones 1, 2 y 3, lo que nos servirá para interpretar el principio de Huygens.

La solución fundamental es la solución que se obtiene al dar un impulso inicial muy localizado en [math] x=0 [/math]. Matemáticamente resuelve el sistema

[math] \left\{ \begin{array}{ll} u_{tt}-c^2\Delta u=0, \hspace{6mm} \hspace{32mm} x\in \mathbb{R}^n, t \gt0\\ u(x,0)=u_0(x), \hspace{3mm} u_t(x,0)=\delta (x) \hspace{6mm} 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_{r\to 0} \frac{1}{|B(0,r)|}\chi_{B(0,r)}(x) \sim \left\{ \begin{array}{ll} \infty \hspace{4mm} x=0,\\ 0 \hspace{4mm} 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.

Su expresión es:

1. 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)= \left\{ \begin{array}{ll} 0 \hspace{4mm} si \hspace{2mm}s\lt0,\\ 1 \hspace{4mm} si \hspace{2mm} s\geq 0,\\ \end{array} \right.[/math] es la función de Heaviside.

2. En dimensión [math] n=2 [/math]:

[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 0 y radio [math] ct [/math].

3. En dimensión [math] n=3 [/math]:

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

2.1 Representaciones

Antes de nada, nótese que la solución fundamental es radial para todas estas dimensiones. La primera lo es puesto que [math] x\in \mathbb{R}[/math], y las otras dos son radiales porque la variable espacial [math] x [/math] aparece siempre en módulo, [math] |x| [/math].

A continuación, dibujamos las dichas soluciones fundamentales:

Solución estacionaria

% Parámetros
c = 1;  % Puedes ajustar el valor de c según lo necesites

% Definir la función de Heaviside
H = @(s) double(s >= 0);

% Definir la función K1(x,t)
K1 = @(x, t) (1 / (2 * c)) * (H(x + c * t) - H(x - c * t));

% Rango de valores para x y t
x = linspace(-10, 10, 1000);  % Ajusta el rango y la resolución según lo necesites
t = linspace(0, 10, 1000);  % Ajusta el rango y la resolución según lo necesites

% Crear mallas para x y t
[X, T] = meshgrid(x, t);

% Calcular K1 para cada par de valores (x, t)
K1_values = arrayfun(@(x_val, t_val) K1(x_val, t_val), X, T);

% Dibujar la función
figure;
surf(X, T, K1_values, 'EdgeColor', 'none');

% Añadir más color
colormap(jet);  % Utilizar el colormap 'jet' para colores más vibrantes
colorbar;
shading interp;  % Interpolación de colores para suavizar la apariencia

% Añadir etiquetas y título
xlabel('x');
ylabel('t');
zlabel('K_1(x,t)');
title(['Gráfico de K_1(x,t) para c = ', num2str(c)]);

% Añadir rejilla y ajustar la vista
grid on;
view(3);  % Vista en 3D

% Ajustar propiedades de la figura
set(gca, 'FontSize', 12);  % Tamaño de fuente más grande para mejorar la legibilidad

Como podemos ver en la expresión analítica, esta función fundamental tiene una singularidad cuando [math] c^2t^2=|x|^2 [/math]. Por ello, para poder representarla añadiremos un sumando [math] \epsilon [/math] al denominador, esto es

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

Su representación gráfica en función de la variable radial es la dada por la siguiente imagen:

Solución estacionaria

% Parámetros
c = 1;  % Puedes ajustar el valor de c según lo necesites

% Definir la función característica de la bola de centro 0 y radio ct
chi_B = @(x, t) double(abs(x) <= c * t);

% Definir la función K2(x,t)
K2 = @(x, t) (1 ./ (0.01 + 2 * pi * c * sqrt(c^2 * t.^2 - x.^2))) .* chi_B(x, t);

% Rango de valores para x y t
x = linspace(0, 1, 1000);  % Ajusta el rango y la resolución según lo necesites
t = linspace(0, 1, 1000);  % Ajusta el rango y la resolución según lo necesites, evitando t=0 para evitar división por cero

% Crear mallas para x y t
[X, T] = meshgrid(x, t);

% Calcular K2 para cada par de valores (x, t)
K2_values = arrayfun(@(x_val, t_val) K2(x_val, t_val), X, T);

% Dibujar la función
figure;
surf(X, T, K2_values, 'EdgeColor', 'none');

% Añadir más color
colormap(jet);  % Utilizar el colormap 'jet' para colores más vibrantes
colorbar;
shading interp;  % Interpolación de colores para suavizar la apariencia

% Añadir etiquetas y título
xlabel('x');
ylabel('t');
zlabel('K_2(x,t)');
title(['Gráfico de K_2(x,t) para c = ', num2str(c)]);

% Añadir rejilla y ajustar la vista
grid on;
view(3);  % Vista en 3D

% Ajustar propiedades de la figura
set(gca, 'FontSize', 12);  % Tamaño de fuente más grande para mejorar la legibilidad

En dimensión [math] n=3 [/math], para poder hacer la representación tendremos que tomar la siguiente aproximación de la delta de Dirac:

[math] \delta(s) \sim \phi _k (s) = \sqrt{\frac{1000}{\pi}e^{-1000s^2}} [/math]

Así, su representación queda como sigue:

Solución estacionaria

% Parámetros
c = 1;  % Puedes ajustar el valor de c según lo necesites
a = 1000;  % Constante del exponente

% Definir la función K3(x,t)
K3 = @(x, t) (sqrt((a / pi) * exp(-a * (abs(x) - c * t).^2))) ./ (4 * pi * c * abs(x));

% Rango de valores para x y t
x = linspace(0, 1, 1000);  % Ajusta el rango y la resolución según lo necesites
t = linspace(0, 1, 1000);  % Ajusta el rango y la resolución según lo necesites, evitando t=0 para evitar división por cero

% Crear mallas para x y t
[X, T] = meshgrid(x, t);

% Calcular K3 para cada par de valores (x, t)
K3_values = arrayfun(@(x_val, t_val) K3(x_val, t_val), X, T);

% Dibujar la función
figure;
surf(X, T, K3_values, 'EdgeColor', 'none');

% Añadir más color
colormap(jet);  % Utilizar el colormap 'jet' para colores más vibrantes
colorbar;
shading interp;  % Interpolación de colores para suavizar la apariencia

% Añadir etiquetas y título
xlabel('x');
ylabel('t');
zlabel('K_3(x,t)');
title(['Gráfico de K_3(x,t) para c = ', num2str(c)]);

% Añadir rejilla y ajustar la vista
grid on;
view(3);  % Vista en 3D

% Ajustar propiedades de la figura
set(gca, 'FontSize', 12);  % Tamaño de fuente más grande para mejorar la legibilidad