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
void calculaVectores_valores(int N, double masa){

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