|
Las matrices cuadradas |
El método de Jacobi es bastante más complicado pero tiene la ventaja de hallar los valores propios y los vectores propios de una matriz simétrica A. El método de Jacobi se basa en que existen matrices ortogonales P, tales que transforman a la matriz A en una matriz cuya diagonal principal está formada por los valores propios.
Estas matrices P tienen la siguiente forma
Son nulos todos los elementos, excepto la diagonal principal y los elementos (k, l) y
su simétrico (l, k). Dado que P es ortogonal la matriz B tal que 
tiene los mismos valores propios que A, y además eligiendo
convenientemente el ángulo q se puede conseguir que los
elementos bkl y blk sean nulos, los valores de dichos
ángulos son:
(1)
donde
(2)
(3)
de esta manera hemos obtenido una matriz B más sencilla que A con sus mismos valores propios.
En el método de Jacobi, se ha de encontrar una sucesión de matrices ortogonales P tal que hagan la matriz A diagonal.
(4)
En la práctica, se determina cada transformación Pi de forma que el elemento de mayor valor absoluto fuera de la diagonal se transforme en cero. El proceso se concluye cuando el máximo valor absoluto de los elementos fuera de la diagonal de B puede ser considerados despreciables frente al valor absoluto de los elementos de la diagonal principal.
Las columnas de la matriz
(5)
son los vectores propios asociados a los valores propios 
que constituyen los elementos de la matriz diagonal B.
Este es el procedimiento para hallar los valores propios reales y sus correspondientes vectores propios de una matriz simétrica A. La tarea ahora será la de codificar dicho procedimiento:
Veamos el procedimiento para hallar el valor absoluto mayor de los elementos que están fuera de la diagonal de una matriz simétrica. Consideremos una matriz cuadrada de dimensión 4,
los elementos marcados con d, son los de la diagonal principal, los elementos marcados con una x, son los que tenemos que examinar para encontrar el elemento (k, l) que tiene el mayor valor absoluto. Partimos del elemento que está situado el segundo en la primera fila k=0, l=1, y recorremos mediante un doble bucle el resto de los términos marcados con la x.
k=0; l=1;
maximo=Math.abs(a.x[k][1]);
for(i=0; i<n-1; i++){
for(j=i+1; j<n; j++){
if(Math.abs(a.x[i][j])>maximo){
k=i; l=j;
maximo=Math.abs(a.x[i][j]);
}
}
Una vez localizados los elementos k y l, se determina la matriz ortogonal P, que inicialmente es la matriz unidad.
Matriz p=new Matriz(n);
//...
for(i=0; i<n; i++){
for(j=0; j<n; j++){
p.x[i][j]=0.0;
if(i==j){
p.x[i][j]=1.0;
}
}
}
Luego, establecemos los valores de los elementos pkk, pkl, plk y pll, de acuerdo con las fórmulas (1) (2) y (3).
y=a.x[k][k]-a.x[l][l];
x=2*a.x[k][l];
z=Math.sqrt(x*x+y*y);
c=Math.sqrt((z+y)/(2*z));
s=signo(x/y)*Math.sqrt((z-y)/(2*z));
p.x[k][k]=c;
p.x[l][l]=c;
p.x[k][l]=s;
p.x[l][k]=-s;
En el caso de que y fuese cero se produciría un error de división por cero en el cálculo del seno. Ese caso se produce cuando akk es igual all, es decir cuando el seno y el coseno del ángulo q tiene el mismo valor, esto es, cuando q es 45º o p/4.
if(Math.abs(y)<CERO){
c=s=Math.sin(Math.PI/4);
Queda ahora definir el criterio por el que decidimos que dos números decimales son
casi iguales, o bien que su diferencia es próxima a cero. Decimos que un valor es cero o
muy próximo a cero cuando es menor que la constante CERO que hemos tomado como 10-8.
La nueva matriz transformada se halla realizando la operación 
: la matriz A por la matriz traspuesta de P, el
resultado se multiplicada por P. El resultado de este doble producto se vuelve a
guardar en la matriz A.
a=Matriz.producto(p, Matriz.producto(a, Matriz.traspuesta(p)));
En cada ciclo se obtiene la matriz ortogonal P, y se transforma A, guardándose de nuevo el resultado en A. Después de un número grande de ciclos, la matriz A estará prácticamente diagonalizada y los elementos de la diagonal principal serán los valores propios de la matriz simétrica A original.
Para obtener los vectores propios hemos de obtener una matriz Q resultado de multiplicar todas las matrices traspuestas P, véase la fórmula (5). Inicialmente Q es la matriz unidad y se multiplica por la traspuesta de P, y se guarda el resultado del producto de nuevo en Q.
q=Matriz.producto(q, Matriz.traspuesta(p));
Un bucle no puede repetirse indefinidamente, tiene que terminarse bajo ciertas condiciones. Tenemos dos condiciones de terminación del bucle, la que primero se cumpla:
En el segundo parámetro maxIter de la función valoresPropios, pasamos el número máximo de iteraciones. En un bucle do ... while la variable contador, que inicialmente es cero, se incrementa en una unidad en cada ciclo. El proceso se interrumpe cuando la variable contador se iguala a maxIter.
En este caso, la función valoresPropios emitirá un mensaje comunicando al usuario que los valores propios y sus correspondientes vectores no se han podido calcular, con la precisión deseada. Usaremos el mecanismo de las excepciones que nos proporciona el lenguaje Java para este propósito
public Matriz valoresPropios(double[] valores, int maxIter) throws ValoresExcepcion{
contador=0;
do{
contador++;
//...
}while(contador<maxIter);
if(contador==maxIter){
throw new ValoresExcepcion("No se han podido calcular los valores propios");
}
//...
}
}
class ValoresExcepcion extends Exception {
public ValoresExcepcion() {
super();
}
public ValoresExcepcion(String s) {
super(s);
}
}
|
La llamada a la función valoresPropios debe de efectuarse en un bucle try...
catch. Por ejemplo, si no se encuentran los valores y los vectores
propios después de 20 iteracciones se lanzará (throw) una excepción
creándose un objeto de la clase ValoresExcepcion, que será capturado por el
bloque catch y se mostrará en la pantalla de texto un mensaje indicando
tal eventualidad.
try{
vect=val.valoresPropios(valores, 20);
}catch(ValoresExcepcion ex){
System.out.println("Al calcular los valores propios se ha producido una excepción\n "
+ex.getClass()+ " con el mensaje "+ ex.getMessage());
}
El máximo valor absoluto de los elementos fuera de la diagonal de A se guarda en la variable maximo. Esta variable se compara con la variable tolerancia, en el caso de que maximo sea menor que tolerancia se sale del proceso cíclico y se muestran los resultados del cálculo.
La tolerancia se define de acuerdo con el criterio siguiente: Se suman los cuadrados de los elementos de la diagonal principal y se guarda el resultado en la variable sumsq. Se halla la raíz cuadrada de dicha suma y se divide por el número n de términos. Se define la tolerancia como la diezmilava parte de dicha cantidad.
do{
sumsq=0.0;
for(i=0; i<n; i++){
sumsq+=a.x[i][i]*a.x[i][i];
}
tolerancia=0.0001*Math.sqrt(sumsq)/n;
if(maximo<tolerancia) break;
contador++;
//...
}while(contador<maxIter);
En primer lugar, una función puede devolver un tipo de dato, cuando se concluye el proceso de cálculo, retornándolo, es decir, mediante la sentencia return, o bien pasando dicho dato a la función. Se pasa por valor una referencia al objeto, de modo que el objeto puede modificarse en el curso de la llamada a al función.
Teniendo en cuenta este principio, podemos elegir las siguientes variantes:
En este caso, se ha elegido la primera variante, dejando al lector la prueba de las otras variantes.
public Matriz valoresPropios(double[] valores, int maxIter){
Matriz q=new Matriz(n);
//...
//valores propios
for(i=0; i<n; i++){
valores[i]=(double)Math.round(a.x[i][i]*1000)/1000;
}
//vectores propios
return q;
}
El código completo de la función valoresPropios es, el siguiente
public Matriz valoresPropios(double[] valores, int maxIter)throws ValoresException{
final double CERO=1e-8;
double maximo, tolerancia, sumsq;
double x, y, z, c, s;
int contador=0;
int i, j, k, l;
Matriz a=(Matriz)clone(); //copia de this
Matriz p=new Matriz(n);
Matriz q=new Matriz(n);
//matriz unidad
for(i=0; i<n; i++){
q.x[i][i]=1.0;
}
do{
k=0; l=1;
maximo=Math.abs(a.x[k][1]);
for(i=0; i<n-1; i++){
for(j=i+1; j<n; j++){
if(Math.abs(a.x[i][j])>maximo){
k=i; l=j;
maximo=Math.abs(a.x[i][j]);
}
}
}
sumsq=0.0;
for(i=0; i<n; i++){
sumsq+=a.x[i][i]*a.x[i][i];
}
tolerancia=0.0001*Math.sqrt(sumsq)/n;
if(maximo<tolerancia) break;
//calcula la matriz ortogonal de p
//inicialmente es la matriz unidad
for(i=0; i<n; i++){
for(j=0; j<n; j++){
p.x[i][j]=0.0;
}
}
for(i=0; i<n; i++){
p.x[i][i]=1.0;
}
y=a.x[k][k]-a.x[l][l];
if(Math.abs(y)<CERO){
c=s=Math.sin(Math.PI/4);
}else{
x=2*a.x[k][l];
z=Math.sqrt(x*x+y*y);
c=Math.sqrt((z+y)/(2*z));
s=signo(x/y)*Math.sqrt((z-y)/(2*z));
}
p.x[k][k]=c;
p.x[l][l]=c;
p.x[k][l]=s;
p.x[l][k]=-s;
a=Matriz.producto(p, Matriz.producto(a, Matriz.traspuesta(p)));
q=Matriz.producto(q, Matriz.traspuesta(p));
contador++;
}while(contador<maxIter);
if(contador==maxIter){
throw new ValoresExcepcion("No se han podido calcular los valores propios");
}
//valores propios
for(i=0; i<n; i++){
valores[i]=(double)Math.round(a.x[i][i]*1000)/1000;
}
//vectores propios
return q;
}
|
Hallar los valores propios y sus correspondientes vectores propios de la matriz simétrica
Fijando en 20 el número de iteraciones, llamamos a la función valoresPropios dentro de un bloque try..catch. Se pasa el array valores propios valores a la función miembro valoresPropios y ésta devuelve la matriz de los vectores propios que se guarda en la matriz cuadrada vectores.
double[][] val1={{7, -1, -1}, {-1, 5, 1},{-1, 1, 5}};
Matriz matriz=new Matriz(val1);
Matriz vectores=new Matriz(matriz.n);
double[] valores=new double[matriz.n];
try{
vectores=matriz.valoresPropios(valores, 20);
}catch(ValoresExcepcion ex){//...}
Se imprimen los valores y los vectores propios
for(int i=0; i<matriz.n; i++){
System.out.print(valores[i]+" , ");
}
System.out.println("");
System.out.println(vectores);
Los resultados se resumen en la Tabla
| Valores propios | Vectores propios | Vectores equivalentes |
| 8 | (0.816, -0.408, -0.408) | (2, -1, -1) |
| 5 | (0.577, 0.577, 0.577) | (1, 1, 1) |
| 4 | (0.0, -0.707, 0.707) | (0, -1, 1) |
