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

Observamos que en [math]t=0[/math] la onda se divide en dos, que viajan hacia cada uno de los extremos espaciales. Cuando llegan a los extremos, cambian de signo y de sentido para volver a [math] x=0.5[/math] y cuando se reencuentran mantienen el signo para repetir el recorrido, pero esta vez, cuando llegan a los extremos vuelven a hacerse positivas y en [math]t=2[/math] se vuelve al estado inicial. Por ello, la solución para estas condiciones iniciales es periódica en el tiempo con periodo [math] T = 2 [/math].

Esto también se puede apreciar si repetimos el dibujo en el intervalo temporal [math] [0,4] [/math]:

Solución ecuación de ondas

Veamos en un GIF cómo se comporta la onda en función del tiempo:

Ondas según el tiempo

1.2.1 Código

clear all
close all
%Resumen: la idea del programa es crear una matriz de pares (tiempo,
%espacio) y en cada uno de esos puntos calcular la aproximación de la
%solución mediante un bucle for que realiza la suma de términos de la serie
%calculada

%creamos la malla en  la que introduciremos las evaluaciones 
evaluaciones2 = zeros(100,100); 
evaluacionesvideo = zeros(1000,100); 
%quedan impuestas los condiciones frontera por ser una matriz de ceros ya
% que posteriormente no modificamos ni la primera ni la última columna

%limite temporal y vectores de tiempos y espacios
T = 2;
t = linspace(0,T,100); 
x = linspace(0,1,100); 
tp = linspace(0,10,1000);

%funcion de la condicion inicial
cond =@(y) exp(-100.*(y-1/2).^2); 

%imponemos la condicion inicial 
evaluaciones2(1,:) = cond(x)';
evaluacionesvideo(1,:) = cond(x)';

%en el siguiente bucle for, para cada k calculamos el coeficiente de
%Fourier y sumamos  el correspondiente término de la serie en cada punto de
% la malla. Para el cálculo de la integral que nos da el coeficiente de
% Fourier usamos la regla del trapecio con 1000 puntos en [0,1]
for k=1:50 
    %puntos
    N=1000;  
    %extremos del intervalo
    a=0; b=1;            
    h=(b-a)/N; 
    %partición del intervalo
    u=a:h:b;           
    %f es la funcion que queremos integrar
    f=(cond(u).*sin(k.*pi.*u))'; 
    %vector de pesos
    w=ones(N+1,1);                 
    w(1)=1/2; w(N+1)=1/2;
    %m es el valor de la integral
    m=h*w'*f; 
    %calculo final de los coeficientes del seno 
    coef = 2*m;   
    %sumamos los 100 primeros términos de la serie en cada punto de la
    % malla con los nuevos coeficientes
    for j =2:99
        for i = 2:100
                evaluaciones2(i,j) = evaluaciones2(i,j) + coef.*sin(k.*pi.*(x(j))).*cos(k*pi*t(i)); 
        end
        for i = 2:1000
                evaluacionesvideo(i,j) = evaluacionesvideo(i,j) + coef.*sin(k.*pi.*(x(j))).*cos(k*pi*tp(i)); 
        end
    end
end

