Ecuación de Ondas (Grupo 2 1/2)

Trabajo realizado por estudiantes
Título	Ecuación de Ondas.
Asignatura	EDP
Curso	2023-24
Autores	Alfredo de Lorenzo, Hugo Sanz, Manuel Fdez.
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

A lo largo de este trabajo se va a trabajar con diferentes soluciones de la ecuación de ondas, interpretando sus comportamientos. También se va a estudiar la solución fundamental en 1, 2 y 3 dimensiones aplicando un impulso inicial en 𝑥=0, gracias a lo cual se podrá interpretar el principio de Huygens. Donde este es:

Principio de Huygens: Todo punto de un frente de onda inicial puede considerarse como una fuente de ondas esféricas secundarias que se extienden en todas las direcciones con la misma velocidad, frecuencia y longitud de onda que el frente de onda del que proceden.

2 Ecuación de ondas I

En esta sección, se va a trabajar con la ecuación de ondas en una dimensión. Para ello, se va 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, se va a suponer que la cuerda está fijada en los extremos y se denota como [math] u_0(x)[/math] y [math] u_1(x)[/math] su posición e impulso iniciales respectivamente.

2.1 Planteamiento del problema

En primer lugar se plantea el sistema de EDP que modeliza el comportamiento de los desplazamientos transversales de la cuerda. Para ello, es suficiente considerar la ecuación de ondas [math] u_{tt} – c^2u_{xx}=0[/math] y se establecen las condiciones iniciales mencionadas anteriormente. De esta manera se obtiene que el sistema EDP a resolver es:

[math]\left \{ \begin{array}{ll} u_{tt} – u_{xx} =0 \\ u(0,t)=u(1,t)=0 \\ u(x,0)=u_0(x) \\ u_t(x,0)=u_1(x) \end{array} \right. [/math]

2.2 Resolución del sistema

La resolución de la ecuación se va realizar mediante el método de separación de variables. Para ello, se considera que la solución es de la forma [math] u(x,t) = T(t) X(x)[/math] .

Realizando el desarrollo correspondiente y aplicando el principio de superposición se obtiene que la solución del sistema es de la siguiente forma:

[math] u(x,t)= \sum _{k=1}^{\infty} \left(c_k cos (k \pi t) + d_k sin (k \pi x)\right) sin(k \pi t) [/math]

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{\int_0^1 sin (k \pi x) u_1 (x) dx}{\int_0 ^1 sin^2(k \pi x) dx}[/math]

Finalmente, si agrupamos todo, la solución de la ecuación de ondas es la siguiente:

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

Esta expresión proporciona la solución completa de la ecuación de onda en función de las condiciones iniciales \( u_0(x) \) y \( u_1(x) \).

2.3 Particularización del problema

En esta parte del documento se van a analizar una serie de ejemplos para comprender de manera precisa lo explicado en la sección anterior.

2.3.1 Caso 1. Solución periódica

El primero de ellos se trata de la representación gráfica de una solución periódica en tiempo. Esta viene dada suponiendo que los datos iniciales son [math]u_0(x)=e^{-100(x-\frac{1}{2})^2}[/math] y [math]u_1(x)=0[/math]. La solución de este problema se va a representar en el intervalo \(t ∈ [0, 2] \).

Representación Solución Caso 1

Como se puede observar, la amplitud máxima de la onda es 1 y se alcanza para valores enteros del tiempo. Además, para x=0 y x=1 se tiene que u(x,t)=0, pues por hipótesis la cuerda está fijada en los extremos.

Representación Solución Caso 1

Además también se ve como la función muestra un comportamiento que se repite regularmente con el tiempo. A estas funciones se les denomina periódicas en tiempo. Para ver esto más fácilmente se introduce el siguiente gif.

Representación Solución en función del tiempo Caso 1

Finalmente se muestra el código con el que se ha realizado dicha representación.

clear
close all
clc

% Se definen las condiciones iniciales
u_0 = @(x) exp(-100*(x-1/2).^2);
u_1 = @(x) 0;

% Definición de los intervalos de tiempo y espacio
t = linspace(0,2,1000);
x = linspace(0,1,1000);

% Número de términos de la serie a dibuja
n = 50;

% definición de matrices para almacenar coeficientes de Fourier
c = zeros(n,1);
d = zeros(n,1);

% Cálculo de coeficientes de Fourier asociados a u_0
for k = 1:n
    f_0 = @(x) u_0(x).*sin(k*pi*x);
    c(k) = 2*trapz(x, f_0(x));
end

% Cálculo de coeficientes de Fourier asociados a u_1
for k = 1:n
    f_1 = @(x) u_1(x).*sin(k*pi*x);
    d(k) =(2/(k*pi))*trapz(x, f_1(x));
end

% Sustituimos en la solución calculada por separación de variables
u = @(xx,tt) 0;
for i = 1:n
    u = @(xx,tt) u(xx,tt) + (c(i)*cos(i*pi*tt) + d(i)*sin(i*pi*tt)) .* sin(i*pi*xx);
end

% Graficamos la solución
[X, T] = meshgrid(x, t);
surf(X, T, u(X, T), 'EdgeColor', 'flat')
title('Solución')
xlabel('x')
ylabel('t')

2.3.2 Caso 2. Onda viajera

Se considera ahora que la onda viaja en un solo sentido, esto es que:

[math] u(x,t)=f(x-t) [/math]

Donde para las condiciones iniciales en este caso se utilizará que [math] u_0 (x)=f(x) [/math] y [math] u_1 (x)=-f'(x) [/math] y particularizando se utilizará la función:

[math] f(x)=e^{-100(x-\frac{1}{2})^{2}} [/math]

Derivando se obtiene:

[math] f'(x)=-200(x-\frac{1}{2})e^{-100(x-\frac{1}{2})^{2}} [/math]

Por lo que con esto ya se tiene todo lo necesario, de tal forma que la ecuación en derivadas parciales es:

[math]\left \{ \begin{array}{ll} u_{tt} – u_{xx} =0 \\ u(0,t)=u(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]

Ahora se procede a resolver el ejercicio, que se hará de manera análoga a los apartados anteriores.

clear
close all
clc

% Se definen las condiciones iniciales
u_0 = @(x) exp(-100*(x-1/2).^2);
u_1 = @(x) 200*(x-1/2).*exp(-100*(x-1/2).^2);

% Definición de los intervalos de tiempo y espacio
t = linspace(0,2,1000);
x = linspace(0,1,1000);

% Número de términos de la serie a dibuja
n = 50;

% definición de matrices para almacenar coeficientes de Fourier
c = zeros(n,1);
d = zeros(n,1);

% Cálculo de coeficientes de Fourier asociados a u_0
for k = 1:n
    f_0 = @(x) u_0(x).*sin(k*pi*x);
    c(k) = 2*trapz(x, f_0(x));
end

% Cálculo de coeficientes de Fourier asociados a u_1
for k = 1:n
    f_1 = @(x) u_1(x).*sin(k*pi*x);
    d(k) =(2/(k*pi))*trapz(x, f_1(x));
end

% Sustituimos en la solución calculada por separación de variables
u = @(xx,tt) 0;
for i = 1:n
    u = @(xx,tt) u(xx,tt) + (c(i)*cos(i*pi*tt) + d(i)*sin(i*pi*tt)).*sin(i*pi*xx);
end

% Graficamos la solución
[X, T] = meshgrid(x, t);
surf(X, T, u(X, T), 'EdgeColor', 'flat')
title('Solución')
xlabel('x')
ylabel('t')

Representación Solución en función del tiempo

Viendo la solución, se puede ver cómo la onda va con una velocidad constante de 1m/s, y tiene un periodo de 2 segundos. Cuando la onda llega a la frontera, cambia de dirección y es negativa. A continuación se ve una gráfica de la solución de perfil, donde se puede apreciar mejor la velocidad, y como cambia el comportamiento de la onda.

Representación Solución de Perfil para apreciar la velocidad

2.3.3 Caso 3. Condición frontera Neumann

Por último se va a trabajar cambiando las condiciones a unas de tipo Neumann, esto es que [math] u_x(0,t)=u_x(1,t)=0 [/math] como condiciones frontera. Con esto la ecuación en derivadas parciales es así:

[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]

Por ello, a lo hora de obtener coeficientes de Fourier para la solución, se calculan de forma distinta, haciendo separación de variables se obtiene:

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

Donde los coeficientes [math] c_k[/math] y [math] d_k [/math] son:

[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]

Ahora se adapta el código cambiando los coeficientes, y el problema se resuelve de manera análoga

clear
close all
clc

% Se definen las condiciones iniciales
u_0 = @(x) exp(-100*(x-1/2).^2);
u_1 = @(x) 200*(x-1/2).*exp(-100*(x-1/2).^2);

% Definición de los intervalos de tiempo y espacio
t = linspace(0,2,1000);
x = linspace(0,1,1000);

% Número de términos de la serie a dibuja
n = 50;

% definición de matrices para almacenar coeficientes de Fourier
c = zeros(n,1);
d = zeros(n,1);

% Cálculo de coeficientes de Fourier asociados a u_0
for k = 1:n
    f_0 = @(x) u_0(x).*cos(k*pi*x);
    c(k) = 2*trapz(x, f_0(x));
end

% Cálculo de coeficientes de Fourier asociados a u_1
for k = 1:n
    f_1 = @(x) u_1(x).*cos(k*pi*x);
    d(k) =2/(k*pi)*trapz(x, f_1(x));
end

% Sustituimos en la solución calculada por separación de variables
u = @(xx,tt) 0;
for i = 1:n
    u = @(xx,tt) u(xx,tt) + (c(i)*cos(i*pi*tt) + d(i)*sin(i*pi*tt)).*cos(i*pi*xx);
end

% Graficamos la solución
[X, T] = meshgrid(x, t);
surf(X, T, u(X, T), 'EdgeColor', 'flat')
title('Solución')
xlabel('x')
ylabel('t')

Representación Solución en función del tiempo

Aquí se puede ver cómo la onda tiene una velocidad constante de 1 casi todo el tiempo, pero cuando la onda llega a la frontera la onda se dispara y disminuye su velocidad, y posteriormente pasa a ir en el sentido contrario.

3 Ecuación de ondas II

En este ejercicio se dibujará la solución fundamental de la ecuación de ondas en dimensiones 1, 2 y 3. De esta forma, se interpretará el principio de Huygens.

3.1 Solución fundamental

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

El problema a resolver resulta el siguiente

[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.

La solución fundamental es:

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

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

Donde H(s) es la función de Heaviside, que vale 0 para [math]s\lt0[/math] y 1 para [math]s\geq0[/math]

Como se puede observar en la expresión, esta es radial.

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 en la bola de centro 0 y radio [math] ct [/math]. Esta expresión vuelve a ser de radial ya que depende del módulo de x.

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

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

La expresión vuelve a ser radial.

3.1.1 Representación

En dimensión 1:

Solución en dimensión 1

% Limpieza del entorno de trabajo
clear;
clc;

% Configuración inicial
x_vals = -1:0.01:1;     % Vector de valores de x
t_vals = 0:0.001:1;     % Vector de valores de t

% Crear la rejilla de puntos
[X_grid, T_grid] = meshgrid(x_vals, t_vals);

% Parámetro de velocidad
prop_speed = 1;

% Cálculo de la función fundamental
K_fun = (X_grid <= prop_speed * T_grid) & (X_grid >= -prop_speed * T_grid);
K_fun = K_fun / (2 * prop_speed);

% Visualización
figure;
surf(X_grid, T_grid, K_fun);
shading interp;
title('Solución fundamental en 1D');
xlabel('x');
ylabel('t');

En dimensión 2:

En esta dimensión hay una singularidad en el origen por lo que donde se ha considerado la regularización:

[math] K^{\varepsilon}_2 (x,t) = \frac{1}{\varepsilon + 2 \pi c \sqrt{c^2t^2 -|x|^2}} \chi _{B(0,1/2)}(x), ~~~~~\varepsilon = 0.01 [/math].

Solución en dimensión 2

% Limpiar el entorno de trabajo
clear;
clc;

% Configuración de los intervalos
r_vals = 0:0.01:1;    % Vector de valores de r
t_vals = 0:0.001:1;    % Vector de valores de t

% Crear la rejilla de puntos
[R_grid, T_grid] = meshgrid(r_vals, t_vals);

% Parámetro de velocidad
prop_speed = 1;

% Cálculo de la función fundamental
K_func = (R_grid <= prop_speed * T_grid) ./ (eps + 2 * pi * prop_speed .* sqrt(prop_speed^2 * T_grid.^2 - R_grid.^2));

% Visualización
figure;
surf(R_grid, T_grid, K_func);
shading interp;
title('Solución fundamental en 2D');
xlabel('r');
ylabel('t');

En dimensión 3:

Para poder representar la solución, se ha sustituido la delta de Dirac por su aproximación

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

con [math]k=1000[/math].

Solución en dimensión 3

% Limpiar el entorno de trabajo
clear;
clc;

% Definición de intervalos
% Configuración de los intervalos
r_vals = 0:0.01:1;    % Vector de valores de r
t_vals = 0:0.001:1;    % Vector de valores de t

% Crear la rejilla de puntos
[R_grid, T_grid] = meshgrid(r_values, t_values);

% Parámetros del modelo
prop_speed = 1;              % Velocidad de propagación
kappa = 1000;                % Parámetro para la aproximación de la delta de Dirac

% Cálculo de la función fundamental
K_function = sqrt(kappa / pi) .* exp(-kappa * (R_grid - prop_speed * T_grid).^2) ./ (4 * pi * prop_speed * R_grid);

% Visualización
figure;
surf(R_grid, T_grid, K_function);
shading interp;
title('Solución fundamental en 3D');
xlabel('r');
ylabel('t');

3.2 Problema en dimensión 2

Se considera el problema en dimensión 2:

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

La solución se obtene por la convolución siguiente

[math] u(x,t)=\int_{\mathbb{R}^2}K_2 (x-y, t)h(y) dy [/math].

Dado que, como se ha visto anteriormente, [math]K_2[/math] y [math]h[/math] son funciones radiales, su convolución también es radial.

Se representa la solución para [math]t=0, 0'5,1,2[/math]

Solución en diferentes tiempos

% Limpiar el entorno de trabajo y cerrar todas las figuras
clear;
clc;
close all;

% Parámetro de velocidad de propagación
prop_speed = 1;

% Definición de intervalos y sus discretizaciones
time_steps = [0, 0.5, 1, 2];     % Valores discretos de tiempo para evaluar
theta_vals = linspace(0, 2 * pi); % Discretización de theta
r_vals = linspace(0, 1, 50); % Discretización del radio

[Theta_grid, R_grid] = meshgrid(theta_vals, r_vals); % Crear la rejilla de puntos

% Inicializar la matriz de solución U
U = zeros(size(Theta_grid));

% Conversión a coordenadas cartesianas para la visualización
X_cartesian = R_grid .* cos(Theta_grid);
Y_cartesian = R_grid .* sin(Theta_grid);

% Función para la solución fundamental en coordenadas polares
K_fun = @(a, l, r, t) (sqrt(r.^2 + l.^2 - 2*r.*l.*cos(a)) < prop_speed * t) ./ ...
    (0.01*(prop_speed^2 * t.^2 - r.^2 - l.^2 + 2*r.*l.*cos(a) < 10^-16) + ...
    2 * pi * prop_speed * sqrt(abs(prop_speed^2 * t.^2 - r.^2 - l.^2 + 2*r.*l.*cos(a))));

% Cálculo de la solución para cada valor de tiempo en time_steps
for t_idx = 1:length(time_steps)
    for r_idx = 1:length(r_vals)
        U(r_idx, :) = integral2(@(a, l) K_fun(a, l, r_vals(r_idx), time_steps(t_idx)), ...
                               0, 2 * pi, 0, 0.5) * ones(1, length(theta_vals)); % Evaluar la integral
    end
    
    % Visualización de la solución
    figure;
    surf(X_cartesian, Y_cartesian, U);
    shading interp;
    title(['Solución en 2D, t = ', num2str(time_steps(t_idx))]);
    xlabel('x');
    ylabel('y');
end

3.3 Conclusión

A medida que se ha ido observando las gráficas, se ha podido apreciar la propagación de las ondas. Esta, es una superficie continua que se propaga hacia afuera y cuya forma es radial. La gráfica muestra claramente cómo la onda se expande y cómo su forma se mantiene coherente con el paso del tiempo, ilustrando el principio de Huygens. Por lo que de esta forma, se ha podido obtener una mejor interpretación de la ecuación de ondas a lo largo de la sección.