Ecuación de ondas

Trabajo realizado por estudiantes
Título	Ecuación de Ondas. Grupo ALA
Asignatura	EDP
Curso	2023-24
Autores	Lucía Amores, Aitana Guill y Andrea Navarro
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

Este trabajo analiza la ecuación de ondas en una cuerda fija en los extremos del intervalo [0,1], con velocidad de propagación 𝑐=1. Incluiremos el sistema que describe el movimiento, así como su solución por separación de variables. Además, consideraremos ejemplos específicos de simulaciones para observar fenómenos como la periodicidad y ondas viajeras. Finalmente, compararemos resultados cambiando las condiciones de frontera a Neumann.

A parte de esto, estudiaremos la solución fundamental de la ecuación de ondas en dimensiones 1, 2 y 3, aplicando un impulso inicial en 𝑥=0 y dibujando estas soluciones en la variable radial para comprender el principio de Huygens.

2 Contexto Histórico

La ecuación de ondas, fundamental en física y matemáticas, se formuló en el siglo XVIII por Jean le Rond d'Alembert y Leonhard Euler, quienes investigaron la vibración de cuerdas. En el siglo XIX, Joseph Fourier revolucionó su análisis al desarrollar la teoría de series de Fourier, facilitando la solución de ecuaciones diferenciales parciales mediante la separación de variables.

El principio de Huygens, propuesto por Christiaan Huygens en el siglo XVII, establece que cada punto de un frente de onda actúa como una fuente de ondas secundarias. Este principio es esencial para comprender la propagación de ondas y las soluciones de la ecuación en diversas dimensiones.

En el siglo XX y XXI, la teoría de la ecuación de ondas se ha ampliado a campos como la acústica, óptica y telecomunicaciones. Los métodos numéricos modernos permiten simular la propagación de ondas en medios complejos, mejorando nuestra comprensión y aplicación en problemas prácticos.

Hoy, la ecuación de ondas es clave para modelar fenómenos desde vibraciones de cuerdas hasta ondas sísmicas y electromagnéticas, siendo esencial en física teórica, ingeniería y tecnologías modernas. Este trabajo se centra en su estudio y soluciones, aplicando estos principios a situaciones prácticas y avanzadas.

3 Modelización de los desplazamientos transversales

Para modelar el comportamiento de los desplazamientos transversales de la cuerda vibrante, utilizamos la ecuación de ondas en una dimensión. Dado que la cuerda está fija en los extremos y tiene una densidad [math] d [/math] y tensión constante [math] \tau_0 [/math] de manera que la velocidad de propagación es [math] c = \tau_0/d = 1 [/math], la ecuación de ondas se simplifica.

La ecuación de ondas en una dimensión para los desplazamientos transversales [math] u(x,t) [/math] de la cuerda es:

[math]\frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2} [/math]

Dado que [math] c = 1 [/math], la ecuación se reduce a:

[math] \frac{\partial^2 u}{\partial t^2} = \frac{\partial^2 u}{\partial x^2} [/math]

Esta ecuación diferencial parcial (EDP) debe acompañarse de las condiciones de contorno y las condiciones iniciales para estar completamente especificada.

3.1 Condiciones de contorno

Dado que la cuerda está fija en los extremos, las condiciones de contorno son:

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

3.2 Condiciones iniciales

Las condiciones iniciales especifican la posición inicial de la cuerda [math] u_0(x) [/math] y su velocidad inicial o impulso [math] u_1(x) [/math]:

[math]u(x, 0) = u_0(x) [/math] [math] \frac{\partial u}{\partial t}(x, 0) = u_1(x) [/math]

3.3 Sistema de EDP

Juntando la ecuación de ondas con las condiciones de contorno e iniciales, el sistema completo que modela el comportamiento de los desplazamientos transversales de la cuerda es:

[math] \begin{cases} \frac{\partial^2 u}{\partial t^2} = \frac{\partial^2 u}{\partial x^2}, & 0 \lt x \lt 1, \ t \gt 0 \\ u(0, t) = 0, & t \geq 0 \\ u(1, t) = 0, & t \geq 0 \\ u(x, 0) = u_0(x), & 0 \leq x \leq 1 \\ \frac{\partial u}{\partial t}(x, 0) = u_1(x), & 0 \leq x \leq 1 \end{cases} [/math]

