Diferencia entre revisiones de «Ecuación de Ondas (CGomJRod)»

De MateWiki
Saltar a: navegación, buscar
(→‎Solución Fundamental)
(→‎Ejemplo en dimensión 2)
Línea 552: Línea 552:
 
<br/>
 
<br/>
  
−
[[Archivo:CGomJRodOndasEjConvTiempos.png|300px|thumb|right|Solución Estacionaria]]
+
[[Archivo:CGomJRodOndasEjConvTiempos.png|800px|thumb|right|Solución Estacionaria]]
 
{{matlab|codigo=
 
{{matlab|codigo=
 
clear all
 
clear all
Línea 616: Línea 616:
 
<br/>
 
<br/>
  
−
[[Archivo:CGomJRodOndasEjConvVideo.gif|600px|thumb|right|Solución Estacionaria]]
+
[[Archivo:CGomJRodOndasSolFundVideo.gif|600px|thumb|right|Solución Estacionaria]]
 
{{matlab|codigo=
 
{{matlab|codigo=
 
clear all
 
clear all

Revisión del 14:20 27 may 2024

1 Introducción

En el siguiente documento se tiene como objetivo estudiar las soluciones de la ecuación de ondas. En primer lugar, la estudiaremos en dimensión en el compacto [math][0,1][/math]. Además, comprobaremos las diferencias entre las soluciones al imponer diferentes condiciones frontera, concretamente Dirichlet y Neumann [math]\textbf{INTERPRETACIÓN?}[/math]. En segundo lugar, pasaremos al estudio de las soluciones fundamentales de esta ecuación en dimensiones 1, 2 y 3. Estas son soluciones de la ecuación de ondas cuando dejamos de resolverla en un conjunto acotado y tratamos de hallarla para todo [math]\mathbb{R}^n[/math]. Su importancia radica en que podemos calcular soluciones de otros problemas a partir de esta y mediante una convolución. En este trabajo nos centraremos únicamente en sacar las gráficas de estas funciones más que en su desarrollo teórico.

2 Conocimientos Previos

Como se ha comentado en la sección anterior, en este documento estudiaremos la ecuación de ondas, por lo que será necesario conocer ciertos resultados en relación con esta famosa ecuación. En primer lugar, aunque su nombre ya es muy representativo, veamos que modeliza para tener algo de intuición sobre esta. Imaginemos una cuerda vibrante que ocupa el intervalo [math][0, L][/math]. Consideremos que tiene una densidad [math]d[/math] y tensión [math]\tau[/math] constantes tales que la velocidad de propagación [math]c = \tau/d[/math]. Se puede demostrar que la ecuación que modeliza este comportamiento es la ecuación de ondas:


[math]u_{tt} – c^2 u_{xx} =0 [/math]

Esta se trata de una ecuación en derivadas parciales de segundo orden por lo que necesitaremos alguna información extra para poder resolver el problema que tenemos entre manos. Supongamos que tiene los extremos fijos, por lo que hay que añadirle las condiciones frontera de tipo Dirichlet [math]u(0,t)=u(L,t)=0[/math]. Finalmente, hay que tener en cuenta las condiciones iniciales del problema. Llamaremos [math]u_0(x)[/math] y [math]u_1(x)[/math] su posición e impulso iniciales respectivamente. Así, obtenemos el problema a resolver:

[math]\left \{ \begin{array}{ll} u_{tt} – c^2 u_{xx} =0 \\ u(0,t)=u(L,t)=0 \\ u(x,0)=u_0(x) \\ u_t(x,0)=u_1(x) \end{array} \right. [/math]

Este problema se puede resolver por separación de variables, y apoyándonos en el principio de superposición y en la teoría de series de Fourier [math]\textbf{LINK TRABAJO SERIES DE FOURIER}[/math] se puede demostrar el siguiente resultado:

[math]\textbf{Teorema (Existencia de soluciones):}[/math] Sea un problema como el anterior, entonces existe una solución de la forma:

[math] u(x,t)= \sum _{k=1}^{\infty} (\frac{L}{k\pi c} u_{1,k} sin(\frac{k \pi c}{L}t)+ u_{0,k}cos(\frac{k\pi c}{L}t)) sin(\frac{k\pi}{L}x) [/math]

donde [math]u_{0,k}[/math] y [math]u_{1,k}[/math] son los coeficientes de Fourier que se calculan como sigue:

[math] u_{0,k}= \frac{\int_0^L sin (k \pi x) u_0 (x) dx}{\int_0 ^L sin^2(k \pi x) dx}[/math]


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

Otro resultado importante que debemos conocer previo al estudio que haremos en este trabajo es la fórmula de D'Alembert. Esta nos da una expresión de la solución del problema de Cauchy global, es decir que el abierto sobre el que trabajamos no es acotado; todo el espacio. Utilizaremos su expresión para dimensión 1.

[math]\textbf{Teorema (Fórmula de D'Alembert):}[/math] Sea el problema de Cauchy global:

[math]\left \{ \begin{array}{ll} u_{tt} – c^2 u_{xx} =0 \\ u(x,0)=u_0(x) \\ u_t(x,0)=u_1(x) \end{array} \right. [/math]

Entonces la solución de la ecuación viene dada por la expresión:

[math]u(x,t)=\frac{1}{2}(u_0(x-ct)+u_0(x+ct))+\frac{1}{2c}\int_{x-ct}^{x+ct}u_1(y)dy[/math]

A partir de esta expresión es posible obtener una cierta intuición sobre el funcionamiento de las soluciones de esta ecuación. Como se puede ver en la fórmula anterior, para calcular [math]u(x,t)[/math] solo necesitamos conocer los valores de [math]u_0[/math] en los puntos [math]x\pm ct[/math] y los de [math]u_1[/math] en [math][x-ct,x+ct][/math]. Es decir, la solución solo depende de cómo es la condición inicial en dicho intervalo. A esto se le conoce como dominio de dependencia de [math]u(x,t)[/math]. Por otro lado, podemos obtener otra interpretación que veremos más adelante pero cuya intuición se puede ver en el primer sumando de la fórmula; la solución está compuesta por dos ondas iguales que se mueven en la misma dirección y sentido contrario a la misma velocidad.

[math]\textbf{Principio de Huygens:}[/math] El principio de Huygens fue formulado originalmente por el físico y matemático holandés Christiaan Huygens en 1678. Nos dice que todo punto de un frente de onda inicial puede considerarse como una fuente de ondas esféricas secundarias que se extienden en todas las direcciones con la misma velocidad, frecuencia y longitud de onda que el frente de onda del que proceden.

3 Ejemplos de resolución

3.1 Condiciones Dirichlet

Una vez hemos visto muy por encima algunos resultados y conceptos teóricos previos, podemos pasar a resolver algunos ejercicios a modo de ejemplo. Supongamos una cuerda vibrante como la anterior con los extremos fijos y donde [math]L=1[/math]. Además, como condición inicial vamos a imponer [math]u_0(x)=e^{-100(x-\frac{1}{2})^2}[/math] y [math]u_1(x)=0[/math]. Es decir, queremos resolver el siguiente problema:

[math]\left \{ \begin{array}{ll} u_{tt} – u_{xx} =0 \\ u(0,t)=u(1,t)=0 \\ u(x,0)=e^{-100(x-\frac{1}{2})^2} \\ u_t(x,0)=0 \end{array} \right. [/math]


Aplicando el teorema de existencia de la sección anterior, obtenemos la siguiente solución:

[math]\textbf{METER FORMULONCIO JAVI}[/math]


A continuación se muestra la gráfica de la solución para [math]t \in [0,2][/math] tomando los 50 primeros términos de la serie.

Solución Estacionaria
clear all
close all
%Definicion de intervalos espacial y temporal
xini=0; xfin=1; nx=150;
tini=0; tfin=2; nt=150;
xx=linspace(xini,xfin,nx)'; tt=linspace(tini,tfin,nt)';
%Numero de terminos de la serie de Fourier
nterms=50;
%Condiciones iniciales
u0=@(x)(exp(-100*((x-1/2).^2)));
u1=@(x)(zeros(height(x),width(x)));
%Definicion de la discretizacion para la regla del trapecio
nfour=200; vectint=linspace(xini,xfin,nfour)';
%Serie de Fourier
U=zeros(nt,nx);
for k=1:nterms
    u0k=trapz(vectint,u0(vectint).*sin(k*pi*vectint))/trapz(vectint,sin(k*pi*vectint).^2);
    u1k=trapz(vectint,u1(vectint).*sin(k*pi*vectint))/trapz(vectint,sin(k*pi*vectint).^2);
    U=U+(u0k*cos(k*pi*tt)+u1k*sin(k*pi*tt)/(k*pi))*sin(k*pi*xx)';
end
%Graficado de la solucion
[XX,TT]=meshgrid(xx,tt);
mesh(XX,TT,U)
zlim([-2,2])
xlabel('x',Interpreter='latex')
ylabel('t',Interpreter='latex')
saveas(gcf,'CGomJRodOndasEjer3.png')



En ella se aprecia como aparecen dos ondas que se mueven en sentidos opuestos hasta que al llegar al límite del intervalo rebotan y ``cambian su orientación. Estas dos se vuelven a unir en el centro y posteriormente continúan su camino hasta que vuelven a rebotar cambiando de nuevo la orientación y encontrándose en el centro, empezando de nuevo el proceso; la solución es periódica en tiempo. Esto se aprecia mejor en el siguiente GIF

Solución Estacionaria
clear all
close all
%Definicion de intervalos espacial y temporal
xini=0; xfin=1; nx=150;
tini=0; tfin=2; nt=150;
xx=linspace(xini,xfin,nx)'; tt=linspace(tini,tfin,nt)';
%Numero de terminos de la serie de Fourier
nterms=50;
%Condiciones iniciales
u0=@(x)(exp(-100*((x-1/2).^2)));
u1=@(x)(zeros(height(x),width(x)));
%Definicion de la discretizacion para la regla del trapecio
nfour=200; vectint=linspace(xini,xfin,nfour)';
%Serie de Fourier
U=zeros(nt,nx);
for k=1:nterms
    u0k=trapz(vectint,u0(vectint).*sin(k*pi*vectint))/trapz(vectint,sin(k*pi*vectint).^2);
    u1k=trapz(vectint,u1(vectint).*sin(k*pi*vectint))/trapz(vectint,sin(k*pi*vectint).^2);
    U=U+(u0k*cos(k*pi*tt)+u1k*sin(k*pi*tt)/(k*pi))*sin(k*pi*xx)';
end
%Hacer el video
pelicula=VideoWriter('CGomJRodOndasEjer3Video.avi');
pelicula.FrameRate=20;
open(pelicula);
for i=1:nt
    plot(xx,U(i,:)','b')
    ylim([-2,2])
    xlabel('x',Interpreter='latex')
    title('$t='+string(tt(i))+'$',Interpreter="latex")
    figura=figure(1);
    imagen=getframe(figura);
    writeVideo(pelicula,imagen)
end
close(pelicula);



Hemos visto que nada más empezar el proceso, la onda se divide en dos y se mueven en sentidos opuestos. Sin embargo, podemos eliminar este efecto cambiando las condiciones iniciales del problema a estudiar. Si consideramos [math]u_0=f(x)[/math] y [math]u_1(x)=-f'(x)[/math] con [math]e^{-100(x-\frac{1}{2})^2}[/math], podemos comprobar a continuación que el comportamiento que muestra la solución anterior es el mismo salvo que en este caso no encontramos dos ondas, sino una.

Solución Estacionaria
clear all
close all
%Definicion de intervalos espacial y temporal
xini=0; xfin=1; nx=150;
tini=0; tfin=2; nt=150;
xx=linspace(xini,xfin,nx)'; tt=linspace(tini,tfin,nt)';
%Numero de terminos de la serie de Fourier
nterms=50;
%Condiciones iniciales
u0=@(x)(exp(-100*((x-1/2).^2)));
u1=@(x)(200*(x-1/2).*exp(-100*((x-1/2).^2)));
%Definicion de la discretizacion para la regla del trapecio
nfour=200; vectint=linspace(xini,xfin,nfour)';
%Serie de Fourier
U=zeros(nt,nx);
for k=1:nterms
    u0k=trapz(vectint,u0(vectint).*sin(k*pi*vectint))/trapz(vectint,sin(k*pi*vectint).^2);
    u1k=trapz(vectint,u1(vectint).*sin(k*pi*vectint))/trapz(vectint,sin(k*pi*vectint).^2);
    U=U+(u0k*cos(k*pi*tt)+u1k*sin(k*pi*tt)/(k*pi))*sin(k*pi*xx)';
end
%Graficado de la solucion
[XX,TT]=meshgrid(xx,tt);
mesh(XX,TT,U)
zlim([-2,2])
xlabel('x',Interpreter='latex')
ylabel('t',Interpreter='latex')
saveas(gcf,'CGomJRodOndasEjer4.png')



Solución Estacionaria
clear all
close all
%Definicion de intervalos espacial y temporal
xini=0; xfin=1; nx=150;
tini=0; tfin=2; nt=150;
xx=linspace(xini,xfin,nx)'; tt=linspace(tini,tfin,nt)';
%Numero de terminos de la serie de Fourier
nterms=50;
%Condiciones iniciales
u0=@(x)(exp(-100*((x-1/2).^2)));
u1=@(x)(200*(x-1/2).*exp(-100*((x-1/2).^2)));
%Definicion de la discretizacion para la regla del trapecio
nfour=200; vectint=linspace(xini,xfin,nfour)';
%Serie de Fourier
U=zeros(nt,nx);
for k=1:nterms
    u0k=trapz(vectint,u0(vectint).*sin(k*pi*vectint))/trapz(vectint,sin(k*pi*vectint).^2);
    u1k=trapz(vectint,u1(vectint).*sin(k*pi*vectint))/trapz(vectint,sin(k*pi*vectint).^2);
    U=U+(u0k*cos(k*pi*tt)+u1k*sin(k*pi*tt)/(k*pi))*sin(k*pi*xx)';
end
%Hacer el video
pelicula=VideoWriter('CGomJRodOndasEjer4Video.avi');
pelicula.FrameRate=20;
open(pelicula);
for i=1:nt
    plot(xx,U(i,:)','b')
    ylim([-2,2])
    xlabel('x',Interpreter="latex")
    title('$t='+string(tt(i))+'$',Interpreter="latex")
    figura=figure(1);
    imagen=getframe(figura);
    writeVideo(pelicula,imagen)
end
close(pelicula);



La pregunta que surge ahora es por qué pasa esto, ya que parece que las nuevas condiciones escogidas salen de la nada. Sin embargo si recordamos la fórmula de D'Alembert que presentamos en [math]\textbf{INSERTAR Link intro}[/math] y hacemos los siguientes cálculos podemos observar que estas condiciones no son tan caprichosas y en realidad obedecen a lo que dicta esta fórmula. Sustituyendo en esta las condiciones del nuevo problema y haciendo uso del teorema fundamental del cálculo:

[math]u(x,t)= = \frac{1}{2}(f(x-ct)+f(x+ct))-\frac{1}{2c}\int_{x-ct}^{x+ct}f'(y)dy=\frac{1}{2}(f(x-ct)+f(x+ct))+\frac{1}{2c} (cf(x-ct)-cf(x+ct))=f(x-ct)[/math]

Particularizando a este caso, como [math]c=1[/math] :

[math]u(x,t)= f(x-t)[/math]

Es decir, que la solución de la ecuación en realidad es una onda transversal con la forma de la gráfica de [math]f[/math]. Este es el motivo por el cual al imponer estas condiciones obtenemos una de las ondas del caso anterior. Nótese que si en vez de imponer [math]u_1(x)=-f'(x)[/math], imponemos [math]u_1(x)=f'(x)[/math] la onda que obtenemos es la que se mueve en sentido contrario. También es necesario recalcar que en el primer caso, al dividirse la onda en dos, como la energía ha de conservarse las dos ondas que obtenemos son más "pequeñas"; se aplanan. Esto es porque la energía debe dividirse entre ambas y por partes iguales, no debería haber una preferencia por las ondas que viajan en uno u otro sentido. Por otro lado, si observamos el GIF del segundo caso se aprecia como esto no sucede, sino que simplemente hay una traslación de la gráfica de la función al no dividirse. Finalmente, una cosa que queda por verificar es la velocidad de propagación de la onda, que como indicamos debe ser [math]c=1[/math]. Para verlo, hemos tomado dos imágenes de la onda en diferentes tiempos, espaciados en 0,2 unidades de tiempo, y las hemos superpuesto. En esta imagen que se muestra a continuación podemos observar como la gráfica se ha desplazado 0,2 unidades en el eje espacial. Por tanto, la velocidad es efectivamente [math]c=1[/math].

Solución Estacionaria
clear all
close all
%Definicion de intervalos espacial y temporal
xini=0; xfin=1; nx=50;
tini=0; tfin=2; nt=50;
xx=linspace(xini,xfin,nx)'; tt=linspace(tini,tfin,nt)';
%Numero de terminos de la serie de Fourier
nterms=50;
%Condiciones iniciales
u0=@(x)(exp(-100*((x-1/2).^2)));
u1=@(x)(200*(x-1/2).*exp(-100*((x-1/2).^2)));
%Definicion de la discretizacion para la regla del trapecio
nfour=200; vectint=linspace(xini,xfin,nfour)';
%Serie de Fourier
U=zeros(nt,nx);
for k=1:nterms
    u0k=trapz(vectint,u0(vectint).*sin(k*pi*vectint))/trapz(vectint,sin(k*pi*vectint).^2);
    u1k=trapz(vectint,u1(vectint).*sin(k*pi*vectint))/trapz(vectint,sin(k*pi*vectint).^2);
    U=U+(u0k*cos(k*pi*tt)+u1k*sin(k*pi*tt)/(k*pi))*sin(k*pi*xx)';
end
%Graficado de la solucion
plot(xx,U(1,:),'b')
hold on
plot(xx,U(5,:),'c')
quiver(xx(21),U(1,21),tt(5),0,1,'r')
xlabel('x',Interpreter='latex')
ylim([-0,1.5])
legend(['Solucion en $t='+string(tt(1))+'$','Solucion en $t='+string(tt(5))+'$'],Interpreter="latex")
saveas(gcf,'CGomJRodOndasEjer41.png')
hold off



3.2 Condiciones Neumann

Una vez hemos entendido qué pasaba en el ejemplo anterior, vamos a considerar el mismo problema de antes pero con condiciones Neumann. Sea el siguiente problema:

[math]\left \{ \begin{array}{ll} 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)=0 \end{array} \right. [/math]

A continuación mostramos las mismas imágenes que antes pero para este nuevo problema:

Solución Estacionaria
clear all
close all
%Definicion de intervalos espacial y temporal
xini=0; xfin=1; nx=150;
tini=0; tfin=2; nt=150;
xx=linspace(xini,xfin,nx)'; tt=linspace(tini,tfin,nt)';
%Numero de terminos de la serie de Fourier
nterms=50;
%Condiciones iniciales
u0=@(x)(exp(-100*((x-1/2).^2)));
u1=@(x)(200*(x-1/2).*exp(-100*((x-1/2).^2)));
%Definicion de la discretizacion para la regla del trapecio
nfour=200; vectint=linspace(xini,xfin,nfour)';
%Serie de Fourier
U=zeros(nt,nx);
for k=1:nterms
    u0k=trapz(vectint,u0(vectint).*cos(k*pi*vectint))/trapz(vectint,cos(k*pi*vectint).^2);
    u1k=trapz(vectint,u1(vectint).*cos(k*pi*vectint))/trapz(vectint,cos(k*pi*vectint).^2);
    U=U+(u0k*cos(k*pi*tt)+u1k*sin(k*pi*tt)/(k*pi))*cos(k*pi*xx)';
end
%Graficado de la solucion
[XX,TT]=meshgrid(xx,tt);
mesh(XX,TT,U)
xlabel('x',Interpreter='latex')
ylabel('t',Interpreter='latex')
saveas(gcf,'CGomJRodOndasEjer5.png')



Solución Estacionaria
clear all
close all
%Definicion de intervalos espacial y temporal
xini=0; xfin=1; nx=150;
tini=0; tfin=2; nt=150;
xx=linspace(xini,xfin,nx)'; tt=linspace(tini,tfin,nt)';
%Numero de terminos de la serie de Fourier
nterms=50;
%Condiciones iniciales
u0=@(x)(exp(-100*((x-1/2).^2)));
u1=@(x)(200*(x-1/2).*exp(-100*((x-1/2).^2)));
%Definicion de la discretizacion para la regla del trapecio
nfour=200; vectint=linspace(xini,xfin,nfour)';
U=zeros(nt,nx);
for k=1:nterms
    u0k=trapz(vectint,u0(vectint).*cos(k*pi*vectint))/trapz(vectint,cos(k*pi*vectint).^2);
    u1k=trapz(vectint,u1(vectint).*cos(k*pi*vectint))/trapz(vectint,cos(k*pi*vectint).^2);
    U=U+(u0k*cos(k*pi*tt)+u1k*sin(k*pi*tt)/(k*pi))*cos(k*pi*xx)';
end
%Hacer el video
pelicula=VideoWriter('CGomJRodOndasEjer5Video.avi');
pelicula.FrameRate=20;
open(pelicula);
for i=1:nt
    plot(xx,U(i,:)','b')
    ylim([-2,2])
    xlabel('x',Interpreter='latex')
    title('$t='+string(tt(i))+'$',Interpreter="latex")
    figura=figure(1);
    imagen=getframe(figura);
    writeVideo(pelicula,imagen)
end
close(pelicula);



[math]\textbf{NOTA. EL COMENTARIO ESTÁ EN MAYUSCULAS PARA QUE SE VEA BIEN, NO TIENE COMO OBJETO OFENDER: MENTIRA, NO SE DIVIDE EN DOS, ADEMÁS LAS CONDICIONES INICIALES ESTÁN MAL PUESTAS, LO DE QUE REBOTA SI QUE ESTÁ BIEN. SOLO NOS QUEDAMOS CON LO SEGUNDO QUE HAS PUESTO.}[/math] En estas se puede ver un comportamiento similar al anterior, vemos una onda que se divide en dos que después de rebotar con los extremos del intervalo terminan combinándose en el centro de nuevo para volver a empezar. Sin embargo, ahora encontramos dos diferencias fundamentales respecto a la situación anterior: la onda ni cambia de orientación ni cambia de sentido inmediatamente. Al cambiar de sentido esta parece acumularse en los extremos hasta que llegado un punto termina por cambiar de sentido. Por otro lado, si modificamos las condiciones iniciales como hicimos en el apartado anterior obtenemos el mismo resultado. En vez de tener dos ondas más pequeñas que se mueven en sentidos contrarios, aparece una única onda que se mueve en un sentido. Nótese que el comportamiento de esta nueva onda es el mismo que en el caso inmediatamente anterior, pues no hemos cambiado las condiciones frontera. A continuación se muestran las imágenes correspondientes a este caso en las que se aprecia lo que acabamos de comentar :


INSERTAR IMÁGENES (CARECE DE SENTIDO)

Por tanto, vemos que cambiar las condiciones de frontera lo único que afecta al comportamiento de la onda es en realidad cómo "rebota". Antes, cambiando su orientación y sentido de forma inmediata, y ahora, acumulándose hasta que termina modificando su trayectoria. Por lo que la pregunta que surge de manera natural es por qué pasa esto. A diferencia de lo anterior, la cuerda que estaríamos no se encuentra fija, sino que es móvil. Podríamos imaginar que los extremos se encuentran atados a un objeto que se mueve por una vía sin resistencia en el eje vertical. Esto da lugar a que efectivamente la tangente a la cuerda en el extremo tenga siempre pendiente nula. Esto además explicaría el porqué no se invierte la orientación.

Imaginemos ahora la onda desplazándose hacia uno de los extremos, según se va acercando este va ganando altura poco a poco. Por conservación de la energía es lógico pensar que la altura máxima alcanzada por el extremo debería ser igual al máximo de la gráfica. Una vez lo alcanza, la onda "rebota" y cambia de sentido para repetir el proceso. Sin embargo, la altura alcanzada es mayor que la comentada. Esto se puede deber a que al llegar ahí se produce un fenómeno parecido a la resonancia. Cuando el frente de la onda golpea el extremo del intervalo este cambia de sentido y se encuentra con la parte que aún sigue desplazándose en el sentido anterior. Al tener la misma fase, ambas se suman provocando que el extremo de la cuerda termine superando la altura máxima de la onda inicial. Es precisamente esto lo que da la sensación de que la onda se acumula en esa zona.

4 Solución Fundamental

Como hemos comentado en [math]\textbf{Link Introducción}[/math] ahora trataremos de dibujar las gráficas de la solución fundamental [math]K_n[/math] para dimensión 1,2 y 3. Estas soluciones se obtienen al dar un impulso inicial localizado en [math]x=0[/math]. Matemáticamente es la solución del problema:

[math]\left \{ \begin{array}{ll} u_{tt} – \Delta u =0 \\ u(x,0)=0 \\ u_t(x,0)=\delta(x) \end{array} \right. [/math]

donde [math]x \in \mathbb{R}^n[/math] y [math]\delta (x)[/math] la delta de Dirac.

Particularizando para las distintas dimensiones se puede demostrar que:

[math]\textbf{dim=1}[/math]

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

donde [math]H[/math] es la función de Heavyside:

[math]H(x)=\left \{ \begin{array}{ll} 0 \hspace{0.3cm}\text{si}\hspace{0.1cm} x\lt0 \\ 1 \hspace{0.3cm} \text{si}\hspace{0.1cm}x \geq 0 \end{array} \right. [/math]

[math]\textbf{dim=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].

[math]\textbf{dim=3}[/math]

[math]K_3(x,t)=\frac{\delta(|x|-ct)}{4\pi c|x|}[/math]
Solución Estacionaria
close all
clear all
%Definicion de intervalos espacial y temporal
tini=0; tfin=2; nt=50;
xini=-2; xfin=2; nx=50;
tt=linspace(tini,tfin,nt)'; xx=linspace(xini,xfin,nx)';
%Definición de la solución fundamental
K1=@(x,t)((heaviside(x+t)-heaviside(x-t))/2);
%Graficado de la solucion
[XX,TT]=meshgrid(xx,tt);
mesh(XX,TT,K1(XX,TT))
xlabel('x',Interpreter='latex')
ylabel('t',Interpreter='latex')
saveas(gcf,'CGomJRodOndasSolFund1.png')



Solución Estacionaria
clear all
close all
%Definicion de intervalos espacial y temporal
e=0.01;
tini=0; tfin=2; nt=500;
rini=-2; rfin=2; nr=500;
tt=linspace(tini,tfin,nt)'; rr=linspace(rini,rfin,nr)';
%Construcción de la solución fundamental
U=zeros(nt,nr);
K2=@(e,r,t)((caract(r,t)./(e+2*pi*sqrt(t.^2-r.^2))))';
for i=1:length(tt)
    U(i,:)=K2(e,rr,tt(i));
end
%Graficado de la solucion
[RR,TT]=meshgrid(rr,tt);
mesh(RR,TT,real(U))
xlabel('r',Interpreter='latex')
ylabel('t',Interpreter='latex')
saveas(gcf,'CGomJRodOndasSolFund2.png')



Solución Estacionaria
clear all
close all
%Función Aproximación de la delta
phi=@(k,s)(sqrt(k/pi)*exp(-k*s.^2));
%Definicion de intervalos espacial y temporal
tini=0; tfin=1; nt=500;
rini=-2; rfin=2; nr=500;
tt=linspace(tini,tfin,nt)'; rr=linspace(rini,rfin,nr)';
%Constante de la función phi
k=1000;
%Definición de la solución fundamental
K3=@(k,r,t)(phi(k,r-t))./(4*pi*r);
%Graficado de la solucion
[RR,TT]=meshgrid(rr,tt);
mesh(RR,TT,K3(k,RR,TT))
xlabel('r',Interpreter='latex')
ylabel('t',Interpreter='latex')
saveas(gcf,'CGomJRodOndasSolFund3.png')



4.1 Ejemplo en dimensión 2

Solución Estacionaria
clear all
close all
%Definición de la solución fundamental en polares preparada para la
%integral
% e es epsylon, (rx,thx) son las coordenadas de x en polares, (ry,thy) son
% las coordenadas de y en polares, como argumento x de la función siempre
% se introduce x-y. Por último t es el t de la definición de la función
K2=@(e,rx,ry,thx,thy,t)(caract(sqrt(rx.^2+ry.^2-2.*rx.*ry.*cos(thy-thx)),t)./(e+2*pi*sqrt(t.^2-rx.^2-ry.^2+2.*rx.*ry.*cos(thy-thx))));
%Definicion de intervalos espacial y temporal
tini=0; tfin=2; nt=101;
rini=-2; rfin=2; nr=101;
tt=linspace(tini,tfin,nt)'; rr=linspace(rini,rfin,nr)';
%Epsylon de la solución fundamental
e=0.01;

%Construcción de la solución por convolución
U=zeros(nr,nt);
%Discretización radial y angular para la integral
nrint=200; rrint=linspace(0,1/2,nrint); 
nthint=200; ththint=linspace(0,2*pi,nthint);
[RRINT,THTHINT]=meshgrid(rrint,ththint);
for i=1:nr
    for j=1:nt
        %Convolución de la solución fundamental con la condición inicial
        U(i,j)=trapz(ththint,trapz(rrint,K2(e,rr(i),RRINT,0,THTHINT,tt(j)),2));
    end
end
U=real(U');
%Graficado de la solución (ejes radial y temporal)
thini=0; thfin=2*pi; nth=50; thth=linspace(thini,thfin,nth);
[RR,TT]=meshgrid(rr,tt);
mesh(RR,TT,U)
xlabel('r',Interpreter='latex')
ylabel('t',Interpreter='latex')
saveas(gcf,'CGomJRodOndasEjConvRadial.png')



Solución Estacionaria
clear all
close all
%Definición de la solución fundamental en polares preparada para la
%integral
% e es epsylon, (rx,thx) son las coordenadas de x en polares, (ry,thy) son
% las coordenadas de y en polares, como argumento x de la función siempre
% se introduce x-y. Por último t es el t de la definición de la función
K2=@(e,rx,ry,thx,thy,t)(caract(sqrt(rx.^2+ry.^2-2.*rx.*ry.*cos(thy-thx)),t)./(e+2*pi*sqrt(t.^2-rx.^2-ry.^2+2.*rx.*ry.*cos(thy-thx))));
%Definicion de intervalos espacial y temporal
tini=0; tfin=2; nt=101;
rini=-2; rfin=2; nr=101;
tt=linspace(tini,tfin,nt)'; rr=linspace(rini,rfin,nr)';
%Epsylon de la solución fundamental
e=0.01;

%Construcción de la solución por convolución
U=zeros(nr,nt);
%Discretización radial y angular para la integral
nrint=200; rrint=linspace(0,1/2,nrint); 
nthint=200; ththint=linspace(0,2*pi,nthint);
[RRINT,THTHINT]=meshgrid(rrint,ththint);
for i=1:nr
    for j=1:nt
        %Convolución de la solución fundamental con la condición inicial
        U(i,j)=trapz(ththint,trapz(rrint,K2(e,rr(i),RRINT,0,THTHINT,tt(j)),2));
    end
end
U=real(U');
%Paso a coordenadas cartesianas
thini=0; thfin=2*pi; nth=50; thth=linspace(thini,thfin,nth);
[RR,THTH]=meshgrid(rr,thth);
XX=RR.*cos(THTH);
YY=RR.*sin(THTH);
%Graficado de la solución (ejes cartesianos, distintos tiempos)
subplot(1,4,1)
mesh(XX,YY,ones(nth,1)*real(U(find(tt==0),:)))
zlim([0,3])
xlabel('x',Interpreter='latex')
ylabel('y',Interpreter='latex')
subtitle('$t=0$',Interpreter='latex')
subplot(1,4,2)
mesh(XX,YY,ones(nth,1)*real(U(find(tt==0.5),:)))
zlim([0,3])
xlabel('x',Interpreter='latex')
ylabel('y',Interpreter='latex')
subtitle('$t=0.5$',Interpreter='latex')
subplot(1,4,3)
mesh(XX,YY,ones(nth,1)*real(U(find(tt==1),:)))
zlim([0,3])
xlabel('x',Interpreter='latex')
ylabel('y',Interpreter='latex')
subtitle('$t=1$',Interpreter='latex')
subplot(1,4,4)
mesh(XX,YY,ones(nth,1)*real(U(find(tt==2),:)))
zlim([0,3])
xlabel('x',Interpreter='latex')
ylabel('y',Interpreter='latex')
subtitle('$t=2$',Interpreter='latex')



Solución Estacionaria
clear all
close all
%Definición de la solución fundamental en polares preparada para la
%integral
% e es epsylon, (rx,thx) son las coordenadas de x en polares, (ry,thy) son
% las coordenadas de y en polares, como argumento x de la función siempre
% se introduce x-y. Por último t es el t de la definición de la función
K2=@(e,rx,ry,thx,thy,t)(caract(sqrt(rx.^2+ry.^2-2.*rx.*ry.*cos(thy-thx)),t)./(e+2*pi*sqrt(t.^2-rx.^2-ry.^2+2.*rx.*ry.*cos(thy-thx))));
%Definicion de intervalos espacial y temporal
tini=0; tfin=2; nt=101;
rini=-2; rfin=2; nr=101;
tt=linspace(tini,tfin,nt)'; rr=linspace(rini,rfin,nr)';
%Epsylon de la solución fundamental
e=0.01;

%Construcción de la solución por convolución
U=zeros(nr,nt);
%Discretización radial y angular para la integral
nrint=200; rrint=linspace(0,1/2,nrint); 
nthint=200; ththint=linspace(0,2*pi,nthint);
[RRINT,THTHINT]=meshgrid(rrint,ththint);
for i=1:nr
    for j=1:nt
        %Convolución de la solución fundamental con la condición inicial
        U(i,j)=trapz(ththint,trapz(rrint,K2(e,rr(i),RRINT,0,THTHINT,tt(j)),2));
    end
end
U=real(U');
%Paso a coordenadas cartesianas
thini=0; thfin=2*pi; nth=50; thth=linspace(thini,thfin,nth);
[RR,THTH]=meshgrid(rr,thth);
XX=RR.*cos(THTH);
YY=RR.*sin(THTH);
%Hacer el video
pelicula=VideoWriter('CGomJRodOndasEjConvVideo.avi');
pelicula.FrameRate=20;
open(pelicula);
for i=1:nt
    mesh(XX,YY,ones(nth,1)*real(U(i,:)))
    zlim([0,3])
    xlabel('x',Interpreter='latex')
    ylabel('y',Interpreter='latex')
    title('$t='+string(tt(i))+'$',Interpreter="latex")
    figura=figure(1);
    imagen=getframe(figura);
    writeVideo(pelicula,imagen)
end
close(pelicula);



Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Ecuación_de_Ondas_(CGomJRod)&oldid=72958»