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í:
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]:
Veamos en un GIF cómo se comporta la onda en función del 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:
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:
Y además hemos creado un GIF para visualizar la onda respecto al tiempo:
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:
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
2.2 Solución del problema en dimensión 2
La solución del problema:
viene dada por la convolución:
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].
A continuación se muestran los dibujos de las soluciones para [math] t =0, 0.5, 1, 2[/math]
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 es constantemente cero. 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 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.
