Matlab Ecuacion de ondas

[irg1994]

Bueno, antes de nada pido perdón porque no sabía donde colocar este hilo, y estoy seguro que me equivoque de lugar. Ruego que no lo elimineis porfa =S.

Os comento, estoy en primero de fisica, y para la semana que viene tengo que terminar un trabajo de computación en el cual me pedian dibujar algo en 3d con funciones matemáticas (no tenia porque ser muy complicado).

El caso esque a mi se me ocurrio hacer una animación (claramente con mi nivel, será muy sencilla y muy alejada de la realidad ) en la que se viera caer una gota de agua y las ondas que produce.

Entonces haciendo el programa, logre encontrar una ecuacion de onda que me produce más o menos algo que me puede servir, pero no consigo entenderla, y de nada me sirve hacer un programa que no entiendo =S. No le veo la relación con la ecuación sencilla que hemos dado hasta el momento de "y(x,t) = A sen(wt-kx+phi)".

Os pongo el script del programa y os lo explico:

% Formación de ondas en el agua.
clear all
clc

% Creamos una retícula.
x=-50:0.5:50;
y=-50:0.5:50;
[X,Y] = meshgrid(x,y);
% Radio de la Onda
R=sqrt(X.^2+Y.^2);
% Numero de onda (k = 2*pi/longitud de onda)
k=0.1;
% Fase inicial
phi=0;
count=1

% Introducimos la ecuación de la onda y la dibujamos.

for freq=0.1:0.1:2*pi

% Ecuación de la onda:
Z=sin((2*pi-freq)*k*R+phi);

surf(X,Y,Z,'Facecolor','blue','Edgecolor','none');
% Usamos Facecolor para seleccionar el color de nuestra
% onda, en nuestro caso azul. Y usamos Edgecolor para
% el color de la retícula de la onda, en nuestro caso
% usamos none para que no haya.

axis equal;
% Mediante este comando logramos igualar los tamaños
% de los ejes.

% Sombra y luz para visualizar mejor la onda.
camlight right;
lighting phong;

% Para crear esta animación en vez de usar el comando movie, uso
% getframe, que automaticamente convierte las imagenes creadas en
% nuestro bucle for agrupandolas en una imagen en movimieto.
S(count)=getframe;
count=count+1;
end


Bueno, ese es el programa, y la ecuación que no entiendo es esta:

Z=sin((2*pi-freq)*k*R+phi); donde freq toma los valores des 0.1 a 2pi, el valor que he dado a k = 0.1, R es la distancia al centro de mi plano y phi es una fase.
No entiendo la relacion con la ecuación que yo conocía =S

Haber si vosotros la veis y me podéis ayudar.

Muchas gracias y un saludo!

- - - Actualizado - - -

Nadie entiende o consigue descifrar que es cada cosa en esta ecuacion? Z=sin((2*pi-freq)*k*R+phi)

[arivasm]

Re: Ecuacion de ondas

Se trata de un modelo para una onda de amplitud 1 y de la forma $Z=\sin(\vec k\cdot\vec r-\omega t+\varphi)$. Lo que sucede es que el autor de esa línea no se ha "comido mucho el coco" y ha apañado las cosas de manera que 2*pi*k representaría el módulo del vector de propagación $\vec k$, y ha metido la evolución temporal de esa manera tan poco elegante.

Ya puestos, estaría mejor usar $\dst Z=A\sin\left[2\pi\left(\frac{r}{\lambda}-\frac{t}{T}\right)+\varphi\right]$, donde el usuario debería introducir los valores para la longitud de onda y el período. Eso sí, como la fase inicial es arbitraria, pues depende de la elección del t=0, yo la eliminaría. Por otra parte, si quieres garantizar la periodicidad, el tiempo lo modelaría por fracciones del período T.

Y siendo aún más perfeccionistas, la amplitud no debería ser constante, sino que debería depender de la distancia al punto de impacto. Aunque rigurosamente en un modelo ideal debería ser $\dst A=\frac{C}{r}$, donde C es una constante que depende de la energía aportada al medio, tendrías el problema del infinito en r=0, que puedes evitar con $\dst A=\frac{C}{r+\delta}$ donde $\delta$ es una cantidad que introduces "truculentamente" para que el programa no casque. Eso sí, el resultado te quedará aún más bonito que el que has dibujado.

[irg1994]

Re: Ecuacion de ondas

Ok, lo hare con esa ecuación mejor que si que la entiendo jaja gracias!

Pero aun así tengo unas cuantas dudas =S

¿La longitud de onda y el periodo, al estar dibujando el caso de una gota que cae al agua, no deberian ser mejor unos valores en vez de inputs?

