Oscilaciones

Sistema formado por tres partículas N=3

Sistema formado por N partículas

Actividades

Relación de dispersión

Referencias

Código fuente

En esta página, vamos a calcular los modos normales de vibración de una cadena diatómica lineal formada por moléculas de masas m  y M unidas por muelles elásticos de constante, k tal como se muestra en la figura

Sistema formado por dos partículas N=2

Ecuaciones del movimiento 

Buscamos una solución de la forma

x₁=A₁cos(ωt), x₂=A₂cos(ωt)

El sistema de ecuaciones se expresa en forma matricial

Como vemos la matriz ha dejado de ser simétrica

Las frecuencias son los valores propios de la matriz cuadrada  |M- ω²I|=0, las raíces de la ecuación de segundo grado en ω² denominado polinomio característico

La amplitudes para cada uno de los modos de vibración se calculan resolviendo el sistema de ecuaciones homogéneo

• Para el modo normal ω₁, se calcula A₂₁ en función de A₁₁

• Para el modo normal ω₂, se calcula A₂₂ en función de A₁₂

El movimiento general de las partículas es una combinación lineal de los dos modos normales de vibración

x₁=A₁₁cos(ω₁t)+A₁₂cos(ω₂t)
x₂=A₂₁cos(ω₁t)+A₂₂cos(ω₂t)

Los coeficientes A₁₁ y A₁₂ se determinan a partir de las condiciones iniciales: en el instante t=0, la posición inicial de las partículas es x₁₀ y x₂₀

Sistema formado por tres partículas N=3

Ecuaciones del movimiento 

Buscamos una solución de la forma

x₁=A₁cos(ωt), x₂=A₂cos(ωt), x₃=A₃cos(ωt)

El sistema de ecuaciones se expresa en forma matricial

Como vemos la matriz ha dejado de ser simétrica

Las frecuencias son los valores propios de la matriz cuadrada, las raíces de la ecuación cúbica en ω²denominado polinomio característico

La amplitudes para cada uno de los modos de vibración se calculan resolviendo el sistema de ecuaciones homogéneo

• Para el modo normal ω₁, se calcula A₂₁ y A₃₁ en función de A₁₁

• Para el modo normal ω₂, se calcula A₂₂ y A₃₂ en función de A₁₂

• Para el modo normal ω₃, se calcula A₂₃ y A₃₃ en función de A₁₃

El movimiento general de las partículas es una combinación lineal de los tres modos normales de vibración

x₁=A₁₁cos(ω₁t)+A₁₂cos(ω₂t)+A₁₃cos(ω₃t)
x₂=A₂₁cos(ω₁t)+A₂₂cos(ω₂t)+A₂₃cos(ω₃t)
x₃=A₃₁cos(ω₁t)+A₃₂cos(ω₂t)+A₃₃cos(ω₃t)

Los coeficientes A₁₁, A₁₂ y A₁₃ se determinan a partir de las condiciones iniciales: en el instante t=0, la posición inicial de las partículas es x₁₀, x₂₀ y x₃₀

Sistema formado por N partículas

Para N partículas calculamos los valores y los vectores propios de la matriz

Como la matriz no es simétrica no podemos aplicar el procedimiento de Jacobi. Como alternativa, obtenemos los coeficientes del polinomio característico aplicando el procedimiento de Leverrier y calculamos las raíces del polinomio aplicando el procedimiento del punto medio para raíces múltiples.

Calculamos las amplitudes A_i mediante la relación de recurrencia

Donde m₁=m, m₂=M, m₃=m, m₄=M,…

• Para el modo normal ω₁ se calculan las amplitudes A₂, A₃…. A_N en función de A₁, que denominaremos A₁₁, A₂₁, A₃₁… A_N1

• Para el modo normal ω_j se calculan las amplitudes A₂, A₃…. A_N en función de A_1, que denominaremos A_1j, A_2j, A_3j… A_Nj

• Para el modo normal ω_N se calculan las amplitudes A₂, A₃…. A_N en función de A_1, que denominaremos A_1N, A_2N, A_3N… A_NN.

Una vez calculadas las amplitudes A_ij para cada frecuencia ω_j se normalizan

Actividades

Se ha fijado

• La masa, m=1 kg

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

Se introduce

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

