PROPAGACIÓN DE ONDAS EN UN CRISTAL UNIDIMENSIONAL

1. PROPAGACIÓN DE ONDAS EN UN CRISTAL UNIDIMENSIONAL

2. ¿EN QUÉ CONSISTE NUESTRO PROYECTO? • En el estudio de la propagación de ondas mediante una simulación por medio de muelles y masas. • Resuelto por dos métodos: • Euler menos preciso. • Verletmás preciso.

3. ¿CUÁL ES SU UTILIDAD? • Estudio de la propagación de ondas longitudinales y mecánicas en distintos medios.

4. OBJETIVOS. • Calcular la posición en función del tiempo. • Dos masas. • N masas. • Solución analítica (dos masas). • Representación gráfica de las posiciones de las masas. • Conservación de la energía.

5. EXPLICACIÓN ANALÍTICA. • Ley de Hooke: F=-kΔx • Dos masas tres muelles con k igual. F₁=-k₁Δx₁+k₂(Δx₂-Δx₁) F₂= -k₃ Δx₂- k₂ (Δx₂-Δx₁) • Se busca una solución con la siguiente forma: x1=A1cos(ωt), x2=A2cos(ωt)

6. La expresión matricial es la siguiente. • Se calculan los valores y vectores propios de la matriz. λ₁=-1 y λ₂=-3. • ω²

7. Para hallar las condiciones iniciales, t=0 α=β= x₀/2. • Matriz caso N partículas,

8. EXPLICACIÓN DEL CÓDIGO. • Crearemos una función con la que el usuario pueda modificar los parámetros iniciales. • PARA DOS PARTÍCULAS. • PARA N PARTÍCULAS.

9. DOS PARTÍCULASEULER • Definición de parámetros y condiciones para t=0. • Bucle que corre el tiempo y calcula la fuerza, velocidad y la posición en los distintos t. • Representación gráfica de la posición de las dos masas.

10. function [r1,r2]=cristal_euler(m1,m2,r1,r2,k1,k2,k3,N) • %Hemos convertido las condiciones y parametros iniciales en vectores • r1(1)=r1; • r2(1)=r2; • k(1)=k1; • k(2)=k2; • k(3)=k3; • v1=0; • v2=0; • %Tiempo inicial • t(1)=0; • %Diferencial de tiempo • dt=0.01; • %Comenzamos el bucle • for i=1:N, • t(i+1)=t(i)+dt; • %Hacemos que corran las fuerzas, velocidades y posiciones • F1=(-k(1)*r1(i))+(k(2)*(r2(i)-r1(i))); • F2=(-k(3)*r2(i))-(k(2)*(r2(i)-r1(i))); • v1(i+1)=v1(i)+(dt*F1)/m1; • v2(i+1)=v2(i)+(dt*F2)/m2; • r1(i+1)=r1(i)+dt*v1(i); • r2(i+1)=r2(i)+dt*v2(i); • end • %Creamos las graficas • plot(t,r1,'k','linewidth',5) • xlabel('Time (s)') • ylabel('Position (m)') • axis([0 100 -2 2]) • holdon • plot(t,r2,'m','linewidth',2); • hold off

11. DOS PARTÍCULAS VERLET • Definición de los parámetros y condiciones iniciales. • Cálculo de la fuerza en el primer instante. • Cálculo de la segunda posición. • Bucle que corre el tiempo y calcula las fuerzas y las posiciones en todos los t. • Representación gráfica.

