Katea diatomiko lineal baten oszilazioak

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

Oszilazioak

 
Osziladore akoplatuak
marca.gif (847 bytes)Katea diatomikoa
 Bi partikulaz osatutako sistema, N=2

Hiru partikulaz osatutako sistema, N=3

N partikulaz osatutako sistema

Saiakuntza

Dispertsio erlazioa

Erreferentziak

Programazio-kodea

 

Orri honetan katea diatomiko lineal baten oszilazioak aztertuko ditugu. Demagun bi atomo mota ditugula, m eta M masadunak, k konstantedun malguki batez lotuta daudela, eta bikotea behin eta berriz errepikatzen dela.

 

Bi partikulaz osatutako sistema, N=2

Higiduraren ekuazioak:

Bila ditzagun itxura honetako soluzioak:

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

Matrize gisa, honela adierazten da ekuazio-sistema hori:

Honako kasuan, matrizea ez da simetrikoa.

Bibrazioen modu normalak kalkulatzen dira ondoko matrizearen determinantea zerora berdinduz: |M- ω2I|=0, eta hemen I matrize unitatea da. ω2-ren menpekoa den bigarren graduko polinomioaren soluzioak (polinomio karakteristikoa deiturikoa) M matrizearen balio propioak dira. Bektore propioak berriz, maiztasun horietako bakoitzarentzat, Ai koefizienteen balioak dira.

ωj  maiztasun normal bakoitzarentzat, modu normalen Ai anplitudeak honako ekuazio-sistema homogeneoa ebatziz kalkula daitezke:

  • ω1  modu normalarentzat, A21 kalkulatzen da A11-en menpe.

  • ω2  modu normalarentzat, A22 kalkulatzen da A12-ren menpe.

Partikula bien higidura orokorra bi oszilazio-modu normalen konbinazioa izango da.

x1=A11cos(ω1t)+A12cos(ω2t)
x2=A21cos(ω1t)+A22cos(ω2t)

A11 eta A12 koefizienteak hasierako baldintzetatik determinatzen dira, alegia, t=0 aldiunean bi partikulen posizioak x10 eta x20 dira.

Hiru partikulaz osatutako sistema, N=3

Higiduraren ekuazioak:

Bila ditzagun itxura honetako soluzioak:

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

Matrize gisa, honela adierazten da ekuazio-sistema hori:

Honako kasuan, matrizea ez da simetrikoa.

Bibrazioen modu normalak lortzen dira matrize karratuaren balio propioak kalkulatzen (polinomio karakteristikoa deiturikoa). Kasu honetan, ω2-ren menpeko hirugarren graduko ekuazioa (kubikoa) da.

Maiztasun horietako bakoitzarentzat, Ai koefizienteak kalkulatzen dira honako ekuazio-sistema homogeneoa ebatziz:

  • ω1  modu normalarentzat, A21 eta A31 kalkulatzen dira A11-en menpe.

  • ω2  modu normalarentzat, A22 eta A32 kalkulatzen dira A12-ren menpe.

  • ω3  modu normalarentzat, A23 eta A33 kalkulatzen dira A13-ren menpe.

Partikulen higidura orokorra, ez da modu normal bat, hiru bibrazio-moduen konbinazio lineala baizik.

x1=A11cos(ω1t)+A12cos(ω2t)+A13cos(ω3t)
x2=A21cos(ω1t)+A22cos(ω2t)+A23cos(ω3t)
x3=A31cos(ω1t)+A32cos(ω2t)+A33cos(ω3t)

A11, A12 eta A13 koefizienteak hasierako baldintzetatik determinatzen dira, alegia, t=0 aldiunean, hiru partikulen posizioak x10 , x20 eta x30 dira.

N partikulaz osatutako sistema

N partikulekin honako matrizearen balio propioak eta bektore propioak kalkulatu behar dira:

Matrizea ez da simetrikoa, eta horregatik ezin da Jacobi-ren prozedura erabili, baina polinomio karakteristikoaren koefizienteak kalkula ditzakegu Leverrier-en prozedura aplikatuz eta, ondoren, polinomioaren erro anitzak kalkula ditzakegu prozedura numeriko batez, besteak beste, erdiko puntuaren prozedura.