Respecto a lo de la amplitud, ¿Como se cual esa costante aportada por el medio, no la encuentro =S? y la fase esa para cuando la distancia es 0, ¿Tendría que tantear cual sería la que me sirve mejor no?

Por ultimo, la variable t en funcion del periodo, no se si habria alguna forma mejor de introducirla a como lo he hecho, ya que la onda me hace algo raro, como si fuera muy lenta o algo nose, y no se si sera por eso o no le influye para nada =S

Bueno, pongo el script, muchas gracias de antemano! :




% Formación de ondas en el agua.
clear all
clc

% Creamos una retícula.
x=-50:0.5:50;
y=-50:0.5:50;
[X,Y] = meshgrid(x,y);
% Radio de la Onda
r=sqrt(X.^2+Y.^2);
% Amplitud
A=?

N=length(r);
count = 1;
% Inputs:

long = input('Introducir la longitud de onda: ');
T = input('Introducir el periodo de la onda:');

for t = linspace(0,T,N);
Z = A.*sin(2*pi.*((r./long)-(t/T)));
surf(X,Y,Z,'Facecolor','blue','Edgecolor','none');


% Mediante este comando logramos igualar los tamaños
% de los ejes.

axis equal;

% Sombra y luz para visualizar mejor la onda.
camlight right;
lighting phong;

% Para crear esta animación en vez de usar el comando movie, uso
% getframe, que automaticamente convierte las imagenes creadas en
% nuestro bucle for agrupandolas en una imagen en movimieto.
M(count)=getframe;
count=count+1;
end

- - - Actualizado - - -

Ademas, yo quiero que la onda se vaya disipando hasta no haber nada, y con este script no lo consigo claramente =S jaja

- - - Actualizado - - -

Bueno, que pare ya lo he conseguido haciendo esto:

for t = linspace(0,T,N);
for freq = 0.1:0.1:2*pi
Z = sin((2*pi-freq).*((r./long)-(t/T)));


...

Pero me falta el problema de que no se como meter la amplitud y el tema de si pongo fijo los valores de la longitud y el periodo??

Ya que una vez que tenga la onda me gustaria intentar meter una esfera cayendo antes de producirse la onda sin complicarme mucho, haber si lo consigo =S

[arivasm]

Re: Ecuacion de ondas

Entiendo que no pretendes modelar un golpe real, sino simplemente hacer una animación que la recuerde. Por tanto, yo tomaría la longitud de onda como input. Respecto del período, lo que te sugería es que los valores de t fuesen tomando valores, por ejemplo, 0, 0.02 T, 0.04 T, etc. Es decir, t/T incrementándose según cierta cantidad, que conviene que "cierre" para mostrar bien el aspecto cíclico.

La fase inicial puedes elegirla. Es más, si quieres que comience en un máximo para r=0 y t=0 toma pi/2 (pues usas senos).

Las constantes de la amplitud fíjalas tú, de tal manera que la figura no se te salga fuera de escala (C) y que represente adecuadamente el resultado sin tener problemas con números enormes para r cerca de cero ($\delta$).

La esfera cayendo quizá no sea buena idea, porque en un golpe real tienes una onda de choque que avanza que te complicará mucho la simulación.

[irg1994]

Re: Ecuacion de ondas

Ok, pero sigo sin entender como me intentas decir que escriba lo de los valores de t =S

Y loo de (2*pi-freq) que añadí para que la onda se disipase lo debería quitar?

- - - Actualizado - - -

Y un ejemplo de amplitud que podría ser A = 16./(r+8)?? no se cuales son los ordenes de magnitud o que debo escoger mas o menos

[arivasm]

Re: Ecuacion de ondas

Antes de nada, yo no sé Mathlab (pero sí programo en otros lenguajes). Por supuesto podrás hacer un bucle for tal que cierta variable (p. ej. la llamaré it) recorra valores desde 0 hasta N-1, siendo N el número de intervalos que tomarás para el tiempo a lo largo de un período. De esa manera t/T es lo mismo que it/N (¡ojo, porque tiene que ser una división de coma flotante, no entera!) con lo que la parte sinusoidal de la onda se codificará como

$\dst\sin\left[2\pi\left(\frac{r}{\lambda}-\frac{t}{T}\right)+\varphi\right]=\sin\left[2\pi\left(\frac{r}{\lambda}-\frac{i_t}{N}\right)+\varphi\right]$

Por cierto que el problema de la división entera no lo tendrás si codificas adecuadamente

$\dst\sin\left(\frac{2\pi r}{\lambda}-\frac{2\pi i_t}{N}+\varphi\right)$