• La masa M de las otras partículas en el control de edición titulado Masa

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 ω₁.

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

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

Las frecuencias de los modos de vibración se guardan en el control área de texto situado a la izquierda del applet.

Se pulsa el botón titulado Gráfica, para representar las frecuencias de los modos de vibración (en el eje vertical) en función de su orden.

Relación de dispersión

Consideremos una cadena lineal de moléculas diatómicas separados una distancia a en la situación de equilibrio, tal como se muestra en la figura.

Consideremos el movimiento de dos átomos contiguos.

• El desplazamiento del átomo i de masa M denominamos y_i

• El desplazamiento del átomo i+1 de masa m denominamos x_i+1

Ecuaciones del movimiento

Buscamos una solución de la forma

y_i=A_icos(ωt), con A_i=Asen(Kia)

x_i+1=A_i+1cos(ωt), con A_i+1=Bsen(K(i+1)a)

donde K se denomina número de onda

Introducimos estas soluciones en las dos ecuaciones diferenciales y obtenemos un sistema homogéneo de dos ecuaciones con dos incógnitas A y B.

A(-ω²M+2k)sen(Kia)-Bk(sen(K(i+1)a)+sen(K(i-1)a))=0
-Ak(sen(K(i+2)a)+sen(Kia))+B(-ω²m+2k) sen(K(i+1)a)=0

Teniendo en cuenta que

sen(K(i+1)a)+sen(K(i-1)a)=2sen(Kia)cos(Ka)
sen(K(i+2)a)+sen(Kia)= 2sen(K (i+1)a)cos(Ka)

y eliminamos A y B en el sistema homogéneo de dos ecuaciones

(-ω²M+2k) (-ω²m+2k)sen(Kia)sen(K(i+1)a)-4k²sen(Kia)cos²(Ka)sen(K(i+1)a))=0

Simplificando el factor común sen(Kia)sen(K(i+1)a)

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

En la figura, se representa se representa la frecuencia angular ω en función del número de onda K en el intervalo (-π/a, +π/a). La curva superior (signo +) se denomina rama óptica, y la inferior (signo -) se denomina rama acústica.

public class Matriz {
	public int n; //dimensión
	private double[][] x;

public Matriz(int n) {
	this.n=n;
	x=new double[n][n];
	for(int i=0; i<n; i++){
		for(int j=0; j<n; j++){
			x[i][j]=0.0;
		}
	}
}
public Matriz(double[][] x) {
	this.x=x;
	n=x.length;
}

double traza(){
	double tr=0.0;
	for(int i=0; i<n; i++){
		tr+=x[i][i];
	}
	return tr;
}
//suma de dos matrices
static Matriz suma(Matriz a, Matriz b){
	Matriz resultado=new Matriz(a.n);
	for(int i=0; i<a.n; i++){
		for(int j=0; j<a.n; j++){
			resultado.x[i][j]=a.x[i][j]+b.x[i][j];
		}
	}
	return resultado;
}
//producto de dos matrices
static Matriz producto(Matriz a, Matriz b){
	Matriz resultado=new Matriz(a.n);
	for(int i=0; i<a.n; i++){
		for(int j=0; j<a.n; j++){
			for(int k=0; k<a.n; k++){
				resultado.x[i][j]+=a.x[i][k]*b.x[k][j];
			}
		}
	}
	return resultado;
}
//polinomio característico

public double[] polCaracteristico(){
	Matriz pot=new Matriz(n);
//matriz unidad
	for(int i=0; i<n; i++){
		pot.x[i][i]=1.0;
	}
	double[] p=new double[n+1];
	double[] s=new double[n+1];
	for(int i=1; i<=n; i++){
		pot=Matriz.producto(pot, this);
		s[i]=pot.traza();
	}
	p[0]=1.0;
	p[1]=-s[1];
	for(int i=2; i<=n; i++){
		p[i]=-s[i]/i;
		for(int j=1; j<i; j++){
			p[i]-=s[i-j]*p[j]/i;
		}
	}
	return p;
}

public String toString(){
	String texto="\n";
	for(int i=0; i<n; i++){
		for(int j=0; j<n; j++){
			texto+="\t "+(double)Math.round(1000*x[i][j])/1000;
		}
		texto+="\n";
	}
	texto+="\n";
	return texto;
}

}