Este sistema describe completamente la evolución temporal de los desplazamientos transversales de una cuerda vibrante con los extremos fijos, dada su posición e impulso iniciales.

3.4 Resolución del sistema por separación de variables en términos de los coeficientes de Fourier de los datos iniciales

Primero, planteamos la solución de la ecuación de onda mediante separación de variables. Asumimos que la solución \( u(x,t) \) se puede escribir como un producto de dos funciones independientes, una que depende solo de \( x \) y otra que depende solo de \( t \):

[math]u(x, t) = X(x) T(t) [/math]

Sustituimos esta forma en la ecuación de onda:

[math]\frac{\partial^2 (X(x) T(t))}{\partial t^2} = \frac{\partial^2 (X(x) T(t))}{\partial x^2} [/math]

Esto se simplifica a:

[math]X(x) \frac{d^2 T(t)}{dt^2} = T(t) \frac{d^2 X(x)}{dx^2} [/math]

Dividimos ambos lados por [math]X(x) T(t) [/math]:

[math]\frac{1}{T(t)} \frac{d^2 T(t)}{dt^2} = \frac{1}{X(x)} \frac{d^2 X(x)}{dx^2} [/math]

Dado que el lado izquierdo depende solo de \( t \) y el lado derecho solo de \( x \), ambos deben ser iguales a una constante que llamaremos \(-\lambda\). Esto nos lleva a dos ecuaciones ordinarias:

[math]\frac{d^2 T(t)}{dt^2} + \lambda T(t) = 0 [/math] [math]\frac{d^2 X(x)}{dx^2} + \lambda X(x) = 0 [/math]

Para satisfacer las condiciones de contorno \( u(0, t) = 0 \) y \( u(1, t) = 0 \), tenemos:

[math]X(0) = 0 [/math] [math]X(1) = 0 [/math]

La ecuación para \( X(x) \) tiene la solución general:

[math]X(x) = A \sin(\sqrt{\lambda} x) + B \cos(\sqrt{\lambda} x) [/math]

Aplicando las condiciones de contorno, obtenemos que \( B = 0 \) y que \( \sqrt{\lambda} \) debe ser un múltiplo de \( \pi \):

[math] \sqrt{\lambda} = n \pi \quad \text{para} \quad n = 1, 2, 3, \ldots [/math]

Por lo tanto, \( \lambda = (n \pi)^2 \) y las funciones propias son:

[math] X_n(x) = \sin(n \pi x) [/math]

La ecuación para \( T(t) \) se convierte en:

[math] \frac{d^2 T_n(t)}{dt^2} + (n \pi)^2 T_n(t) = 0 [/math]

Con la solución general:

[math] T_n(t) = C_n \cos(n \pi t) + D_n \sin(n \pi t) [/math]

Combinando \( X(x) \) y \( T(t) \), la solución general para \( u(x,t) \) es:

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

Donde \( A_n \) y \( B_n \) son coeficientes a determinar a partir de las condiciones iniciales.

1. 1. 1. Determinación de los coeficientes \( A_n \) y \( B_n \)

Las condiciones iniciales son:

[math]u(x, 0) = u_0(x) [/math] [math] \frac{\partial u}{\partial t}(x, 0) = u_1(x) [/math]

Para \( t = 0 \):

[math] u(x, 0) = \sum_{n=1}^{\infty} A_n \sin(n \pi x) = u_0(x) [/math]

Usamos la ortogonalidad de las funciones seno para determinar \( A_n \):

[math] A_n = 2 \int_0^1 u_0(x) \sin(n \pi x) \, dx [/math]

Para la derivada temporal:

[math] \frac{\partial u}{\partial t}(x, 0) = \sum_{n=1}^{\infty} B_n n \pi \cos(n \pi \cdot 0) \sin(n \pi x) = u_1(x) [/math] [math] B_n n \pi = 2 \int_0^1 u_1(x) \sin(n \pi x) \, dx [/math]

Por lo tanto:

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

Finalmente, la solución de la ecuación de onda en términos de los coeficientes de Fourier de los datos iniciales es:

[math]u(x, t) = \sum_{n=1}^{\infty} \left( \left(2 \int_0^1 u_0(x) \sin(n \pi x) \, dx \right) \cos(n \pi t) + \left(\frac{2}{n \pi} \int_0^1 u_1(x) \sin(n \pi x) \, dx \right) \sin(n \pi t) \right) \sin(n \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) \).

3.4.1 Ejemplo de solución periódica en tiempo

En esta parte del documento para comprender lo anteriormente explicado de manera óptima se va a proceder a analizar una serie de ejemplo. 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]. Es por ello que la representación de la solución en el intervalo t ∈ [0, 2] es la siguiente:

centro

Para la obtención de esta representación se ha creado el siguiente código en Matlab:

% 
clear
close all
clc

% Definimos nuestros datos iniciales:
u_0=@(x) (exp(-100.*(x-1/2).^2));
u_1=@(x) (0);

% Definimos nuestro intervalo de tiempo t:
t=0:0.001:2;

% Definimos nuestro intervalo x:
x=0:0.001:1;

% Definimos el número de términos de la serie a dibujar:
n=50;

% Calculamos los términos de fourier asociados a u_0:
ak_0=zeros(n,1);
bk_0=zeros(n,1);
ck_0=zeros(n,1);

f_sen_0=@(x,k) u_0(x).*sin(k*pi*x);
f_cos_0=@(x,k) u_0(x).*cos(k*pi*x);
f_medio_0=@(x,k) u_0(x).*1/2;

for k=1:n
    % Vector con los valores de f_sin:
    y_sen_0=f_sen_0(x,k);
            
    % Vector con los valores de f_cos:
    y_cos_0=f_cos_0(x,k);
            
    % Vector con los valores de f_sin:
    y_medio_0=f_medio_0(x);

    ak_0(k)=trapz(x,y_sen_0);
    bk_0(k)=trapz(x,y_cos_0);
    ck_0(k)=trapz(x,y_medio_0);
end

% Calculamos los términos de fourier asociados a u_1:
ak_1=zeros(n,1);
bk_1=zeros(n,1);
ck_1=zeros(n,1);

f_sen_1=@(x,k) u_1(x).*sin(k*pi*x);
f_cos_1=@(x,k) u_1(x).*cos(k*pi*x);
f_medio_1=@(x,k) u_1(x).*1/2;

for k=1:n
    % Vector con los valores de f_sin:
    y_sen_1=f_sen_1(x,k);
            
    % Vector con los valores de f_cos:
    y_cos_1=f_cos_1(x,k);
            
    % Vector con los valores de f_sin:
    y_medio_1=f_medio_1(x);

    ak_1(k)=trapz(x,y_sen_1);
    bk_1(k)=trapz(x,y_cos_1);
    %ck_1(k)=trapz(x,y_medio_1);
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) + (ak_0(i).*cos(i*pi*tt)+bk_1(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 ecuación de ondas')
xlabel('x')
ylabel('t')

Como se observa en la imagen anterior la función muestra un comportamiento que se repite regularmente con el tiempo. A este tipo de funciones se les denomina periódicas en tiempo. Esto se puede observar de manera óptima para su comprensión en la siguiente gráfica, en la cual se observa el comportamiento explicado a anteriormente cuyo periodo es 2.

centro

Para la obtención de esta animación se ha hecho uso del siguiente cógido:

% 
%% Ejercicio 3 animación
clear
close all
clc

% Definimos nuestros datos iniciales:
u_0=@(x) (exp(-100.*(x-1/2).^2));
u_1=@(x) (0);

% Definimos nuestro intervalo de tiempo t:
t=0:0.001:2;

% Definimos nuestro intervalo x:
x=0:0.001:1;

% Definimos el número de términos de la serie a dibujar:
n=50;

% Calculamos los términos de fourier asociados a u_0:
ak_0=zeros(n,1);

f_sen_0=@(x,k) u_0(x).*sin(k*pi*x);

for k=1:n
    % Vector con los valores de f_sin:
    y_sen_0=f_sen_0(x,k);

    ak_0(k)=2.*trapz(x,y_sen_0);
end

% Calculamos los términos de fourier asociados a u_1:
bk_1=zeros(n,1);

f_cos_1=@(x,k) u_1(x).*cos(k*pi*x);

for k=1:n
            
    % Vector con los valores de f_cos:
    y_cos_1=f_cos_1(x,k);
            
    bk_1(k)=2/(k*pi)*trapz(x,y_cos_1);
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) + (ak_0(i).*cos(i*pi*tt)+bk_1(i).*sin(i*pi*tt)).*sin(i*pi*xx);
end

% Graficamos la solución:
% Crear un objeto de VideoWriter para el video
video = VideoWriter('onda_ejercicio_3.avi');
video.FrameRate = 100;
open(video);

% Crear la figura
figure

for i = 1:length(t)

    plot(x,u(x,t(i)))
    xlim([0,1])
    ylim([-1.5,1.5])
    xlabel('x')
    ylabel('u(x,t)')
    title(['Solución de la ecuación de onda en t=', num2str(t(i))])
    
    
    % Capturar el frame para el video
    frame = getframe(gcf);
    writeVideo(video, frame);
    
    % Limpiar la figura para la siguiente iteración, excepto los ejes
    clf; % Clear current figure
    
    
end


% Cerrar el objeto de video
close(video);

3.4.2 Ejemplo de la onda viajera

3.4.2.1 Condiciones Dirichlet

En este caso se va a tratar el ejemplo de una onda viajera en un solo sentido denotada por [math] u(x,t)=f(x-t)[/math]. Es por ello que se toman 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-\frac{1}{2})^2}[/math]

