Oscilaciones  de una cadena monoatómica lineal

prev.gif (1231 bytes)home.gif (1232 bytes)next.gif (1211 bytes)

Oscilaciones

Osciladores acoplados
Dos acoplados
Tres acoplados
marca.gif (847 bytes)Cadena monoatómica
Varilla que pende 
 de dos muelles
Péndulo de Wilberforce
El péndulo doble
Péndulo-muelle
De las oscilaciones
a las ondas
Péndulos no acoplados
de distinta longitud
Ecuaciones del movimiento

Modos normales de vibración

Casos particulares

Actividades

Relación de dispersión

Código fuente

 

Vamos a estudiar los modos normales de vibración de un sistema formado por muelles y partículas, como continuación y generalización del sistema formado por tres osciladores acoplados. Este ejemplo nos ayudará a comprender los modos normales de vibración de una cuerda fija por sus extremos, también denominados ondas estacionarias.

Supongamos un conjunto de partículas de masas m1, m2, m3…mN-1, mN unidas por muelles elásticos de constantes k0, k1, k2, k3, … kN-1, kN.

 

Ecuaciones del movimiento

En las figuras se muestran las fuerzas sobre la primera partícula de masa m1, la partícula i y la última partícula N de masa mN

Resumiendo, las ecuaciones del movimiento de las partículas son

Buscamos una solución de la forma

xi=Aicos(ωt)

Al ser la fase inicial π/2, las partículas parten del reposo en el instante t=0, de la posición x0i

Teniendo en cuenta que la derivada segunda de xi respecto del tiempo t es

Obtenemos el siguiente sistema de ecuaciones homogéneo

Que podemos expresar en forma matricial de la forma

De forma simbólica podemos expresar M·X= ω2X,  donde X es el vector columna de los coeficientes A1, A2... AN.

Las frecuencias de los modos normales de vibración se calculan haciendo que el determinante de la matriz |M- ω2I|=0, donde I es la matriz unidad. Las raíces del polinomio de grado N en ω2 son los valores propios de la matriz M. Los vectores propios son los valores de los coeficientes Ai para cada una de dichas frecuencias.

También, podemos calcular las amplitudes Ai de los modos normales de vibración para cada una de las frecuencias ωj, empleando la relación de recurrencia

  • Para el modo normal de frecuencia ω1 se calculan las amplitudes A2, A3…. AN en función de A1, que denominaremos A11, A21, A31… AN1, que son a su vez, las componentes del vector propio correspondiente al valor propio ω21

  • Para el modo normal de frecuencia ωj se calculan las amplitudes A2, A3…. AN en función de A1 que denominaremos A1j, A2j, A3j… ANj, que son a su vez, las componentes del vector propio correspondiente al valor propio ω2j

  • Para el modo normal de frecuencia ωN se calculan las amplitudes A2, A3…. AN en función de A1 que denominaremos A1N, A2N, A3N… ANN., que son a su vez, las componentes del vector propio correspondiente al valor propio ω2N

Movimiento de las partículas

El movimiento de cada partícula es una superposición de todos los modos de vibración

x1=A11cos(ω1·t)+A12cos(ω2·t)+...A1N·cos(ωN·t)
x2
=A21cos(ω1·t)+A22cos(ω2·t)+....A2N·cos(ωN·t)
...........
xN
=AN1cos(ω1·t)+AN2cos(ω2·t)+....ANN·cos(ωN·t)

Se determinan los coeficientes A11,A12.... A1N  a partir de las condiciones iniciales: en el instante t=0, las posiciones de las partículas son x10, x20, x30,...xN0.

 

Modos normales de vibración

  • El primer modo normal de vibración de frecuencia ω1 se establece cuando A12=A13=…A1N=0 y A11≠0

El movimiento de las partículas es

x1=A11cos(ω1·t)
x2
=A21cos(ω1·t)
...........

xN
=AN1cos(ω1·t)

Donde A11, A21, ...AN1 son las componentes del vector propio correspondientes correspondiente al valor propio ω21

  • El modo normal de vibración de frecuencia ωj se establece cuando A11=A13=…A1N=0 y A1j≠0

El movimiento de las partículas es

x1=A1jcos(ωj·t)
x2
=A2jcos(ωj·t)
...........

xN
=ANjcos(ωj·t)

Donde A1j, A2j, ...ANj son las componentes del vector propio correspondientes correspondiente al valor propio ω2j

Hasta aquí hemos formulado el problema general, pasamos a hora a resolver diversos casos particulares

 

Casos particulares

Un caso particular importante, es una cadena monoatómica lineal cuando las constantes de los muelles elásticos tienen el mismo valor k0=k1=k2=k3, …=kN-1=kN=k,  y las partículas tienen la misma masa m0=m1=m2=m3, …=mN-1=mN=m

Calculamos los valores propios ω2 y los vectores propios de la matriz simétrica matriz simétrica M por el procedimiento de Jacobi (véase el código fuente)

Ejemplos:

Ecuaciones del movimiento