public abstract class Ecuacion {
	protected final double CERO=1e-10;
	protected final double ERROR=0.0001;
	protected final int MAXITER=100;
	protected final int MAXRAICES=20;
	protected double raices[]=new double[MAXRAICES];
	protected int iRaiz=0;

protected double puntoMedio(double a, double b)throws RaizExcepcion{
	double m, ym;
	int iter=0;
	do{
		m=(a+b)/2;
		ym=f(m);
		if(Math.abs(ym)<CERO) break;
		if(Math.abs((a-b)/m)<ERROR) break;

		if((f(a)*ym)<0) b=m;
		else a=m;
		iter++;
	}while(iter<MAXITER);
	if(iter==MAXITER){
		throw new RaizExcepcion("No se ha encontrado una raíz");
	}
	return m;
}

protected void explorar(double xIni, double xFin, double dx){
	double y1, y2;
	iRaiz=0;
	y1=f(xIni);
	for(double x=xIni; x<xFin; x+=dx){
		y2=f(x+dx);
//Uno de los extremos del intervalo es raíz
		if(Math.abs(y1)<CERO && iRaiz<MAXRAICES){
			raices[iRaiz++]=x;
			y1=y2;
			continue;
		}
//no hay raíz en este intervalo
		if(y1*y2>=0.0){
			y1=y2;
			continue;
		}
//hay una raíz en este intervalo
		if(iRaiz<MAXRAICES){
			try{
				raices[iRaiz]=puntoMedio(x, x+dx);
				iRaiz++;
			}catch(RaizExcepcion ex){
				System.out.println(ex.getMessage());
			}
		}
		y1=y2;
	}
}
public double[] hallarRaices(double ini, double fin, double paso){
	explorar(ini, fin, paso);
	double solucion[]=new double[iRaiz];
	for(int i=0; i<iRaiz; i++){
		solucion[i]=(double)Math.round(raices[i]*1000)/1000;
	}
	return solucion;
}

abstract public double f(double x);
}

class RaizExcepcion extends Exception {
	public RaizExcepcion(String s) {
	super(s);
}
}

public class Funcion1 extends Ecuacion{
	double[]coef;

Funcion1 (double[]coef){
	this.coef=coef;
}
public double f(double x){
	double y=0.0;
//sucesivas potencias de x, se puede utilizar tambien la funcion Math.pow
	double[] pot_x=new double[coef.length];
	pot_x[0]=1.0;
	for(int i=1; i<coef.length; i++){
		pot_x[i]=pot_x[i-1]*x;
	}
//valores de los sucesivos términos
	for(int i=0; i<coef.length; i++){
		y+=coef[i]*pot_x[coef.length-i-1];
	}
	return y;
}

}

//calcula los valores y vectores propios
//matriz M
	double[] k=new double[N+1];
	double[] m=new double [N+1];
	for(int i=0; i<k.length; i++){
		k[i]=1.0;
		if(i%2==0) m[i]=masa;
		else 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];
//polinomio característico
	Matriz p=new Matriz(matrix);
	double pol[]=p.polCaracteristico();
	System.out.println("Polinomio característico");
	for(int i=0; i<pol.length; i++){
		System.out.print((double)Math.round(pol[i]*1000)/1000+" , ");
	}
//valores propios
	double[] valores=new Funcion1(pol).hallarRaices(0.0, 10, 0.05);
	double[][] vectores=new double[N][N];
 	System.out.print("Raíces de una ecuación");
	for(int i=0; i<N; i++){
		System.out.print(" "+valores[i]);
	}
	System.out.println(""); 
//vectores propios
	double[] aux=new double[N+1];
	double temp;
	for(int j=0; j<N; j++){
		aux[0]=0.0;
		aux[1]=1.0;
		for(int i=1; i<N; i++){
			aux[i+1]=-k[i-1]*aux[i-1]/k[i]+(k[i-1]/k[i]+1-valores[j]*m[i])*aux[i];
		}
		temp=0.0;
		for(int i=1; i<=N; i++){
			temp+=aux[i]*aux[i];
		}
		for(int i=0; i<N; i++){
			vectores[i][j]=aux[i+1]/Math.sqrt(temp);
		}
	}