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}sen(k\pi t)+u_{0,k}cos(k\pi t)\right]sen(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) sen(k\pi x) [/math]

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

DIBUJO CON EJES BUENOS

La solución para estas condiciones iniciales es periódica en el tiempo:

COMPROBACIÓN

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}sen(k\pi t)+u_{0,k}cos(k\pi t)\right]sen(k\pi x) [/math]

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]

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