Buscamos una solución de la forma

x1=A1cos(ωt), x2=A2cos(ωt)

En forma matricial, el sistema de ecuaciones se expresa

Las frecuencias ω2 son los valores propios de la matriz cuadrada y los vectores propios, las amplitudes de los modos normales de vibración

Valores propios ω2

Vectores propios

k/m

(A11, A11)

3k/m

(A12, -A12)

Los coeficientes A11 y A12 se determinan a partir de las condiciones iniciales

Las ecuaciones del movimiento son:

Buscamos una solución de la forma

x1=A1cos(ωt), x2=A2cos(ωt), x3=A3cos(ωt)

En forma matricial, el sistema de ecuaciones se expresa

Las frecuencias ω2 son los valores propios de la matriz cuadrada y los vectores propios, las amplitudes de los modos normales de vibración

Valores propios ω2

Vectores propios

(A12, 0, -A12)

Los coeficientes A11, A12 y A13 se determinan a partir de las condiciones iniciales

  • Sistema formado por cuatro partículas N=4

Buscamos una solución de la forma

x1=A1cos(ωt), x2=A2cos(ωt), x3=A3cos(ωt), x4=A4cos(ωt)

En forma matricial, el sistema de ecuaciones se expresa

Las frecuencias ω2 son los valores propios de la matriz cuadrada y los vectores propios, las amplitudes de los modos normales de vibración

Valores propios ω2

Vectores propios

Los coeficientes A11, A12, A13 y A14 se determinan a partir de las condiciones iniciales

 

Actividades

Se ha fijado

  • La masa de las partículas, m=1 kg

  • La constante elástica de los muelles, k=1N/m

Se introduce

  • El número N de partículas, en el control de selección titulado Número de partículas

Se pulsa el botón titulado Nuevo

Observamos el primer modo normal de vibración, en la parte superior izquierda se proporciona la frecuencia angular ω1.

Se pulsa el botón titulado Siguiente>>,  observamos los siguientes modos normales ω2, … etc.

Se pulsa el botón titulado Anterior<< para observar el modo normal de vibración anterior.

 

                  

 

Relación de dispersión

Como vemos en la figura tenemos N partículas de masa m unidas a N+1 muelles iguales de constante k, cuyos extremos están fijos. La separación de equilibrio entre las partículas es a.

Supongamos que en un instante dado t, la partícula 1 se desplaza x1, la partícula 2 se desplaza x2, ... la partícula i se desplaza xi, etc.

La ecuación del movimiento para la partícula i será entonces

Supongamos, que el sistema vibra en un modo de frecuencia w. Cada partícula describirá un M.A.S. de la misma frecuencia w , pero cuya amplitud Ai vamos a calcular.

xi=Ai·cos(w t)

Introduciendo esta expresión en la ecuación diferencial que describe el movimiento de cada partícula, obtenemos, la relación entre las amplitudes de los M.A.S. de las partículas i+1, i, e i-1.

Buscamos una solución a esta ecuación de la forma

Ai=A·sen(K·ia)

donde K es el número de onda K=2p/l

Asen(Kia+Ka)+ Asen(Kia-Ka)= Asen(Kia)(2-2/k)

Esta ecuación que relaciona la frecuencia angular w con el número de onda K, se denomina relación de dispersión.

En la figura, se representa la frecuencia angular ω en función del número de onda K en el intervalo (-π/a, +π/a).

 

Código fuente

