Procedimientos numéricos

Valores y vectores propios (Método de Jacobi)

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) {
	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;
	for (ip=0;ip<n;ip++) { //Initialize b and d to the diagonal of a.
//This vector will accumulate terms of the form tapq as in equation } 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 void calculaValores_vectores(int N){ 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++){
	for(int i=0; i<N; i++){
		for(int j=0; j<N; j++)
}      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++){
				for(int k=0; k<x.length; k++)
				for(int k=0; k<x.length; k++)
				for(int k=0; k<x.length; 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