12. function [r1,r2]=cristal_verlet(m1,m2,r1,r2,k1,k2,k3,N) • %Hemos convertido las condiciones y parametros iniciales en vectores • r1(1)=r1; • r2(1)=r2; • k(1)=k1; • k(2)=k2; • k(3)=k3; • v1=0; • v2=0; • %Tiempo inicial • t(2)=0; • %Diferencial de tiempo • dt=0.01; • %Condiciones iniciales para verlet, definimos fuerza y posicion • F1=(-k(1)*r1(1))+(k(2)*(r2(1)-r1(1))); • F2=(-k(3)*r2)-(k(2)*(r2(1)-r1(1))); • r1(2)=((1/2)*(F1/m1)*dt.^2)+(v1*dt)+(r1(1)); • r2(2)=((1/2)*(F2/m2)*dt.^2)+(v2*dt)+(r2(1)); • %Crearemos un bucle para que vaya corriendo el tiempo y se vayan • %modificando las fuerzas y las posiciones iniciales. • for i=2:N-1, • t(i+1)=t(i)+dt; • F1=(-k(1)*r1(i))+(k(2)*(r2(i-1)-r1(i-1))); • F2=(-k(3)*r2(i))-(k(2)*(r2(i-1)-r1(i-1))); • r1(i+1)=(2*r1(i))-r1(i-1)+(F1/m1)*dt.^2; • r2(i+1)=(2*r2(i))-r2(i-1)+(F2/m2)*dt.^2; • end • %pintamos las dos graficas • plot(t,r1) • xlabel('Time (s)') • ylabel('Position (m)') • holdon • plot(t,r2,'r') • hold off

13. RESULTADOSm,k y r iguales

14. m, k iguales, r1=-1, r2=1

15. m1=1, m2=3, r y k iguales

16. Batimiento

17. N PARTÍCULAS • Resuelto por Euler. • Diferencias con respecto a dos partículas: • Bucle para definir condiciones iniciales. • Dos bucles anidados: • 1º bucle que corre el tiempo, determina la fuerza, velocidad y posición de la primera y última partícula en todos los t. • 2º bucle que corre el nº de partículas y calcula lo mismo que el primer bucle anidado.

18. function x=G6(N,r,m,k,T) %usage G6(3,[1 1 1],[1 1 1],[1 1 1 1],20000) for i=1:N, x(i,1)=r(i); v(i,1)=0; end t(1)=0; %definimos el tiempo incial dt=0.001; %definimos el diferencial de tiempo Ec=zeros(N,T); Ep=zeros(N,T); for i=1:T, %tiempo t(i+1)=t(i)+dt; F(1,i)=-(k(1)+k(2))*x(1,i)+k(2)*x(2,i); F(N,i)=k(N)*x(N-1,i)-(k(N)+k(N+1))*x(N,i); v(1,i+1)=v(1,i)+(dt*F(1,i)/m(1)); v(N,i+1)=v(N,i)+(dt*F(N,i)/m(N)); x(1,i+1)=x(1,i)+dt*v(1,i); x(N,i+1)=x(N,i)+dt*v(N,i); for j=2:N-1, %particulas F(j,i)=k(j)*x(j-1,i)-(k(j)+k(j+1))*x(j,i)+k(j+1)*x(j+1,i); v(j,i+1)=v(j,i)+(dt*F(j,i)/m(j)); x(j,i+1)=x(j,i)+dt*v(j,i); Ec(j,i)=0.5*m(j)*v(j,i).^2; Ep(j,i)=0.5*k(j)*x(j,i).^2+0.5*k(j+2)*x(j+1,i).^2+0.5*k(j+1)*(x(j+1,i)-x(j,i)).^2; end Ec(1,i)=0.5*m(1)*v(1,i).^2; Ec(N,i)=0.5*m(N)*v(N,i).^2; Ep(1,i)=0.5*k(1)*x(1,i).^2+0.5*k(3)*x(2,i).^2+0.5*k(2)*(x(2,i)-x(1,i)).^2; Ep(N,i)=0.5*k(N+1)*x(N,i).^2+0.5*k(N-1)*x(N-1,i).^2+0.5*k(N)*(x(N,i)-x(N-1,i)).^2; end ET=Ec+Ep; Et=sum(ET); figure (1) plot(Et) Title('Energy') for p=1:N, figure (2) subplot(1,N,p), plot(x(p,:)); title('Masa') xlabel('Time') ylabel('Position') end

19. RESULTADOS