Ecuación de ondas. Grupo Eau De Parfum(EDP)

Trabajo realizado por estudiantes
Título	Ecuación de Ondas. Grupo Eau De Parfum (EDP)
Asignatura	EDP
Curso	2023-24
Autores	Lestau Torres, Pablo López Rojo, Celia Muñoz Guijarro, Sofía
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

En este trabajo, representaremos e interpretaremos diferentes soluciones de la ecuación de ondas en una dimensión . Además, interpretaremos el principio de Huygens, estudiando la solución fundamental de la ecuación de ondas en 1, 2 y 3 dimensiones.

2 Contexto histórico

Jean le Rond D'Alembert

La ecuación de onda es una ecuación diferencial en derivadas parciales lineal de segundo orden, crucial para describir la propagación de diversas ondas, como las ondas sonoras, las ondas de luz y las ondas en el agua. Su relevancia se extiende a múltiples campos, incluyendo la acústica, el electromagnetismo, la mecánica cuántica y la dinámica de fluidos. Históricamente, Jean le Rond d'Alembert fue el primero en estudiar el problema de una cuerda vibrante, como las utilizadas en instrumentos musicales, en 1746.

3 Conocimientos previos

En esta sección se van a describir algunos conceptos a tener en cuenta durante todo el trabajo.

-Principio de Huygens: Afirma que cada punto en un frente de onda actúa como una fuente de ondas secundarias que se propagan hacia adelante con la misma velocidad que la propia onda. Estas ondas secundarias se combinan para formar un nuevo frente de onda, que es tangente a todas las ondas secundarias. Es decir, cada punto del frente de onda original puede considerarse como un punto de partida para nuevas ondas que se mueven en la misma dirección que la perturbación inicial. Estas ondas secundarias tienen la misma velocidad y frecuencia que la onda original.

Cabe destacar que vamos a utilizar Series de Furier

4 Ecuación de ondas I

Estamos considerando una cuerda vibrante que se extiende en el intervalo [0,1], con una densidad [math]d[/math] y una tensión constante [math]\tau_0[/math], lo que resulta en una velocidad de propagación [math]c=\frac{\tau_0}{d}=1[/math]. Además suponemos que la cuerda está fija en los extremos y llamamos [math]u_0(x)[/math] y [math]u_1(x)[/math] a su posición e impulso iniciales respectivamente.

El sistema de ecuaciones que modeliza el comportamiento de los desplazamientos transversales de la cuerda es:

[math] \left\{ \begin{aligned} &u_{tt}(x,t) = u_{xx}(x,t) & 0 \lt x \lt 1, t \gt 0, \\ &u(0, t) = u(1, t)=0, & t \gt 0, \\ &u(x, 0) =u_0(x) \\ &u_t(x, 0) = u_1(x), \end{aligned} \right. [/math]

La solución general de la ecuación de ondas empleando separación de variables se puede expresar en términos de la serie de Fourier de las condiciones iniciales [math]u_0(x)[/math], [math]u_1(x)[/math]. La solución tiene la forma

[math]u(x,t)=\sum_{n=1}^{\infty} [A_n \cos(n\pi t)+B_n \sin(n \pi t)]sin(n \pi x) [/math]

Donde los coeficientes [math]A_n[/math] y [math]B_n[/math] se determinan a partir de las condiciones iniciales [math]u_0(x)[/math], [math]u_1(x)[/math].Por lo que la solución completa de a ecuación de ondas para la cuerda vibrante con las condiciones iniciales dadas es:

[math]u(x,t)=\sum_{n=1}^{\infty}[ (2\int_{0}^{1} u_0(x) \sin(n \pi x) dx) \cos(n \pi t) + (\frac{2}{n\pi} \int_{0}^{1} u_0(x) \sin(n \pi x) dx)\sin(n\pi t)] \sin(n\pi x)[/math]

4.1 Ejemplo de la solución de la ecuación

Para resolver este problema, primero debemos calcular los coeficientes de Fourier [math] A_n [/math] y [math] B_n [/math] usando los datos iniciales [math] u_0(x)=e^{-100(x-1/2)^2} [/math] y [math] u_1(x)=0 [/math]. Luego, usaremos estos coeficientes para construir la solución aproximada [math] u(x,t) [/math] tomando los primeros 50 términos de la serie.

Los coeficientes de Fourier son:

[math] A_n = 2 \int_0^1 u_0(x) \sin(n \pi x) \, dx = 2 \int_0^1 e^{-100(x-1/2)^2} \sin(n \pi x) \, dx. [/math]

[math] B_n = \frac{2}{n \pi} \int_0^1 u_1(x) \sin(n \pi x) \, dx = 0 \quad (\text{porque } u_1(x) = 0). [/math]

Por lo tanto, la solución se simplifica a:

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

Procedemos a representar la solución obteniendo:

import numpy as np
import matplotlib.pyplot as plt
 
# Definimos las funciones iniciales u0(x) y u1(x)
def u0(x):
    return np.exp(-100 * (x - 0.5) ** 2)
 
def u1(x):
    return np.zeros_like(x)
 
# Parámetros
L = 1.0
c = 1.0  # Velocidad de propagación de las ondas
N = 100  # Número de puntos para el método del trapecio
x = np.linspace(0, L, N)
 
# Intervalo de tiempo y número de términos de Fourier
T = 2.0
n_terms = 50
 
# Resolución temporal
t_steps = 200
t = np.linspace(0, T, t_steps)
 
# Inicialización de u(x,t)
u_xt = np.zeros((N, t_steps))
 
# Calculamos la solución para cada tiempo t
for j, tj in enumerate(t):
    u_t = np.zeros_like(x)
    for n in range(1, n_terms + 1):
        # Calculamos los coeficientes A_n y B_n
        integrand_f = u0(x) * np.sin(n * np.pi * x / L)
        A_n = 2 * np.trapz(integrand_f, x)
 
        integrand_g = u1(x) * np.sin(n * np.pi * x / L)
        B_n = 2 / (n * np.pi * c) * np.trapz(integrand_g, x)
 
        # Suma de términos en la serie de Fourier
        u_t += (A_n * np.cos(n * np.pi * c * tj) + B_n * np.sin(n * np.pi * c * tj)) * np.sin(n * np.pi * x / L)
    u_xt[:, j] = u_t
 
# Graficamos la solución en diferentes tiempos
plt.figure(figsize=(10, 6))
for j in range(0, t_steps, t_steps // 10):
    plt.plot(x, u_xt[:, j], label=f't={t[j]:.2f}')
 
plt.xlabel('x')
plt.ylabel('u(x,t)')
plt.title('Desplazamientos transversales de la cuerda')
plt.legend()
plt.show()
 
# Comprobamos la periodicidad en el tiempo
plt.figure(figsize=(10, 6))
for i in range(0, N, N // 10):
    plt.plot(t, u_xt[i, :], label=f'x={x[i]:.2f}')
 
plt.xlabel('t')
plt.ylabel('u(x,t)')
plt.title('Periodicidad de la solución en el tiempo')
plt.legend()
plt.show()

La gráfica muestra la periodicidad en el tiempo de la solución [math]u(x,t) [/math]. La función claramente repite su patrón a medida que [math] t [/math] varía, lo que confirma que la solución es periódica en el tiempo con un período de 2. Este período corresponde al término más lento [math] n=1 [/math], ya que [math]\cos(\pi t)[/math] tiene un período de 2.

Representándolo en tres dimensiones obtenemos

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, PillowWriter
 
 
# Graficamos la solución en 3D
X, T = np.meshgrid(x, t)
 
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection='3d')
 
# Transponemos u_xt para que las dimensiones coincidan con las de X y T
ax.plot_surface(X, T, u_xt.T, cmap='viridis')
 
ax.set_xlabel('x')
ax.set_ylabel('t')
ax.set_zlabel('u(x,t)')
ax.set_title('Desplazamientos transversales de la cuerda')
 
plt.show()

Mostramos a continuación una animación que muestra los desplazamientos trasversales de la cuerda en función del instante de tiempo. Donde además confirmamos el periodo 2 de la solución obtenida.

# Configuración de la animación
fig, ax = plt.subplots()
line, = ax.plot(x, u_xt[:, 0])
ax.set_xlim(0, L)
ax.set_ylim(-1.10, 1.10)
ax.set_xlabel('x')
ax.set_ylabel('u(x,t)')
ax.set_title('Desplazamientos transversales de la cuerda')
 
# Agregar un texto que muestre el instante de tiempo
time_text = ax.text(0.02, 0.95, '', transform=ax.transAxes)
 
def update(frame):
    line.set_ydata(u_xt[:, frame])
    time_text.set_text(f't = {t[frame]:.2f}s')
    return line, time_text
 
ani = FuncAnimation(fig, update, frames=t_steps, blit=True)
 
# Guardar la animación como GIF
writer = PillowWriter(fps=20)
ani.save("wave_equation.gif", writer=writer)
 
plt.show()

4.2 Onda en un solo sentido

Veamos ahora lo que ocurre cuando estudiamos una onda que viaja en un solo sentido, es decir, [math] u(x, t) = f (x − t)[/math]. Que particularizaremos tomando 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-1/2)^2}[/math]. Por tanto el sistema con el que trabajamos es

