Oscilaciones |
Osciladores acoplados |
Ecuaciones del movimiento | |||||||||||||||||||||||||||||
| 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 movimientoEn 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
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) 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
Hasta aquí hemos formulado el problema general, pasamos a hora a resolver diversos casos particulares
Casos particularesUn 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 M por el procedimiento de Jacobi (véase el código fuente)
Ejemplos:
ActividadesSe ha fijado
Se introduce
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ónComo 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.
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-mω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).
|
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
Cadena monoatómica
