Ecuación de ondas. Otelo, Yan y Mika
Contenido
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:
Cuya solución encontrada gracias al método de separación de las variables, el principio de superposición y las condiciones iniciales es
donde:
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
Y el dibujo de la solución en el intervalo de tiempo [math] [0,2] [/math] queda así:
Podemos observar que solución para estas condiciones iniciales es periódica en el tiempo con periodo [math] T = 2 [/math] ya que en la gráfica se repite el patrón inicial.
Esto también se puede apreciar si repetimos el dibujo en el intervalo temporal [math] [0,4] [/math]:
1.2.1 Código
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:
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:
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:
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
En el video anterior se puede observar como tarda una unidad de tiempo en recorrer una unidad de espacio ya que el pico de la onda en t = 0 empieza en x = 0.5 y vuelve a x = 0.5 en t = 1. Sin embargo, una vez alcanzada la frontera de x = 1, la onda cambia de signo y de sentido. Este hecho vuelve a ocurrir cuando alcanza la próxima frontera x = 0, donde vuelve a tomar el signo positivo y el sentido inicial. Esto se repite cada vez que alcanza una frontera. Además, se puede observar que en t = 2, la onda vuelve a ser igual a la inicial, comprobando así que el periodo T es igual a 2.
1.3.1 Código
1.4 Diferencia con la condición Neumann
Por último, veremos qué ocurre si repetimos lo anterior cambiando las condiciones Dirichlet por Neumann:
Con ello, el sistema pasa a ser:
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:
donde, ahora:
En este caso, la gráfica de la solución queda así:
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.
1.4.1 Código
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
Donde [math] \delta (x) [/math] es la delta de Dirac. Formalmente, se define como el siguiente límite
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]:
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]:
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]:
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:
% 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
Su representación gráfica en función de la variable radial es la dada por la siguiente imagen:
% 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:
Así, su representación queda como sigue:
% 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