[math] \left\{ \begin{aligned} &u_{tt}(x,t) = u_{xx}(x,t) & 0 \lt x \lt 1, t \gt 0, \\ &u(0, t) = u(1, t)=0, & t \gt 0, \\ &u(x, 0) =e^{-100(x-1/2)^2} \\ &u_t(x, 0) = 200(x-1/2)e^{-100(x-1/2)^2}, \end{aligned} \right. [/math]

Su solución viene dada mediante el método de separación de variables en términos de los coeficientes de Fourier:

[math] A_n = 2 \int_0^1 u_0(x) \sin(n \pi x) \, dx = 2 \int_0^1 e^{-100(x-1/2)^2} \sin(n \pi x) \, dx. [/math]

[math] B_n = \frac{2}{n \pi} \int_0^1 u_1(x) \sin(n \pi x) \, dx =\frac{2}{n \pi} \int_0^1 (200(x-1/2)e^{-100(x-1/2)^2})\sin(n \pi x) \, dx . [/math]

Por lo que la solución completa de la ecuación de ondas para la cuerda vibrante con las condiciones iniciales dadas es:

[math]u(x,t)=\sum_{n=1}^{\infty}[ ( 2 \int_0^1 e^{-100(x-1/2)^2} \sin(n \pi x) \, dx) \cos(n \pi t) + (\frac{2}{n \pi} \int_0^1 (200(x-1/2)e^{-100(x-1/2)^2})\sin(n \pi x) \, dx)\sin(n\pi t)] \sin(n\pi x)[/math]

La representación gráfica de esta solución es:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, PillowWriter
 
# Definimos la función f(x)
def f(x):
    return np.exp(-100 * (x - 0.5) ** 2)
 
# Definimos la derivada de f(x)
def primef(x):
    return -200 * (x - 0.5) * np.exp(-100 * (x - 0.5) ** 2)
 
# Parámetros
L = 1.0
c = 1.0  # Velocidad de propagación de las ondas
N = 300  # Número de puntos para el método del trapecio
x = np.linspace(0, L, N)
 
# Intervalo de tiempo y número de términos de Fourier
T = 2.0
n_terms = 50
 
# Resolución temporal
t_steps = 600
t = np.linspace(0, T, t_steps)
 
# Inicialización de u(x,t)
u_xt = np.zeros((N, t_steps))
 
# Calculamos la solución para cada tiempo t
for j, tj in enumerate(t):
    u_t = np.zeros_like(x)
    for n in range(1, n_terms + 1):
        # Calculamos los coeficientes A_n y B_n
        integrand_f = f(x) * np.sin(n * np.pi * x / L)
        A_n = 2 * np.trapz(integrand_f, x)
 
        integrand_g = primef(x) * np.sin(n * np.pi * x / L)
        B_n = 2 / (n * np.pi * c) * np.trapz(integrand_g, x)
 
        # Suma de términos en la serie de Fourier
        u_t += (A_n * np.cos(n * np.pi * c * tj) + B_n * np.sin(n * np.pi * c * tj)) * np.sin(n * np.pi * x / L)
   
    u_xt[:, j] = u_t
 
# Configuración de la animación
fig, ax = plt.subplots()
line, = ax.plot(x, u_xt[:, 0])
ax.set_xlim(0, L)
ax.set_ylim(-1.10, 1.10)
ax.set_xlabel('x')
ax.set_ylabel('u(x,t)')
ax.set_title('Desplazamientos transversales de la cuerda')
 
# Agregar un texto que muestre el instante de tiempo
time_text = ax.text(0.02, 0.95, '', transform=ax.transAxes)
 
def update(frame):
    line.set_ydata(u_xt[:, frame])
    time_text.set_text(f't = {t[frame]:.2f}s')
    return line, time_text
 
ani = FuncAnimation(fig, update, frames=t_steps, blit=True)
 
# Guardar la animación como GIF
writer = PillowWriter(fps=20)
ani.save("wave_equation_dirichlet.gif", writer=writer)
 
plt.show()

Podemos comprobar que el perfil de la solución "viaja" con velocidad [math] c = 1 [/math] observando la animación proporcionada. La propagación de la onda se evidencia por el desplazamiento de la forma de la onda a lo largo del eje [math] x [/math] a medida que avanza el tiempo. Si observamos cualquier punto fijo en la onda, veremos cómo se mueve a lo largo del eje [math] x [/math] con una velocidad constante de una unidad por unidad de tiempo, lo que confirma que la onda se propaga con velocidad [math] c = 1 [/math].

En cuanto a lo que ocurre cuando la onda llega a la frontera, en el caso de condiciones de frontera de tipo Dirichlet, la onda se refleja en la frontera. Cuando la onda alcanza los extremos del dominio (en x=0 y x=1), la función de onda u(x,t) queda fija en cero. Además cuando la onda que llega a la frontera se refleja con amplitud y fase opuestas.

4.3 Cambio en condición frontera

Para comprender la diferencia entre las condiciones de frontera Dirichlet y Neumann en la ecuación de ondas, repetiremos el análisis del apartado anterior, esta vez aplicando condiciones de Neumann. En lugar de fijar los valores de la función en los extremos (condiciones de Dirichlet), estableceremos que las derivadas parciales respecto a [math]x[/math] en los extremos sean cero. Es decir, consideramos las condiciones [math] u_x(0,t)=u_x(1,t)=0[/math].

Por tanto nuestro nuevo problema es

[math] \left\{ \begin{aligned} &u_{tt}(x,t) = u_{xx}(x,t) & 0 \lt x \lt 1, t \gt 0, \\ &u_x(0, t) = u_x(1, t)=0, & t \gt 0, \\ &u(x, 0) =e^{-100(x-1/2)^2} \\ &u_t(x, 0) = 200(x-1/2)e^{-100(x-1/2)^2} \end{aligned} \right. [/math]

Para resolver esta ecuación utilizando el método de Fourier, se expande la función [math] u(x,t) [/math] en una serie de Fourier:

[math] u(x,t) = \sum_{n=1}^{\infty} \left( A_n \cos(\frac{n\pi ct}{L} ) + B_n \sin(\frac{n\pi ct}{L} ) \right) \cos \left( \frac{n\pi x}{L} \right) [/math]

donde [math] L [/math] es la longitud de la cuerda y los coeficientes [math] A_n [/math] y [math] B_n [/math] se calculan mediante integrales utilizando las funciones iniciales [math] f(x) [/math] y [math] f'(x) [/math].

Estas funciones se utilizan para calcular los coeficientes \( A_n \) y \( B_n \), y luego se utilizan en la serie de Fourier para obtener la solución \( u(x,t) \).

[math] A_n = 2 \int_{0}^{L} f(x) \cos \left( \frac{n\pi x}{L} \right) \, dx [/math]

[math] B_n = \frac{2}{n\pi c} \int_{0}^{L} f'(x) \cos \left( \frac{n\pi x}{L} \right) \, dx [/math]

La gran diferencia con la solución por serie de Fourier del apartado anterior, radica en la resolución del problema de autovalores para la variable x.

Comenzamos buscando soluciones de la forma [math]U(x, t) = w(t) v(x)[/math] con [math]v(0) = v(L) = 0[/math]. Insertando [math]U[/math] en la ecuación de onda, encontramos:

[math] 0 = U_{tt} - c^2 U_{xx} = w''(t) v(x) - c^2 w(t) v''(x) [/math]

O, separando las variables,

[math] \frac{1}{c^2} \frac{w''(t)}{w(t)} = \frac{v''(x)}{v(x)} [/math]

Es una identidad entre dos funciones, una dependiendo solo de [math]t[/math] y la otra solo de [math]x[/math]. Por lo tanto, ambos lados de la igualdad deben ser iguales a la misma constante, digamos [math]\lambda[/math]. Así, llegamos a la ecuación

[math] w''(t) - \lambda c^2 w(t) = 0 [/math]

y al problema de los autovalores

[math] v''(x) - \lambda v(x) = 0 [/math]

[math] v'(0) = v'(L) = 0 [/math]

Para la solución del problema de los autovalores. Hay tres posibilidades.

1. Si [math]\lambda = 0[/math], entonces [math]v(x) = A + Bx[/math] esto implica [math]A = B = 0[/math].
2. Si [math]\lambda \gt 0[/math], entonces [math]v(x) = Ae^{-\sqrt\lambda x} + Be^{\sqrt\lambda x}[/math] y nuevamente implica [math]A = B = 0[/math].
3. Si [math]\lambda \lt 0[/math], entonces [math]v(x) = A \sin(\sqrt\lambda x) + B \cos(\sqrt\lambda x)[/math].

De las condiciones obtenemos:

[math] v'(0) = \sqrt\lambda A = 0 [/math]

por tanto [math]A =0 [/math].

[math] v'(1) = A\sqrt{\lambda}\cos(\sqrt{\lambda}) - B\sqrt{\lambda}\sin(\sqrt{\lambda}) = 0 [/math]

de donde

[math] B \text{ arbitrario}, A = 0, \lambda = \frac{m^2\pi^2}{L^2}, n = 1, 2, \ldots [/math]

Por lo tanto, en el caso c) solo encontramos soluciones no triviales, de la forma.

[math] v_m(x) = B_n \cos(\lambda_n x), \quad \lambda_n = \frac{m^2\pi^2}{L^2} [/math]

Resolviendo el otro problema de forma análoga, llegamos a la solución general mencionada anteriormente.

[math] u(x,t) = \sum_{n=1}^{\infty} \left( A_n \cos(\frac{n\pi ct}{L} ) + B_n \sin(\frac{n\pi ct}{L} ) \right) \cos \left( \frac{n\pi x}{L} \right) [/math]

Donde los coeficientes [math] A_n, B_n [/math] ya están definidos.

Particularizando a nuestros datos:

[math] u(x,t) = \sum_{n=1}^{\infty} \left( A_n \cos(n\pi t) + B_n \sin(n\pi t) \right) \cos(n\pi x) \quad A_n = 2 \int_{0}^{1} f(x) \cos(n\pi x) \, dx \quad B_n = \frac{2}{n\pi} \int_{0}^{1} f'(x) \cos(n\pi x) \, dx [/math]

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, PillowWriter
 
# Definimos la función f(x)
def f(x):
    return np.exp(-100 * (x - 0.5) ** 2)
 
# Definimos la derivada de f(x)
def primef(x):
    return -200 * (x - 0.5) * np.exp(-100 * (x - 0.5) ** 2)
 
# Parámetros
L = 1.0
c = 1.0  # Velocidad de propagación de las ondas
N = 300  # Número de puntos para el método del trapecio
x = np.linspace(0, L, N)
 
# Intervalo de tiempo y número de términos de Fourier
T = 2.0
n_terms = 50
 
# Resolución temporal
t_steps = 300
t = np.linspace(0, T, t_steps)
 
# Inicialización de u(x,t)
u_xt = np.zeros((N, t_steps))
 
# Calculamos la solución para cada tiempo t
for j, tj in enumerate(t):
    u_t = np.zeros_like(x)
    for n in range(1, n_terms + 1):
        # Calculamos los coeficientes A_n y B_n
        integrand_f = f(x) * np.cos(n * np.pi * x / L)
        A_n = 2 * np.trapz(integrand_f, x)
 
        integrand_g = primef(x) * np.cos(n * np.pi * x / L)
        B_n = 2 / (n * np.pi * c) * np.trapz(integrand_g, x)
 
        # Suma de términos en la serie de Fourier
        u_t += (A_n * np.cos(n * np.pi * c * tj) + B_n * np.sin(n * np.pi * c * tj)) * np.cos(n * np.pi * x / L)
   
    u_xt[:, j] = u_t
 
# Configuración de la animación
fig, ax = plt.subplots()
line, = ax.plot(x, u_xt[:, 0])
ax.set_xlim(0, L)
ax.set_ylim(-0.25, 2)
ax.set_xlabel('x')
ax.set_ylabel('u(x,t)')
ax.set_title('Desplazamientos transversales de la cuerda con condiciones de Neumann')
 
# Agregar un texto que muestre el instante de tiempo
time_text = ax.text(0.02, 0.95, '', transform=ax.transAxes)
 
def update(frame):
    line.set_ydata(u_xt[:, frame])
    time_text.set_text(f't = {t[frame]:.2f}s')
    return line, time_text
 
ani = FuncAnimation(fig, update, frames=t_steps, blit=True)
 
# Guardar la animación como GIF
writer = PillowWriter(fps=20)
ani.save("wave_equation_neumann.gif", writer=writer)
 
plt.show()

En el caso de condiciones de frontera de tipo Neumann, la onda no queda fija en los extremos. En lugar de eso, cuando la onda llega a uno de los extremos rebota y no cambia de dirección, es decir no se invierte. Esto se debe a que [math] u_x=0 [/math] en los extremos.

5 Ecuación de ondas II

En esta sección, abordaremos la resolución de la ecuación de ondas en dimensiones 1, 2 y 3, una tarea que nos permitirá explorar la solución fundamental y, a su vez, comprender el principio de Huygens.

5.1 Soluciones fundamentales para diferentes dimensiones

Comenzamos estudiando la solución fundamental, que surge como respuesta a un impulso inicial localizado en \( x = 0 \). Formalmente, se define mediante límites y funciones características. El objetivo es observar y graficar las soluciones fundamentales en su forma radial para cada dimensión, lo que nos ayudará a comprender mejor la naturaleza de las ondas en diferentes entornos dimensionales.

5.1.1 Solución fundamental en dimensión 1

La solución fundamental en dimensión 1, \( K_1(x, t) \), se define como:

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

donde \( H(s) \) es la función de Heaviside, dada por:

[math] H(s) = \begin{cases} 0, & \text{si } s \lt 0 \\ 1, & \text{si } s \geq 0 \end{cases} [/math]

Dicha función, \( K_1(x, t) \), se compone de dos funciones de Heaviside desplazadas horizontalmente a lo largo del eje \( x \) en \( ct \) y \( -ct \) y luego divididas por \( 2c \). Esto significa que en el intervalo \( (-\infty, -ct) \), \( K_1(x, t) \) es 0, luego aumenta a 1 en el intervalo \( (-ct, ct) \), y finalmente vuelve a 0 en el intervalo \( (ct, \infty) \). La constante \( c \) controla la velocidad de propagación de la onda, que hemos exigido \(1\). Graficamos la función para tener una primera idea intuitiva de la solución fundamental:

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
 def heaviside(x):
    return np.where(x >= 0, 1, 0)
 def fundamental_solution_1D(x, t, c):
    return (1 / (2 * c)) * (heaviside(x + c * t) - heaviside(x - c * t))
 # Definimos los parámetros
c = 1
t_values = [0.5, 1.0, 1.5]  # Valores de tiempo que queremos graficar
x_range = np.linspace(-5, 5, 100)
 # Graficamos las soluciones fundamentales en 1D, 2D y 3D para los diferentes valores de tiempo
fig, axes = plt.subplots(1, len(t_values), figsize=(15, 5))
 for i, t in enumerate(t_values):
    # Solución en 1D
    axes[i].plot(x_range, fundamental_solution_1D(x_range, t, c))
    axes[i].set_title(f'Solución Fundamental en 1D (t={t})')
    axes[i].set_xlabel('x')
    axes[i].set_ylabel('K(x, t)') 
plt.tight_layout()
plt.show()

Visualmente, la gráfica de \( K_1(x, t) \) parece una función escalón, donde hay un cambio abrupto en el valor de la función en \( x = \pm ct \). Esto refleja la propagación de una onda con velocidad \( c \) a lo largo del eje \( x \). Además, debido a que su valor depende únicamente de la distancia radial desde el origen, es decir, de la magnitud de \( x \), podemos decir que se trata de una función radial y, por tanto, graficarla de esta forma. En la expresión de \( K_1(x, t) \), la variable \( x \) aparece únicamente como \( x \pm ct \), lo que permite su reescritura como \( K_1(r, t) \), donde \( r = |x| \) es la distancia radial.

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
def heaviside(x):
    return np.where(x >= 0, 1, 0)
def fundamental_solution_1D_radial(r, t, c):
    return (1 / (2 * c)) * (np.heaviside(r + c * t, 0) - np.heaviside(r - c * t, 0))
# Definimos los parámetros
c = 1
t_values = np.arange(0, 5, 0.01)  # Valores de tiempo en incrementos de 0.01
r_values = np.linspace(0, 5, 500)
# Creamos una función para actualizar el gráfico en cada cuadro de la animación
def update(frame):
    plt.cla()  # Limpiamos el gráfico actual
    plt.plot(r_values, fundamental_solution_1D_radial(r_values, t_values[frame], c),color = 'blue')
    plt.title(f'Solución Fundamental en 1D - t={t_values[frame]}')
    plt.xlabel('r')
    plt.ylabel('K(r, t)')
    plt.grid(True)
# Creamos la animación
fig = plt.figure(figsize=(8, 6))
ani = animation.FuncAnimation(fig, update, frames=len(t_values), interval=50)
# Guardamos la animación como un GIF
ani.save('fundamental_solution_1D.gif', writer='pillow', fps=30)
plt.show()

Al graficar \( K_1(r, t) \) como función de \( r \) para diferentes valores de \( t \), se observa que la forma de la curva es la misma independientemente de la dirección específica en la que nos movamos desde el origen, lo que confirma su carácter radial.

5.1.2 Solución fundamental en dimensión 2

La solución fundamental en dimensión 2, \( K_2(x, t) \), se define como:

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

donde \(\chi_B(0,ct)(x)\) es la función característica de la bola de centro \(\ 0 \) y radio \(\ c \cdot t \). Esta función describe cómo una perturbación inicial, localizada en el centro de una bola de radio \( ct \), se propaga a lo largo del plano. La función característica \( \chi_{B(0, ct)}(x) \) asegura que la perturbación solo se propague dentro de esta bola. La presencia del término \( \frac{1}{c^2 t^2 - |x|^2} \) indica cómo la amplitud de la onda varía con la distancia desde el centro de la bola y el tiempo. Mostraremos gráficamente esta función para tener una primera idea intuitiva. Se mostrará un video que varía en función del tiempo:

import numpy as np
import matplotlib.pyplot as plt
import imageio
from io import BytesIO
def chi_B(r, ct):
    return np.where(r <= ct, 1, 0)
def fundamental_solution_2D(x, y, t, c, epsilon=0.01):
    r = np.sqrt(x ** 2 + y ** 2)
    discriminant = c ** 2 * t ** 2 - r ** 2
    solution = 1 / (epsilon + 2 * np.pi * c * np.sqrt(discriminant)) * chi_B(r, c * t)
    return np.where(discriminant > 0, np.minimum(solution, 1), 0)
# Definimos los parámetros
c = 1
epsilon = 0.01
t_values = np.linspace(0.1, 5, 50)  # Valores de tiempo a animar
x = np.linspace(-5, 5, 1000)
y = np.linspace(-5, 5, 1000)
# Generamos los frames y los guardamos
frames = []
for t in t_values:
    fig, ax = plt.subplots(figsize=(8, 8))
    X, Y = np.meshgrid(x, y)
    Z = fundamental_solution_2D(X, Y, t, c, epsilon)
    ax.contourf(X, Y, Z, cmap='viridis') #Greys?
    ax.set_title(f'Solución Fundamental en 2D para t={t}')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.grid(True)
    # Convertir la figura a una imagen numérica
    buf = BytesIO()
    fig.savefig(buf, format='png')
    buf.seek(0)
    frames.append(imageio.imread(buf))
    plt.close(fig)
# Guardamos los frames como un gif
imageio.mimsave('fundamental_solution_cartesiana_2D_2D.gif', frames, fps=10)

Podemos notar que la circunferencia amarilla no es uniforme, sino que muestra irregularidades. Esto se debe a que la solución fundamental presenta discontinuidades cuando \(\left| x \right| \rightarrow t\). Para visualizar esto mejor, en lugar de mostrar solo el contorno de la solución, vamos a representarla en tres dimensiones:

import numpy as np
import matplotlib.pyplot as plt
import imageio
from mpl_toolkits.mplot3d import Axes3D
def chi_B(r, ct):
    return np.where(r <= ct, 1, 0) 
def fundamental_solution_2D(x, y, t_values, c, epsilon=0.01):
    solutions = []
    for t in t_values:
        r = np.sqrt(x ** 2 + y ** 2)
        discriminant = c ** 2 * t ** 2 - r ** 2
        solution = 1 / (epsilon + 2 * np.pi * c * np.sqrt(discriminant)) * chi_B(r, c * t)
        solutions.append(np.where(discriminant > 0, np.minimum(solution, 1), 0))
    return solutions
# Definimos los parámetros
c = 1
epsilon = 0.01
t_values = np.linspace(0.1, 5, 50)  # Valores de tiempo a animar
x = np.linspace(-5, 5, 1000)
y = np.linspace(-5, 5, 1000)
X, Y = np.meshgrid(x, y)
# Generamos las soluciones para cada valor de tiempo
solutions = fundamental_solution_2D(X, Y, t_values, c, epsilon)
# Creamos los frames
frames = []
for t, solution in zip(t_values, solutions):
    fig = plt.figure(figsize=(10, 8))
    ax = fig.add_subplot(111, projection='3d')
    ax.plot_surface(X, Y, solution, cmap='viridis')
    ax.set_title(f'Solución Fundamental en 2D para t={t:.2f}')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('K(x, y, t)')
    ax.set_zlim(0, 1)
    plt.tight_layout()
    fig.canvas.draw()  # Dibujamos la figura para actualizarla
    frame = np.array(fig.canvas.renderer._renderer)
    plt.close(fig)
    frames.append(frame)
# Guardamos el GIF
imageio.mimsave('fundamental_solution_cartesiana_2D_3D.gif', frames, fps=10)

Es importante señalar que la discontinuidad no aparece de la misma manera en todos los ángulos, lo que podría sugerir que no es radial. Sin embargo, esta aparente discrepancia se debe a que estamos utilizando un mallado rectangular en lugar de uno polar en el código. A pesar de esto, podemos concluir que la función es claramente radial, lo que nos permite representar una sección transversal de la misma. En otras palabras, podemos visualizarla en función de la variable radial:

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
 
def chi_B(r, ct):
    return np.where(r <= ct, 1, 0)
 
def fundamental_solution_2D_radial(r, t, c, epsilon=0.01):
    discriminant = c ** 2 * t ** 2 - r ** 2
    solution = 1 / (epsilon + 2 * np.pi * c * np.sqrt(discriminant))* chi_B(r, c * t)
    return np.where(discriminant > 0, np.minimum(solution, 1), 0)
 
# Definimos los parámetros
c = 1
t_values = np.arange(0, 5, 0.01)  # Valores de tiempo en incrementos de 0.01
r_values = np.linspace(0, 5, 500)
 
# Creamos una función para actualizar el gráfico en cada cuadro de la animación
def update(frame):
    plt.cla()  # Limpiamos el gráfico actual
    plt.plot(r_values, fundamental_solution_2D_radial(r_values, t_values[frame], c), color='blue')
    plt.ylim(0, 1.01)  # Ajustamos el eje y en el rango de 0 a 1
    plt.title(f'Solución Fundamental en 2D - t={t_values[frame]}')
    plt.xlabel('r')
    plt.ylabel('K(r, t)')
    plt.grid(True)
 
# Creamos la animación
fig = plt.figure(figsize=(8, 6))
ani = animation.FuncAnimation(fig, update, frames=len(t_values), interval=50)
 
# Guardamos la animación como un GIF
ani.save('fundamental_solution_2D.gif', writer='pillow', fps=30)
 
plt.show()

Observamos la presencia de una asíntota en \( K_2(x, t) \). Esta asíntota refleja cómo la amplitud de la onda se acerca a cero cuando el denominador se aproxima a cero. Sin embargo, a medida que el tiempo avanza y el radio \( |x| \) aumenta, la amplitud de la onda disminuye. Esta disminución en la amplitud se interpreta como una disminución en la intensidad o "altura" de la onda (del agua, por ejemplo) a medida que se propaga a través del espacio en dos dimensiones.

5.1.3 Solución fundamental en dimensión 3

La solución fundamental en dimensión 3, \( K_3(x, t) \), se define como:

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

donde \( \delta(|x| - ct) \) representa la función delta de Dirac. Esta función modela cómo una perturbación inicial, localizada en un radio \( ct \) del origen, se propaga en el espacio tridimensional. La presencia de la función delta de Dirac garantiza que la perturbación solo se propague en un radio \( ct \) del origen. Comenzamos, como en los dos apartados anteriores, graficando la función. En este caso, no se puede representar en un mismo video cómo varía en función de \( x,y,z \) y \(t\) de manera clara. De modo que se ha optado por representar en un primer lugar un video de la relación entre \( x,y,z \) fijando \( t=1 \):

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation 
def dirac_approximation(s, k):
    return np.sqrt(k / (np.pi)) * np.exp(-k * s ** 2) 
def fundamental_solution_3D(x, y, z, t, c, k=1000):
    r = np.sqrt(x ** 2 + y ** 2 + z ** 2)
    return dirac_approximation(r - c * t, k) / (4 * np.pi * c * r)
# Parámetros
t = 1.0
c = 1.0
k = 1000
x = np.linspace(-5, 5, 500)
y = np.linspace(-5, 5, 500)
X, Y = np.meshgrid(x, y)
# Valores de z
z_values = np.linspace(-1, 1, 50)
# Crear figura y ejes
fig, ax = plt.subplots(figsize=(8, 8)) 
def update(frame):
    z = z_values[frame]
    ax.clear()
    Z = fundamental_solution_3D(X, Y, z, t, c, k)
    cset = ax.contourf(X, Y, Z, cmap='viridis')
    ax.set_title(f'z = {z:.2f}')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
ani = animation.FuncAnimation(fig, update, frames=len(z_values), interval=100) 
# Guardar la animación como un GIF
ani.save('fundamental_solution_3D_corte_z.gif', writer='pillow', fps=10)
plt.show()

Observamos que se trata de una sección de esfera hueca, debido a la presencia de la función delta de Dirac en el numerador. ¿Qué pasa si fijamos \(z\) y variamos \(t\)?

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation 
def dirac_approximation(s, k):
    return np.sqrt(k / (np.pi)) * np.exp(-k * s ** 2)
def fundamental_solution_3D(x, y, z, t, c, k=1000):
    r = np.sqrt(x ** 2 + y ** 2 + z ** 2)
    return dirac_approximation(r - c * t, k) / (4 * np.pi * c * r)
# Parámetros
z = 0
c = 1.0
k = 1000
x = np.linspace(-5, 5, 500)
y = np.linspace(-5, 5, 500)
X, Y = np.meshgrid(x, y) 
# Valores de z
t_values = np.linspace(0, 5, 50) 
# Crear figura y ejes
fig, ax = plt.subplots(figsize=(8, 8)) 
def update(frame):
    t = t_values[frame]
    ax.clear()
    Z = fundamental_solution_3D(X, Y, z, t, c, k)
    cset = ax.contourf(X, Y, Z, cmap='viridis')
    ax.set_title(f't = {t:.2f}')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
ani = animation.FuncAnimation(fig, update, frames=len(z_values), interval=100) 
# Guardar la animación como un GIF
ani.save('fundamental_solution_3D_dif_t.gif', writer='pillow', fps=10)
plt.show()

Obtenemos también una circunferencia que se va expandiendo con el tiempo. Uniendo ambas ideas podemos ver que en dimensión 3, la solución fundamental es una esfera hueca que con el tiempo va aumentando de radio y disminuyendo la intensidad del primer impulso.

Claramente es una función que depende únicamente del radio, así pues, de manera análoga, la representamos en su variable radial.

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
 
def dirac_approximation(s, k):
    return (k / np.pi) ** 0.5 * np.exp(-k * s ** 2)
 
def fundamental_solution_3D_radial(r, t, c, k=1000):
    return dirac_approximation(r - c * t, k) / (4 * np.pi * c * r)
 
# Definimos los parámetros
c = 1
t_values = np.arange(0, 5, 0.01)  # Valores de tiempo en incrementos de 0.01
r_values = np.linspace(0, 5, 500)
 
# Creamos una función para actualizar el gráfico en cada cuadro de la animación
def update(frame):
    plt.cla()  # Limpiamos el gráfico actual
    plt.plot(r_values, fundamental_solution_3D_radial(r_values, t_values[frame], c), color= 'blue')
    plt.ylim(0, 5)  # Ajustamos el eje y en el rango de 0 a 1
    plt.title(f'Solución Fundamental en 3D - t={t_values[frame]}')
    plt.xlabel('r')
    plt.ylabel('K(r, t)')
    plt.grid(True)
 
# Creamos la animación
fig = plt.figure(figsize=(8, 6))
ani = animation.FuncAnimation(fig, update, frames=len(t_values), interval=50)
 
# Guardamos la animación como un GIF
ani.save('fundamental_solution_3D.gif', writer='pillow', fps=30)
 
plt.show()

Este video puede intuir que un impulso solo es percibido por el individuo en el momento en que la esfera lo alcanza, y con una intensidad que va disminuyendo. Esto es fundamental en la comprensión de la propagación de ondas, conocido como el principio de Huygens.

5.2 Solución fundamental del problema en dimensión 2

En esta sección nos centraremos en resolver el siguiente problema de la ecuación de ondas en dos dimensiones:

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

La solución se obtendrá mediante una convolución entre la solución fundamental, explicada en la sección anterior, y \(u_t(x,0)\). Es decir, la solución se calcula de la siguiente manera:

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

Dado que hemos demostrado previamente que \(K_2(x, t)\) es radial, es evidente que \(K_2(x-y, t)\) también lo es. Además, consideremos \(h(x)=\chi_{B(0,\frac{1}{2})}(x)\), la función característica de una bola centrada en \(0\). Esta función también es radial. Por lo tanto, podemos simplificar nuestros cálculos mediante un cambio a coordenadas polares.

Para resolver la integral de la función \( u(x,t) \) en coordenadas polares, primero necesitamos reformular la integral en estos términos, donde

[math] \left\{ \begin{array}{ll} x = (r_x \cos \theta_x, r_x \sin \theta_x) \\ y = (r_y \cos \theta_y, r_y \sin \theta_y) \end{array} \right. [/math]

La integral se transforma usando el jacobiano de la transformación polar \( r \, dr \, d\theta \):

[math] u(r_x, \theta_x, t) = \int_0^{2\pi} \int_0^{\infty} K_2(r_x, r_y, \theta_x, \theta_y, t) h(r_y \cos \theta_y, r_y \sin \theta_y) r_y \, dr \, d\theta [/math]

Podemos simplificar nuestros cálculos y calcular la integral para una sección transversal, es decir, fijando un ángulo. Esto se debe a que para cualquier ángulo diferente, la representación obtenida seguirá teniendo la misma forma. Podemos exigir, por ejemplo \( \theta_x = 0\) y representar \(x=(r_x ,0)\). De modo, que mi nueva integral a calcular será la siguiente, tomando la velocidad de propagación \(c=1\):

[math] u(r_x, t) = \int_0^{2\pi} \int_0^{\infty} \frac{r_y}{2\pi \sqrt{t^2 - r^2}} \cdot \chi_{B(0,t)}(r) \cdot \chi_{B(0,\frac{1}{2})}(r_y)\, dr_y \, d\theta_y, \, \text{donde} \, r = \sqrt{r_x^2 + r_y^2 - 2r_x r_y \cos(\theta_y)} [/math]

Podemos simplificar nuevamente la integral, pues observamos que el integrando está ponderado por una función característica dependiente de la variable de integración, la cual está definida en todo \(\mathbb{R}^+\). La ventaja de tener una función con soporte compacto implica que se puede reducir el extremo de integración a un intervalo compacto, en este caso \(r_y = 0\) hasta \(\frac{1}{2}\). Por lo tanto, obtenemos la integral equivalente:

[math] u(r_x, t) = \int_0^{2\pi} \int_0^{\frac{1}{2}} \frac{r_y}{2\pi \sqrt{t^2 - r_x^2 - r_y^2 + 2r_x r_y \cos(\theta_y)}} \cdot \chi_{B(0,t)} (\sqrt{r_x^2 + r_y^2 - 2r_x r_y \cos(\theta_y)} ) \, dr_y \, d\theta_y [/math]

A continuación mostramos la solución obtenida para diferentes tiempos \(t\)

import numpy as np
from scipy.integrate import dblquad
import matplotlib.pyplot as plt
 
# Definimos la función del integrando
def K(r_y, theta_y, r_x, t):
    norma = r_x**2 + r_y**2 - 2*r_x*r_y*np.cos(theta_y)
    numerator = r_y
    denominator = np.where(t**2 - norma <= 10e-5, 10e-5, 2*np.pi*np.sqrt(abs(t**2 - norma)))
    return np.where(np.sqrt(norma) <= t, numerator / denominator, 0)
 
# Definimos la función que calcula la integral doble
def integral(rx, t):
    # Límites de integración para ry y theta_y
    lower_limit_ry = 0
    upper_limit_ry = 0.5
    lower_limit_theta_y = 0
    upper_limit_theta_y = 2 * np.pi
    result, error = dblquad(K, lower_limit_theta_y, upper_limit_theta_y,
                            lambda theta_y: lower_limit_ry, lambda theta_y: upper_limit_ry,
                            args=(rx, t))
    return result
 
# Definimos los valores de rx y t
rx_values = np.linspace(0, 1, 100)
T = [0, 0.5, 1, 2]
 
for t in range(len(T)):
    # Calculamos los valores de la integral para cada rx
    integral_values = [integral(rx, T[t]) for rx in rx_values]
 
    # Representamos los resultados
    plt.plot(rx_values, integral_values, label='Integral en función de $r_x$')
    plt.xlabel('$r_x$')
    plt.ylabel('Integral')
    plt.ylim(-0.1,3)
    plt.title(f'Evolución de la sección transversal en 2D para $t$={T[t]}')
    plt.legend()
    plt.grid(True)
    plt.show()