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.

COMPROBAR QUE EL PERFIL DE LA SOLUCIÓN VIAJA CON VELOCIDAD C=1 Y OBSERVAR QUÉ OCURRE EN LA FRONTERA

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

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} t\gt 0, x\in \mathbb{R}^n, \\ 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 por que 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