Se ha creado un código en MATLAB para obtener la siguiente respresentación en el mismo intervalo:

centro

Se observa que el tiempo de llegada de la onda de un extremo al otro es una unidad de tiempo. A su vez puesto que la cuerda mide una unidad de espacio se concluye que la velocidad de propagación es 1.

Para la obtención de esta representación se ha creado el siguiente código en Matlab:

% 

   %% Ejercicio 4
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;



% Cálculo de coeficientes de Fourier asociados a u_0
ak_0 = zeros(n,1);
for k = 1:n
    f_sen_0 = @(x) u_0(x).*sin(k*pi*x);

    ak_0(k) = 2*trapz(x, f_sen_0(x));
end

% Cálculo de coeficientes de Fourier asociados a u_1
bk_1 = zeros(n,1);
for k = 1:n
    f_sen_1 = @(x) u_1(x).*sin(k*pi*x);

    bk_1(k) =(2/(k*pi))*trapz(x, f_sen_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) + (ak_0(i)*cos(i*pi*tt) + bk_1(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 ecuación de ondas')
xlabel('x')
ylabel('t')

A continuación, al igual que en el apartado anterior, se dispone una representación del comportamiento de la onda en función del tiempo.

centro

Como se observa, cuando la onda alcanza la frontera la solución se anula la solución para todo tiempo.

% 
%% Ejercicio 4 animación
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;



% Cálculo de coeficientes de Fourier asociados a u_0
ak_0 = zeros(n,1);
for k = 1:n
    f_sen_0 = @(x) u_0(x).*sin(k*pi*x);

    ak_0(k) = 2*trapz(x, f_sen_0(x));
end

% Cálculo de coeficientes de Fourier asociados a u_1
bk_1 = zeros(n,1);
for k = 1:n
    f_sen_1 = @(x) u_1(x).*sin(k*pi*x);

    bk_1(k) =(2/(k*pi))*trapz(x, f_sen_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) + (ak_0(i)*cos(i*pi*tt) + bk_1(i)*sin(i*pi*tt)).*sin(i*pi*xx);
end

% Crear un objeto de VideoWriter para el video
video = VideoWriter('onda_ejercicio_4.avi');
video.FrameRate = 100;
open(video);

% Crear la figura
figure

for i = 1:length(t)

    plot(x,u(x,t(i)))
    xlim([0,1])
    ylim([-1.5,1.5])
    xlabel('x')
    ylabel('u(x,t)')
    title(['Solución de la ecuación de onda en t=', num2str(t(i))])
    
    
    % Capturar el frame para el video
    frame = getframe(gcf);
    writeVideo(video, frame);
    
    % Limpiar la figura para la siguiente iteración, excepto los ejes
    clf; % Clear current figure
    
    
end


% Cerrar el objeto de video
close(video);

3.4.2.2 Condiciones Neumann

A continuación, para comprender mejor la diferencia respecto a las condiciones fronteras, se va a resolver el mismo problema cambiando las condiciones Dirichlet a condiciones Neumann, es decir [math]u_0(0,t)=u_1(1,t)=0[/math]. De esta manera se pretende resolver el siguiente problema:

\begin{cases} 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{cases}

En este caso, aplicando separación de variables, obtenemos que la solución es:

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

siendo

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

y

[math]b_k = \frac{\int_{0}^{1} \cos(k\pi x) u_1(x) \, dx}{\int_{0}^{1} \cos^2(k\pi x) \, dx} [/math]

La representación obtenida mediante el código que se presenta posteriormente es la siguiente:

centro

%
%% Ejercicio 5
clear
close all
clc

% Definimos f(x)
f = @(x) exp(-100 * (x - 1/2).^2);

% -Derivada de f(x)
df = @(x) 200 * (x - 0.5) .* exp(-100 * (x - 0.5).^2);

% Definimos el número de terminos de la serie
n=50;

% Definimos el intervalo de tiempo y espacio
x = linspace(0, 1, 200);
t_values = linspace(0, 2, 200); % Evaluamos desde t=0 hasta t=2

% Calculamos los coeficientes de Fourier para Neumann
a_k =zeros(n,1);
b_k =zeros(n,1);

for i=1:n
    a_k(i) =2 * trapz(x, f(x) .* cos(i * pi * x));
    b_k(i) =2 / (i * pi) * trapz(x, df(x) .* cos(i * pi * x));
end
% Calculamos la solución de la ecuación de onda
u=@(xx,tt) 0;
for i=1:n
    u=@(xx,tt) u(xx,tt) + (a_k(i).*cos(i*pi*tt)+b_k(i).*sin(i*pi*tt)).*cos(i*pi*xx);
end

% Calculamos u(X, T) para todos los puntos de la malla
[X, T] = meshgrid(x, t_values);


% Representamos gráficamente
figure
surf(X, T, u(X,T), 'EdgeColor', 'flat')
title('Solución de la ecuación de onda con condiciones Neumann')
xlabel('x')
ylabel('t')
zlabel('u(x, t)')

En este caso, para su correcta visualización, se ha vuelto a crear la siguiente representación:

centro

Como se observa el periodo en este caso vuelve a ser 2 y puesto que los extremos no son fijos porque tenemos de datos las derivadas la onda varía en la frontera.

%
%% Ejercicio 5 animación
clear
close all
clc

% Definimos f(x)
f = @(x) exp(-100 * (x - 1/2).^2);

% -Derivada de f(x)
df = @(x) 200 * (x - 0.5) .* exp(-100 * (x - 0.5).^2);

% Definimos el número de terminos de la serie
n=50;

% Definimos el intervalo de tiempo y espacio
x = linspace(0, 1, 200);
t_values = linspace(0, 2, 200); % Evaluamos desde t=0 hasta t=2

% Calculamos los coeficientes de Fourier para Neumann
a_k =zeros(n,1);
b_k =zeros(n,1);

for i=1:n
    a_k(i) =2 * trapz(x, f(x) .* cos(i * pi * x));
    b_k(i) =2 / (i * pi) * trapz(x, df(x) .* cos(i * pi * x));
end
% Calculamos la solución de la ecuación de onda
u=@(xx,tt) 0;
for i=1:n
    u=@(xx,tt) u(xx,tt) + (a_k(i).*cos(i*pi*tt)+b_k(i).*sin(i*pi*tt)).*cos(i*pi*xx);
end

% Crear un objeto de VideoWriter para el video
video = VideoWriter('onda_condiciones_newman_ejercicio_5_buena.avi');
video.FrameRate = 100;
open(video);
% Graficamos la solución en diferentes tiempos
figure

for t = t_values
    u_values = u(x, t);
    plot(x, u_values)
    xlim([0, 1])
    ylim([-0.5, 2])
    title(['Solución de la ecuación de onda con condiciones Neumann en t=', num2str(t)])
    xlabel('x')
    ylabel('u(x, t)')
    % Capturar el frame para el video
    frame = getframe(gcf);
    writeVideo(video, frame);
    clf
end

% Cerrar el objeto de video
close(video);

4 Soluciones fundamentales en dimensiones 1, 2 y 3

La solución fundamental resuelve el sistema,

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

En esta parte del trabajo vamos a dibujar la solución fundamental de la ecuación de ondas en dimensiones 1, 2 y 3. Cuyas expresiones son las siguientes.

Dimensión [math]n=1[/math]: [math] K_1(x,t) =\frac{1}{2c}\textbf{1}_{[-ct,ct]}(x)=\frac{1}{2c}\textbf{1}_{[0,ct]}(|x|). [/math]

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)=\frac{1}{2\pi c \sqrt{c^2t^2-|x|^2}} \textbf{1}_{[0,ct]}(|x|). [/math]

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

Para dibujarlas, es primero necesario hacer unos cambios en los casos de dimensiones 2 y 3.

• En el caso de la dimensión 2, encontramos una singularidad en el denominador de la función de la solución fundamental. Esta ocurre cuando ∣x∣=ct, lo que resulta en una división por cero. Para evitar esta singularidad, se agrega un término de regularización ϵ en el denominador. Esto asegura que la función sea suave y bien comportada en todas partes, incluidos los puntos cercanos a ∣x∣=ct.

[math] K^{ε}_{2} (x,t) = \frac{1}{ε + 2\pi c \sqrt(c^2 t^2 - ∣x∣^2)} ⲭ_{B(0,ct)}(x) [/math]

• Para la representación del caso de dimensión 3 se ha sustituido la delta de Didac por su aproximación,

[math] \delta (s) ≈ \phi_k (s) = \sqrt(\frac{k}{\pi})e^{(-ks^2)} [/math]

Con estos cambios, representamos las 3 soluciones fundamentales observando que estas son radiales.

centro

centro

centro

Que las soluciones fundamentales sean "radiales" significa que son simétricas respecto al origen en coordenadas radiales. Dicha periodicidad es clara en las gráficas anteriores. Esto implica que el valor de la solución en un punto depende únicamente de la distancia radial desde el origen, y no de la dirección en la que se mida esa distancia.

El código empleado para esta representación ha sido el siguiente.

clear
close all
clc

% Dimensión 1
% Definimos los intervalos x,t:
x=-1:0.001:1;
t=0:0.001:1;

% Creación de la rejilla
[X,T] = meshgrid(x,t);

% Velocidad de propagación
c = 1;

% Cálculo de la solución fundamental
K = (X <= c*T) .* (X >= -c*T) / (2*c);

% Visualización de la solución fundamental
figure("Name",'Solución fundamental dimensión 1')
surf(X, T, K,'EdgeColor','interp');
title("Solución fundamental en dimensión 1");
xlabel("x");
ylabel("t");

% Dimensión 2
% Definimos nuestros intervalos
r=0:0.001:1;
t=0:0.001:1;

% Creamos la malla
[R,T]=meshgrid(r,t);

% Velocidad de propagación
c=1;

% Cálculo de la solución fundamental
K=(R<=c*T)./(eps + 2*pi*c.*sqrt(c^2*T.^2-R.^2));

figure("Name",'Solución fundamental dimensión 2')
surf(R,T,K,'EdgeColor','flat')
title("Solución fundamental de dimensión 2")
xlabel("r")
ylabel("t")

% Definimos nuestros intervalos
r=0:0.001:1;
t=0:0.001:1;


[R,T]=meshgrid(r,t);

c=1;                     % Velocidad de propagación 
k=1000;                  

% Calculamos la solución fundamental
K=sqrt(k/pi).*exp(-k * (R-c*T).^2)./(4*pi*c*R);

figure("Name",'Solución fundamental dimensión 3')
surf(R,T,K,'EdgeColor','flat')
title("Solución fundamental de dimensión 3")
xlabel("r")
ylabel("t")

4.1 Ejemplo en dimesión 2

Veamos a continuación el sistema que resuelve la solución fundamental en dimensión 2 en una bola de radio [math] \frac{1}{2}[/math].

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

En este caso, la solución fundamental viene dada por la convolución,

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

Dibujamos las soluciones para t = 0, 0.5, 1, 2.

centro

centro

centro

centro

A diferencia de lo visto en el caso general, aquí no podemos apreciar periodicidad. Sin embargo, el sistema se puede expresar en función de radio únicamente, es decir,

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

El código empleado es el siguiente,

PONE CÓDIGO