%plot en 3D
figure(1)
surf(t,x,evaluaciones2') 
xlabel('t')
ylabel('x')
zlabel('u(x,t)')

% Crear la animación y guardarla como GIF
% Nombre del archivo GIF
filename = 'Ejercicio63.gif';
 
% Configuración de la animación
fps = 5; % Fotogramas por segundo
delayTime = 1 / fps; % Tiempo de retraso entre fotogramas

for i = 1:10:1000 % Reducir el número de fotogramas
    figure(2);
    plot(x,evaluacionesvideo(i,:))
    axis([0 1 -1.2 1.2])
    xlabel("Espacio")
    ylabel("u(x,t)")
    title("t = " + mat2str(tp(i)))
     
    % Guardar cada frame como imagen en formato GIF
    frame = getframe(gcf);
    im = frame2im(frame);
    [imind,cm] = rgb2ind(im,64); % Reducir el número de colores
     
    % Escribir en el archivo GIF
    if i == 1
        imwrite(imind,cm,filename,'gif', 'Loopcount',inf, 'DelayTime', delayTime);
    else
        imwrite(imind,cm,filename,'gif','WriteMode','append', 'DelayTime', delayTime);
    end
end

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 con longitud de onda [math]\lambda = 2 [/math] metros. Por tanto la velocidad es [math]c = \frac{\lambda}{T} = \frac{2}{2} = 1 m/s [/math]

La onda parte de [math]x = 0.5 [/math]m y se dirige hacia [math]x=1 [/math]. En [math]t = 0.5 [/math]s llega a esa frontera y allí cambia de signo y de dirección, para dirigirse hacia el extremo [math]x=0[/math], que tarda en alcanzar [math]1 [/math]s. Por último, cuando llega a este extremo, pasa de nuevo a tomar valores positivos y toma la dirección inicial para volver a la posición de partida [math]x = 0.5 [/math]. De nuevo probando. que el periodo es de [math]2 [/math] segundos.

Todo esto se ve mejor observando la onda desde otro ángulo:

Onda que viaja en un solo sentido

Y además hemos creado un GIF para visualizar la onda respecto al tiempo:

Onda que viaja en un solo sentido

1.3.1 Código

clear all 
close all


% Resumen: la idea del programa es crear una matriz de pares (tiempo,
% espacio) y en cada uno de esos puntos calcular la aproximación de la
% solución mediante un bucle for que realiza la suma de términos de la serie
% calculada

% Creamos la malla en  la que introduciremos las evaluaciones 
evaluaciones2 = zeros(1000,100); 

% Límite temporal y vectores de tiempos y espacios
T = 2;
t = linspace(0,T,100); 
x = linspace(0,1,100); 
tp = linspace(0,10,1000);

% Función de la condicion inicial
cond0 =@(y) exp(-100*((y-(1/2)).^2)); 
cond1 = @(y) 200.*(y-(1/2)).*exp(-100*((y-(1/2)).^2));

% Imponemos la condicion inicial 
evaluaciones2(1,:) = cond0(x)';

% En el siguiente bucle for, para cada k calculamos el coeficiente de
% Fourier y sumamos  el correspondiente término de la serie en cada punto de
% la malla. Para el cálculo de la integral que nos da el coeficiente de
% Fourier usamos la regla del trapecio con 1000 puntos en [0,1]
for k=1:100 
    % Puntos
    N=1000;  
    % Extremos del intervalo
    a=0; b=1;            
    h=(b-a)/N; 
    % Partición del intervalo
    u=a:h:b;           
    % f es la funcion que queremos integrar
    f0 = (cond0(u).*sin(k.*pi.*u))';
    f1 = (cond1(u).*sin(k.*pi.*u))';
    % Vector de pesos
    w=ones(N+1,1);                 
    w(1)=1/2; w(N+1)=1/2;
    % m es el valor de la integral
    m0 = h*w'*f0; 
    m1 = h*w'*f1;
    % Calculo final de los coeficientes del seno 
    coef0 = 2*m0;
    coef1 = 2*m1;
    % Sumamos los 100 primeros términos de la serie en cada punto de la
    % malla con los nuevos coeficientes
    for j =1:99
        for i = 2:1000
                evaluaciones2(i,j) = evaluaciones2(i,j) + coef0.*sin(k.*pi.*(x(j))).*cos(k.*pi.*(tp(i))) + ...
                (1/(k*pi)).*coef1.*sin(k.*pi.*(x(j))).*sin(k.*pi.*(tp(i)));
        end
    end
end

% Crear la animación y guardarla como GIF
% Nombre del archivo GIF
filename = 'Ejercicio4.gif';

% Configuración de la animación
fps = 5; % Fotogramas por segundo
delayTime = 1 / fps; % Tiempo de retraso entre fotogramas

for i = 1:10:1000 % Reducir el número de fotogramas
    figure(1);
    plot(x,evaluaciones2(i,:))
    axis([0 1 -1.2 1.2])
    xlabel("Espacio")
    ylabel("u(x,t)")
    title("t = " + mat2str(tp(i)))
    
    % Guardar cada frame como imagen en formato GIF
    frame = getframe(gcf);
    im = frame2im(frame);
    [imind,cm] = rgb2ind(im,64); % Reducir el número de colores
    
    % Escribir en el archivo GIF
    if i == 1
        imwrite(imind,cm,filename,'gif', 'Loopcount',inf, 'DelayTime', delayTime);
    else
        imwrite(imind,cm,filename,'gif','WriteMode','append', 'DelayTime', delayTime);
    end
end

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

En este caso, la velocidad de la onda es la misma y el desplazamiento ocurre en el mismo sentido y dirección. Lo único que cambia es que cuando la onda llega a los extremos del intervalo espacial, en vez de cambiar de signo, aumenta hasta el máximo y regresa al mismo valor con el cambio de dirección. Esto hace que en los extremos se cumpla la condición Neumann que hemos impuesto.

Veamos en un GIF cómo se comporta la onda en función del tiempo:

Onda con condiciones Neumann

1.4.1 Código

clear all
close all
%Resumen: la idea del programa es crear una matriz de pares (tiempo,
%espacio) y en cada uno de esos puntos calcular la aproximación de la
%solución mediante un bucle for que realiza la suma de términos de la serie
%calculada

%creamos la malla en  la que introduciremos las evaluaciones 
evaluaciones2 = zeros(100,100);
evaluacionesvideo = zeros(1000,100);
%quedan impuestas los condiciones frontera por ser una matriz de ceros ya
% que posteriormente no modificamos ni la primera ni la última columna

%limite temporal y vectores de tiempos y espacios
T = 2;
t = linspace(0,T,100); 
x = linspace(0,1,100); 
tp = linspace(0,10,1000);

%f y f'
cond =@(y) exp(-100.*(y-1/2).^2);
cond2 =@(y) (200.*y-100).*exp(-100.*(y-1/2).^2);


%en el siguiente bucle for, para cada k calculamos el coeficiente de
%Fourier y sumamos  el correspondiente término de la serie en cada punto de
% la malla. Para el cálculo de la integral que nos da el coeficiente de
% Fourier usamos la regla del trapecio con 1000 puntos en [0,1]
for k=1:50 
    %puntos
    N=1000;  
    %extremos del intervalo
    a=0; b=1;            
    h=(b-a)/N; 
    %partición del intervalo
    u=a:h:b;           
    %f es la funcion que queremos integrar
    f=(cond(u).*cos(k.*pi.*u))'; 
    g=(cond2(u).*cos(k.*pi.*u))';
    %vector de pesos
    w=ones(N+1,1);                 
    w(1)=1/2; w(N+1)=1/2;
    %m es el valor de la integral
    m=h*w'*f; 
    n=h*w'*g;
    %calculo final de los coeficientes del seno 
    coef = 2*m;   
    coef2 = 2*n;
    %sumamos los 100 primeros términos de la serie en cada punto de la
    % malla con los nuevos coeficientes
    for j =1:100
        for i = 1:100
                evaluaciones2(i,j) = evaluaciones2(i,j) + (coef.*cos(k*pi*t(i)) + coef2/(k*pi).*sin(k*pi*t(i))).*cos(k.*pi.*(x(j))); 
        end
        for i = 1:1000
                evaluacionesvideo(i,j) = evaluacionesvideo(i,j) + (coef.*cos(k*pi*tp(i)) + coef2/(k*pi).*sin(k*pi*tp(i))).*cos(k.*pi.*(x(j))); 
        end
    end
end

%plot en 3D
figure(1)
surf(t,x,evaluaciones2') 
xlabel('t')
ylabel('x')
zlabel('u(x,t)')

% Crear la animación y guardarla como GIF
% Nombre del archivo GIF
filename = 'Ejercicio65.gif';
 
% Configuración de la animación
fps = 5; % Fotogramas por segundo
delayTime = 1 / fps; % Tiempo de retraso entre fotogramas

for i = 1:10:1000 % Reducir el número de fotogramas
    figure(2);
    plot(x,evaluacionesvideo(i,:))
    axis([0 1 -1.2 1.2])
    xlabel("Espacio")
    ylabel("u(x,t)")
    title("t = " + mat2str(tp(i)))
     
    % Guardar cada frame como imagen en formato GIF
    frame = getframe(gcf);
    im = frame2im(frame);
    [imind,cm] = rgb2ind(im,64); % Reducir el número de colores
     
    % Escribir en el archivo GIF
    if i == 1
        imwrite(imind,cm,filename,'gif', 'Loopcount',inf, 'DelayTime', delayTime);
    else
        imwrite(imind,cm,filename,'gif','WriteMode','append', 'DelayTime', delayTime);
    end
end

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} \hspace{32mm} x\in \mathbb{R}^n, t \gt0\\ 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 porque 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

2.2 Solución del problema en dimensión 2

La solución del problema:

[math] \left\{ \begin{array}{ll} u_{tt}-c^2\Delta u=0, \hspace{6mm} \hspace{32mm} x\in \mathbb{R}^2, t \gt0\\ u(x,0)=u_0(x), \hspace{3mm} u_t(x,0)= h(x) = \chi_{B(0,1/2)}(x)\hspace{6mm} x\in \mathbb{R}^2 . \\ \end{array} \right. [/math]

viene dada por la convolución:

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

Aprovecharemos que la solución es radial para calcularla sólo en los puntos de la forma [math](x_1,0)[/math] con [math]x_1 \gt 0[/math].Además, para evitar la singularidad de la solución fundamental volveremos a tomar el [math]epsilon\lt\math\gt en denominador. A continuación se muestran los dibujos de las soluciones para \ltmath\gt t =0, 0.5, 1, 2[/math]

Solución a partir de la convolución

t = [0;0.5;1;2];                 % Vector de tiempos donde evaluamos la solución
x = linspace(0,3,100);          % Vector de discretización del espacio
sols = zeros(length(t),length(x));    % Matriz de valores de la solución en la discretización en los diferentes tiempos

% Cálculo de valores de la solución y representación gráfica
figure(1)
for j=1:length(t)   
    for i = 1:length(x)
        sols(j,i) = g(x(i),t(j));
    end
    subplot(2,2,j)
    plot(x,sols(j,:))
    title("Solución para t =" + num2str(t(j)))
    xlabel('x1')
    ylabel('u((x1,0),t)')
end

% Función que devuelve los valores de la solución por convolución
function valor = g(x,t)
    % Definir la función característica de la bola de centro 0 y radio ct
    c = 1;
    chi_B = @(y) double(abs(y) <= c .* t); % Función característica en la bola de centro 0 y radio ct
    chi = @(y) double(abs(y) <= 1/2); % Función característica en la bola de centro 0 y radio 1/2
    % Definir la función K2(x,t)
    K2 = @(y) (1 ./ (0.01 + 2 .* pi .* c .* sqrt(c.^2 * t.^2 - (y.^2)))) .* chi_B(y);  
    valor = integral( @(y) (K2(x-y).*(chi(y))) ,0,1/2);  % Valor de la integral
end

Nótese que la representación gráfica es en 2D puesto que al ser la solución radial hemos representado los puntos de la forma (x1,0) por lo que tanto la segunda coordenada como el tiempo están fijados. Además, podemos observar que dependiendo del tiempo, los valores de la solución son constantemente cero a partir de distintos puntos. Esta dependencia se debe a las funciones características, pues la característica de la solución fundamental tiene un radio que depende del tiempo a diferencia de la característica de la convolución que tiene radio fijo igual a 1/2. En definitiva, la solución será constantemente igual a cero fuera de la bola de radio ct+1/2. Excepto cuando t es igual a cero que la convolución es igual a cero ya que el integrando es igual a cero. También, cabe destacar que el pico de la solución se alcanza cuando el radio es igual a ct, pues la solución fundamental toma el mayor valor posible cuando el denominador es lo menor posible, o sea, en los puntos de radio ct.