Ai anplitudeak kalkulatzeko erabil dezagun honako errepikapen erlazioa edo ekuazio-sistema homogeneoa:

Hemen m1=m , m2=M , m3=m , m4=M,…

  • ω1 modu normalarentzat, A2 , A3…, AN anplitudeak kalkulatzen dira A1-en menpe, baina honela izendatuko ditugu: A11, A21, A31… AN1.

  • ωj modu normalarentzat, A2 , A3…, AN anplitudeak kalkulatzen dira A1-en menpe, baina honela izendatuko ditugu: A1j, A2j, A3j… ANj.

  • ωN  modu normalarentzat, A2 , A3…, AN anplitudeak kalkulatzen dira A1-en menpe, baina honela izendatuko ditugu: A1N, A2N, A3N… ANN.

Anplitude guztiak kalkulatu ondoren (Aij) normalizatu egin behar dira: (N2 anplitude)

 

Saiakuntza

Finkotzat hartzen dira:

  • Partikula "txikien" masa, m=1 kg.

  • Malgukien konstante elastikoa, k=1 N/m.

Eta aukeran idatz daitezke:

  • Partikula kopurua, N, dagokion kontrolean idatziz.

  • Partikula "handien" Masa, M, dagokion kontrolean idatziz.

Berria botoia sakatu.

Lehenik, partikula-multzoa bibratzen ikusten da, lehen bibrazio-modu normalean. Goiko eta ezkerreko erpinean maiztasun angeluarra idatziz erakusten da: ω1.

Hurrengoa>> izeneko botoia sakatuz, hurrengo bibrazio-modu normala ikusten da, ω2, eta behin eta berriz sakatuz, ω3 ,ω4 , … etab.

Aurrekoa<< izeneko botoia sakatuz, aurreko bibrazio-modu normala ikusten da.

Modu normalen kopurua eta partikula-kopurua beti dira berdinak. Leihatilaren ezkerraldeko testu-zutabean, programak bibrazio-moduen maiztasunak idazten ditu.

Azkenean, Grafikoa botoia sakatuz, bibrazio-moduen maiztasunak adierazten dira ordenaren arabera. 
 

 

Dispertsio erlazioa

Ondorengo irudiak erakusten du N partikulaz osaturiko katea lineal bat, baina bakoitiak m masa dute eta bikoitiek M. Demagun oreka-posizioan partikulen arteko separazioa a dela.

Azter dezagun bi atomo moten higidura:

  • Dei diezaiogun M masadun i-garren atomoaren desplazamenduari, yi

  • Eta dei diezaiogun m masadun (i+1)-garren atomoaren desplazamenduari, xi+1

Hona hemen higiduraren ekuazioak:

Bila ditzagun itxura honetako soluzioak:

yi=Aicos(ωt), eta Ai=A·sin(Kia)

xi+1=Ai+1cos(ωt), eta Ai+1=B·sin(K(i+1)a)

Hemen, K-ri uhin-zenbaki deritzo.

Ordezka ditzagun soluzio horiek ekuazio diferentzialetan eta lortzen da, bi ekuaziodun eta bi ezezaguneko sistema homogeneoa. Ezezagunak A eta B dira:

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

Eta kontutan hartzen bada:

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

Eta eliminatzen baditugu A eta B bi ekuaziodun sistema homogeneoan:

(-ω2M+2k) (-ω2m+2k)sin(Kia)sin(K(i+1)a)-4k2sin(Kia)cos2(Ka)sin(K(i+1)a))=0

Faktore komun bat osorik sinplifikatzen da: sin(Kia)sin(K(i+1)a)

Orduan:

Ekuazio horrek erlazionatzen ditu w maiztasun angeluarra eta K uhin-zenbakia, alegia, dispertsio-erlazioa deritzona.

Ondoko irudiak erakusten du ω maiztasun angeluarra K uhin-zenbakiaren menpe, eta (-π/a, +π/a) tartean. Goiko kurba urdinari (positiboari) adar optiko deritzo eta azpiko gorriari (negatiboari) adar akustiko.

 

 

Erreferentzia

Runk R. B. Stul J. L. Anderson G. L. A laboratory analog for lattice dynamics. Am. J. Phys. (31) 1963, pp. 915-921

 

Programazio-kodea

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);
		}
	}