public class Jacobi {
//el vector d son los valores propios
//Las columnas de la matriz v son los vectores propios
//devuelve el número de iteracciones
static int calcula(double[][] a, int n, double[] d, double[][] v) {
//Computes all eigenvalues and eigenvectors of a real symmetric matrix a[1..n][1..n]. On
//output, elements of a above the diagonal are destroyed. d[1..n] returns the eigenvalues of a.
//v[1..n][1..n] is a matrix whose columns contain, on output, the normalized eigenvectors of
//a. nrot returns the number of Jacobi rotations that were required.

	int j,iq,ip,i;
	double tresh,theta,tau,t,sm,s,h,g,c;
	double[] b=new double[n]; //b=vector(1,n);
	double[] z=new double[n]; //z=vector(1,n);
	for (ip=0;ip<n;ip++) { //Initialize to the identity matrix.
		for (iq=0;iq<n;iq++) v[ip][iq]=0.0;
		v[ip][ip]=1.0;
	}
	for (ip=0;ip<n;ip++) { //Initialize b and d to the diagonal of a.
		b[ip]=d[ip]=a[ip][ip];
		z[ip]=0.0; //This vector will accumulate terms of the form tapq as in equation (11.1.14).
	}
	int nrot=0;
	for (i=1;i<=50;i++) {
		sm=0.0;
		for (ip=0;ip<=n-2;ip++) { //Sum o -diagonal elements.
			for (iq=ip+1;iq<n;iq++)
				sm += Math.abs(a[ip][iq]);
		}
		if (sm == 0.0) { 
//The normal return, which relies on quadratic convergence to machine underflow.
			return nrot;
		}
		if (i < 4)
			tresh=0.2*sm/(n*n);// ...on the rst three sweeps.
		else
			tresh=0.0; //...thereafter.
		for (ip=0;ip<n-1;ip++) {
			for (iq=ip+1;iq<n;iq++) {
				g=100.0*Math.abs(a[ip][iq]);
//After four sweeps, skip the rotation if the o -diagonal element is small.
				if (i > 4 && (float)(Math.abs(d[ip])+g) == (float)Math.abs(d[ip]) 
					&& (float)(Math.abs(d[iq])+g) == (float)Math.abs(d[iq]))
					a[ip][iq]=0.0;
				else if (Math.abs(a[ip][iq]) > tresh) {
					h=d[iq]-d[ip];
					if ((float)(Math.abs(h)+g) == (float)Math.abs(h))
						t=(a[ip][iq])/h; // t = 1=(2 )
					else {
						theta=0.5*h/(a[ip][iq]); //Equation (11.1.10).
						t=1.0/(Math.abs(theta)+Math.sqrt(1.0+theta*theta));
						if (theta < 0.0) t = -t;
					}
					c=1.0/Math.sqrt(1+t*t);
					s=t*c;
					tau=s/(1.0+c);
					h=t*a[ip][iq];
					z[ip] -= h;
					z[iq] += h;
					d[ip] -= h;
					d[iq] += h;
					a[ip][iq]=0.0;
					for (j=0;j<=ip-1;j++) { //Case of rotations 1 j < p.
						g=a[j][ip];
						h=a[j][iq];
						a[j][ip]=g-s*(h+g*tau);
						a[j][iq]=h+s*(g-h*tau);
					}
					for (j=ip+1;j<=iq-1;j++) { //Case of rotations p < j < q.
						g=a[ip][j];
						h=a[j][iq];
						a[ip][j]=g-s*(h+g*tau);
						a[j][iq]=h+s*(g-h*tau);
					}
					for (j=iq+1;j<n;j++) { //Case of rotations q < j n.
						g=a[ip][j];
						h=a[iq][j];
						a[ip][j]=g-s*(h+g*tau);
						a[iq][j]=h+s*(g-h*tau);
					}
					for (j=0;j<n;j++) {
						g=v[j][ip];
						h=v[j][iq];
						v[j][ip]=g-s*(h+g*tau);
						v[j][iq]=h+s*(g-h*tau);
					}
					++nrot;
				}
			}
		}
		for (ip=0;ip<n;ip++) {
			b[ip] += z[ip];
			d[ip]=b[ip]; //Update d with the sum of tapq,
			z[ip]=0.0; //and reinitialize z.
		}
	}
	return nrot;
}
}
//Valores y vectores propios de la matriz  M
	double N=3;	
	double[] k=new double[N+1];
	double[] m=new double [N+1];
	for(int i=0; i<k.length; i++){
		k[i]=1.0;
		m[i]=1.0;
	}
	double[][] matrix=new double[N][N];
	for(int i=0; i<N; i++)
		for(int j=0; j<N; j++)
			matrix[i][j]=0.0;

	matrix[0][0]=(k[0]+k[1])/m[1];
	matrix[0][1]=-k[1]/m[1];
	for(int i=1; i<N-1; i++){
		matrix[i][i-1]=-k[i]/m[i+1];
		matrix[i][i]=(k[i]+k[i+1])/m[i+1];
		matrix[i][i+1]=-k[i+1]/m[i+1];
	}
	matrix[N-1][N-1]=(k[N-1]+k[N])/m[N];
	matrix[N-1][N-2]=-k[N-1]/m[N];
	double[] valores=new double[N];
	double[] vectores=new double[N][N];
	Jacobi.calcula(matrix, N, valores, vectores);
	ordenar(valores, vectores); 
	System.out.println("valores y vectores propios");
	for(int i=0; i<N; i++){
		System.out.print(valores[i]+"\t");
	}
	System.out.println();
	for(int i=0; i<N; i++){
		System.out.println();
		for(int j=0; j<N; j++)
		System.out.print(vectores[i][j]+"\t");
	} 
private void ordenar(double[] x, double[][] mat){
	double aux;
	double[] vAux=new double[x.length];
	for(int i=0; i<x.length-1; i++){
		for(int j=i+1; j<x.length; j++){
			if(x[i]>x[j]){
				aux=x[j];
				for(int k=0; k<x.length; k++)
					vAux[k]=mat[k][j];
				x[j]=x[i];
				for(int k=0; k<x.length; k++)
					mat[k][j]=mat[k][i];
				x[i]=aux;
				for(int k=0; k<x.length; k++)
					mat[k][i]=vAux[k];
			}
		}
	}
}

Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P. Numerical Recipes in C, Second edition,  Eigensystems. Jacobi Transformations of a Symmetric Matrix, Chapter 11º. pp. 463-469. Cambridge University Press. Código en C adaptado por el autor al lenguaje Java