Sobre la disipación, en realidad no hay tal cosa (o es prácticamente despreciable). La atenuación de la amplitud de la onda se debe a que su energía se distribuye a lo largo de una circunferencia cada vez más grande.

Ahora que lo comento, me acabo de dar cuenta de que la amplitud no decaerá con r (como sucederá en un caso 3D), sino con su raíz cuadrada, pues la intensidad se distribuye según $I=P/2\pi r$, con lo que la amplitud, que es proporcional a la raíz cuadrada de la intensidad, deberá seguir el patrón

$\dst A(r)=\frac{C}{\sqrt{r+\delta}}$

(con la $\delta$ por el problema de los infinitos).

Para la elección de C y $\delta$ haz pruebas, pues dependerá de la representación. La C es arbitraria, en el sentido de que refleja la potencia del impacto que causa la onda, pero te limitará la imagen: si es demasiado grande la altura de las ondas será excesiva y no se verá nada, si es demasiado pequeña tendrás una onda muy cutre. Para $\delta$ comienza usando un valor pequeño (en comparación con la longitud de onda, por ejemplo la centésima parte de ésta); si es demasiado pequeño tendrás un pico central enorme.

Terminaré ampliando lo que te dije antes. En un golpe real, de un objeto pequeño contra un líquido inicialmente en reposo, tienes una onda de choque que se propaga, que hará que la amplitud no sólo dependa de la distancia sino también del tiempo. Si te quieres meter en ese follón deberías manejar para la amplitud no una C, sino una C(r,t) que representará la amplitud de dicha onda de choque.

Yo de ti, primero codificaría el tren de ondas que se obtiene con C=cte, y que equivaldrá a algo así como que haya en el origen una varilla oscilante perturbando constantemente el medio. Una vez que lo tengas resuelto puedes animarte a meter en danza una C dependiente de la distancia y del tiempo que represente un pulso que avanza desde r=0. Aunque quizá un perfil gaussiano sería una buena descripción, si no fuese que no es cómodo de codificar, se puede buscar algo que se parezca pero que no sea exactamente igual. Por ejemplo, se podría probar con algo así como $\dst C(r,t)=C_0\ee^{-\alpha|r-vt|}$ donde habría que escoger adecuadamente el C_0 (debería valer el mismo C que se use antes) y los parámetros $\alpha$ y $v$.

[irg1994]

Re: Ecuacion de ondas

puff, en menudo berenjenal me he metido =S jajajaja En si yo solo tenia que dibujar una función en 3d algo curiosa y elaborada para computacion I de primero de fisica, y en buena hora se me ocurrio este ejemplo xD

Intento comprender y asimilar todo lo que me dices, pero me resulta complicado, y mas aun trasladarlo esto a un programa que conocí hace 2 meses xD, creo que tratare de ir probando cosas en el poco tiempo que me quede a fin de conseguir algo que al menos no sea ridiculo xD

Muchas gracias igualmente!

- - - Actualizado - - -

Haber pongo el script que tengo, que ahora mismo lo unico que hace es que te pide la longitud el periodo y la fase inicial, y he definido mediante condicionales la distancia minima para que no falle la amplitud. El caso esque esto, a mi me a costado bastante hacerlo porque soy muy malo jaja pero cualquiera podria hacerlo en dos minutos, entonces es un poco cutre presentar esto, y no se que hacer =S

% Ondas en el agua
clear all
clc

%Inputs

long = input('Introducir la longitud de onda: ');
T = input('Introducir el periodo de la onda: ');
Tiempo= 15*T; % Tiempo que estara oscilando.
phi = input('Introducir fase incial: ');
count = 1; % Esta variable es para usar el getframe para hacer la animación


% Creamos una retícula.
x=(-10*long):0.5:10*(long);
y=(-10*long):0.5:10*(long);
[X,Y] = meshgrid(x,y);

% Radio de la Onda
r=(X.^2+Y.^2).^0.5;

% Definimos la amplitud de la onda como A = C/sqrt(r), donde C es una
% costante dependiente de la energía aportada por el medio y r es la
% distancia al foco. Por ello para evitar que la amplitud se nos vaya a
% infinito es necesario tomar un valor de la amplitud para distancias muy
% cercanas. Yo tome C=1 por lo que cuando la distancia al foco sea menor
% que 0.8, diremos que la Amplitud es 1.
r0=0.8;
n = length(r);




% Programa:


for t = 0:0.1:Tiempo
% Mediante estos dos bucles y el condicional distinguimos los dos casos
% de la distancia al foco para determinar la Amplitud de la onda.
for i = 1:n
for j = 1:n
if (r(i,j))>r0
A = 1./(r(i,j).^0.5);
% Ecuación de la onda:
Z(i,j) = A.*(sin(2*pi.*((r(i,j)./long)-(t/T))));
else
Z(i,j) = sin(2*pi.*((r(i,j)./long)-(t/T)));
end
end
end

surf(X,Y,Z,'Facecolor','blue','Edgecolor','none');

% Mediante este comando logramos igualar los tamaños
% de los ejes.
axis equal;

% Sombra y luz para visualizar mejor la onda.
camlight right;
lighting phong;

% Para crear esta animación en vez de usar el comando movie, uso
% getframe, que automaticamente convierte las imagenes creadas en
% nuestro bucle "for" agrupandolas en una imagen en movimieto.

M(count)=getframe;
count=count+1;

end

[natanael]

Re: Ecuacion de ondas

¡Qué bonita! ¿Sería posible hacer una versión con un par de adaptaciones? En primer lugar, como es un caso bidimensional, hacer que la amplitud decaiga no con r^2, sino con r, es decir, $C=1/\sqrt{x^2+y^2}$ debido a que la potencia se distribuye uniformemente a lo largo de una línea de tamaño proporcional al radio ($2\pi r$), en vez de una superficie proporcional al cuadrado del mismo ($4\pi r^2$), como sucede en el caso 3D; en segundo, modificar levemente esta última expresión para evitar los infinitos, de manera que sea $C=1/\sqrt{x^2+y^2+\delta}$, donde el $\delta$ es una cantidad pequeña que se añade de manera forzada para evitar el infinito en el origen. Lo digo porque un compañero del foro había intentado algo similar a lo que has hecho (en este hilo), pero me temo que con poco éxito.

Añado: como no se trata de marear a nuestro amigo Julián, la petición es para ese otro hilo, no para éste.

si lo hago así como dices me sale esto:



si utlizo el factor ese $\dst C(r,t)=C_0\ee^{-\alpha|r-vt|}$ me sale:



y mira porque tardé tanto en responderte, hay una función, llamada "Función unitaria de Heaviside" (se vá a convertir en una de mis favoritas )

$H(x) = u(x)=\begin{cases} 0, & x < 0 \\ 1, & x > 0 \end{cases}$

pero es muy brusca debido a que se define como

$H(x)=\lim_{k \rightarrow \infty}\frac{1}{2}(1+\tanh kx)=\lim_{k \rightarrow \infty}\frac{1}{1+\mathrm{e}^{-2kx}}$

estonces la suavizo haciendo $k=20$ y entonces multipliqué esa función por la ultima que tiene el factor $\dst C(r,t)=C_0\ee^{-\alpha|r-vt|}$ y quedó algo así, siendo $R=\sqrt{x^2+y^2}$

$\frac{3}{1+{e}^{-2\times 20\left(\lambda\,T\,t-R +\frac{1}{2} \right) } } \frac{{e}^{-0.5\left\vert{}R-\lambda \,T\,t \right\vert{} } }{R+2} \sin \left[2\pi \left( \frac{R}{\lambda } -\frac{t}{T} \right) +\phi \right]$

con $\lambda =\frac{3}{2}$, $T=1$ y $\phi =0$ y $t$ bailando entre 0 y 2, aquí lo que estoy haciendo es limitar la oscilación de la onda a la región contenida dentro de un circulo cuyo radio crece con la velocidad del primer frente de onda, con todo esto el Mathematica escupió esto:



lo más increible es que la misma función, "Función unitaria de Heaviside", me sirvió para hacer algo más divertido, sea $R_1=\sqrt{{\left( x+2 \right)}^2+{\left( y+2 \right)}^2 }$ y $R_2=\sqrt{{\left( x-2 \right)}^2+{\left( y-2 \right)}^2 }$ entonces:

$\frac{3}{1+{e}^{-2\times 20\left(\lambda\,T\,t-R_1 +\frac{1}{2} \right) } } \frac{{e}^{-0.5\left\vert{}R_1-\lambda \,T\,t \right\vert{} } }{R_1+2} \sin \left[2\pi \left( \frac{R_1}{\lambda } -\frac{t}{T} \right) +\phi \right] + ...$
$..+\frac{1}{1+{e}^{-2\times 20\left( t-a \right) } } \frac{3}{1+{e}^{-2\times 20\left(\lambda\,T\left( t-a \right) -R_2 +\frac{1}{2} \right) } } \frac{{e}^{-0.5\left\vert{}R_2-\lambda \,T\left( t-a \right) \right\vert{} } }{R_2+2} \sin \left[2\pi \left( \frac{R_2}{\lambda } -\frac{t-a}{T} \right) +\phi \right]$
con $\lambda =\frac{3}{2}$, $T=a=1$ y $\phi =0$ y $t$ bailando entre 0 y 5, es decir, estoy obligando a que la segunda gota exista a partir de un tiempo determinado, con todo esto el Mathematica escupió esto:



Despues lo pasé al Matlab 7.8.0

clear all;
clc;
x=-4:0.1:4;
y=-4:0.1:4;
[X,Y] = meshgrid(x,y);
lambda=1.5;
T=1;
phi=0.;
R1=sqrt((X-2).^2+(Y+2).^2);
R2=sqrt((X+2).^2+(Y-2).^2);
a=2;
count=1;
set(gca,'nextplot','replacechildren','visible','off');


fig=figure;
aviobj = avifile('example.avi')

for t=-1:0.1:10
Z=(3./(1+exp(-2.*20.*(lambda*T*t-R1+0.5)))).*((exp(-0.5*abs(R1-lambda*T*t)))./(R1+2.)).*(sin(2.*pi*((R1./lambda)-(t./T))+phi))+(1./(1+exp(-2*20*(t-a)))).*(3./(1+exp(-2.*20.*(lambda*T*(t-a)-R2+0.5)))).*((exp(-0.5*abs(R2-lambda*T*(t-a))))./(R2+2.)).*(sin(2*pi*((R2./lambda)-((t-a)./T))+phi));
%surf(X,Y,Z,'Facecolor','blue','Edgecolor','none');
surf(X,Y,Z,double(imread('la web.png'))/255,'FaceColor','texturemap','Edgecolor','none');
camproj perspective;
camva(6);
xlim([-4 4]);
ylim([-4 4]);
Zlim([-7 7]);
camlight right;
lighting phong;
view([(((45-0)/(10+1))*(t+1)+0) (((45-10)/(10+1))*(t+1)+10)]);
F = getframe(fig);
aviobj = addframe(aviobj,F);
end
close(fig)
aviobj = close(aviobj);

y escupió esto:



lo mejor de todos esto, es que esta función es contínua gracias a que la Función unitaria de Heaviside es contínua si $k$ es un número real, lo que quisiera es poder resolver un problema de este tipo:

Empiesa a llover, y en un instante dado se conocen las posiciones $\boldsymbol{r} _i$ de $N$ gotas de agua, con esto, describa el patrón de interferencia despues de un tiempo $t$.

me imagino que la forma de esa función debe ser algó como así:

$\sum_{i = 1}^N \frac{1}{1+{e}^{-2\times 20\left( t-a_i \right) } }\frac{3}{1+{e}^{-2\times 20\left(\lambda\,T\left( t-a_i \right) -\vert\vert \boldsymbol{r} -\boldsymbol{r}_i\vert\vert +\frac{1}{2} \right) } } \frac{{e}^{-0.5\left\vert{}\vert\vert \boldsymbol{r} -\boldsymbol{r}_i\vert\vert-\lambda \,T\left( t-a_i \right) \right\vert{} } }{\vert\vert \boldsymbol{r} -\boldsymbol{r}_i\vert\vert+2} \times ...$
$...\times\sin \left[2\pi \left( \frac{\vert\vert \boldsymbol{r} -\boldsymbol{r}_i\vert\vert}{\lambda } -\frac{\left( t-a_i \right)}{T} \right) +\phi \right]$

donde los $a_i$ son los tiempos que tardan cada gota $i$ en llegar a la superficie del agua. Habría que calcular la distancia mas corta entre un punto (la gota $i$) y un plano (la superficie del agua) y aplicar... nosé, imagino que caída libre. Pero bueno ya vez más o menos a donde es que quiero llegar

[arivasm]

Re: Ecuacion de ondas

Caramba, natanael, esta vez me has dejado sin aliento! El resultado final es perfecto!!!

[pethr2]

Re: Ecuacion de ondas

Hola
He copiado el codigo a Matlab y no he podido reproducir como lo muestras me dice que no existe La web. png

podrias ayudarme soy nuevo en este campo
y en el pprograma

[Lorentz]

Re: Ecuacion de ondas

Si no tienes ningún archivo de imagen llamado la_web.png entonces no te va a salir la imagen de La Web de Física. De todas formas creo que si es la primera vez que utilizas un programa como Matlab, deberías realizar algún ejercicio más sencillo para ir asimilando conceptos y comandos, y posteriormente ya intentar cosas de este estilo.

Un saludo