3. Valores y Vectores Característicos 3.1. Introducción El producto de una matriz cuadrada, A, por un vector (matriz columna), x, es otro vector, cuyas componentes son habitualmente no proporcionales a x. Sin embargo, puede existir un vector φ no nulo tal que: A φ = λ φ
(3.1a)
Se dice entonces que φ es un vector característico (también llamado vector propio, eigenvector o modo) de la matriz A. El correspondiente escalar λ es un valor característico (también llamado valor propio, autovalor o eigenvalor). Nótese que si un vector satisface las ecuaciones (3.1a) también un múltiplo arbitrario (un vector "paralelo") es solución. solución. Sin embargo, embargo, se trata esencialmente esencialmente de la misma solución; los vectores vectores característicos sólo se califican como distintos si sus componentes no son proporcionales. Por ejemplo, 2 − 1
− 1 1
1 = 1 2 1 1
2 − 1
− 1 1 1 = 3 2 − 1 − 1
1 1 en este caso φ1 = y φ 2 = son vectores característicos, a los que corresponden 1 − 1 los valores propios 1 y 3, respectivamente. Otros vectores, no paralelos a los dos antes mencionados, no cumplen la condición (3.1a): 2 − 1
− 1 1 0 = 2 2
3
0 1 El vector no puede expresarse como un múltiplo de . 3 2 El problema clásico de de valores y vectores característicos consiste en la determinación de los vectores φ y los correspondientes escalares λ para los que se cumple (3.1a). Con frecuencia se presenta el problema general: A φ = λ B φ
(3.1b)
En muchas aplicaciones las matrices A y B son simétricas y definidas positivas. positivas. En algunos casos se hacen hipótesis simplificadoras que resultan en B diagonal. El problema clásico, definido en (3.1a), corresponde al caso particular B = I. 3.1.1 Conversión del Problema General a la Forma Clásica
Un problema de la forma general (3.1b) puede convertirse a otro equivalente de la forma clásica (3.1a). Así por ejemplo, si B es no singular puede determinarse: -1
B A φ = λ φ
H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
(3.2)
3-1
Sin embargo, si A y B son simétricas (como es, por ejemplo, el caso en problemas de vibración, en los que esas matrices son respectivamente rigideces y masas) conviene más hacer la descomposición (Cholesky): T
B = R R
(3.3a)
y efectuar entonces el cambio de variables φ = R-1 z
(3.3b)
con lo que se obtiene: (R-1)T A R-1 z = λ z
(3.3c)
Esto es particularmente fácil si B es diagonal. ½
B = B B
½
φ = B-½ z B
-½
Donde hij =
AB a ij
-½
(3.4) z = H z = λ z
.
bi b j
Nótese que los valores característicos son los mismos que los del problema original; los correspondientes vectores característicos se relacionan mediante (3.4b). 3.1.2 Polinomio Característico y Valores Propios
Las ecuaciones A φ = λ B φ pueden también rescribirse como: (A - λ B) φ = 0
(3.5a)
que tiene soluciones no triviales sólo si la matriz (A - λ B) es singular, es decir, si: p(λ ) = det ( A − λB ) = 0
(3.5b)
p(λ ) se denomina polinomio característico . Siendo A y B matrices cuadradas de orden n, p(λ ) es un polinomio de grado n, cuyas raíces son λ1, λ2,
lo que sigue se λn . En lo
supone, sin perder generalidad, que: λ1 ≤ λ 2 ≤ λ 3 ≤ ⋅ ⋅ ⋅ ≤ λ n 3.1.3 Independencia Lineal de los Vectores Característicos
Asociado a cada uno de los n valores característicos λ i se tiene un vector φi . Si λ i es una raíz de multiplicidad m, el correspondiente vector φi puede obtenerse resolviendo el sistema de ecuaciones homogéneas: (A - λ i B) φi = 0 suponiendo m componentes arbitrarias en φi . Los vectores característicos correspondientes a valores característicos distintos son linealmente independientes. Supóngase que éste éste no fuera el caso, pudiéndose obtener uno de los vectores como combinación lineal de otros que sí son linealmente independientes: j
φ s =
∑ c φ
i i
(3.6a)
i =1
Y entonces:
H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
3-2
Sin embargo, si A y B son simétricas (como es, por ejemplo, el caso en problemas de vibración, en los que esas matrices son respectivamente rigideces y masas) conviene más hacer la descomposición (Cholesky): T
B = R R
(3.3a)
y efectuar entonces el cambio de variables φ = R-1 z
(3.3b)
con lo que se obtiene: (R-1)T A R-1 z = λ z
(3.3c)
Esto es particularmente fácil si B es diagonal. ½
B = B B
½
φ = B-½ z B
-½
Donde hij =
AB a ij
-½
(3.4) z = H z = λ z
.
bi b j
Nótese que los valores característicos son los mismos que los del problema original; los correspondientes vectores característicos se relacionan mediante (3.4b). 3.1.2 Polinomio Característico y Valores Propios
Las ecuaciones A φ = λ B φ pueden también rescribirse como: (A - λ B) φ = 0
(3.5a)
que tiene soluciones no triviales sólo si la matriz (A - λ B) es singular, es decir, si: p(λ ) = det ( A − λB ) = 0
(3.5b)
p(λ ) se denomina polinomio característico . Siendo A y B matrices cuadradas de orden n, p(λ ) es un polinomio de grado n, cuyas raíces son λ1, λ2,
lo que sigue se λn . En lo
supone, sin perder generalidad, que: λ1 ≤ λ 2 ≤ λ 3 ≤ ⋅ ⋅ ⋅ ≤ λ n 3.1.3 Independencia Lineal de los Vectores Característicos
Asociado a cada uno de los n valores característicos λ i se tiene un vector φi . Si λ i es una raíz de multiplicidad m, el correspondiente vector φi puede obtenerse resolviendo el sistema de ecuaciones homogéneas: (A - λ i B) φi = 0 suponiendo m componentes arbitrarias en φi . Los vectores característicos correspondientes a valores característicos distintos son linealmente independientes. Supóngase que éste éste no fuera el caso, pudiéndose obtener uno de los vectores como combinación lineal de otros que sí son linealmente independientes: j
φ s =
∑ c φ
i i
(3.6a)
i =1
Y entonces:
H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
3-2
j
Αφ s =
∑
j
c i Αφ i =
i =1
∑c
i
λι
Βφ i
(3.6b)
i =1
Por otro lado, por la definición del problema, (3.1b): j
Αφ s = λs Β φ s =
∑ c
i λs Β φ i
(3.6c)
i =1 j
Restando (3.6b) de (3.6c) se obtiene:
∑c ( i
λs
− λ i ) Βφ i = 0
i =1
Si
λs
≠ λ i debería
entonces tenerse ci = 0 para todo i , lo que se opone a la hipótesis.
Excepcionalmente pueden presentarse valores característicos repetidos. Aún en este caso es factible obtener vectores característicos linealmente independientes. independientes. Sin embargo, el conjunto de vectores asociados a los valores característicos repetidos define un subespacio, tal que cualquier vector del subespacio (es decir una combinación lineal de aquellos tomados como base) es también un vector característico: A φ i = λ i B φi A φ i = λ i B φi
(3.7)
A ( c1 φ1 + c2 φ2 + c3 φ3 + … ) = λ i B ( c1 φ1 + c2 φ2 + c3 φ3 + … )
Teniéndose n vectores característicos linealmente independientes de dimensión n, estos constituyen una base completa. Cualquier otro otro vector de tal dimensión dimensión puede expresarse como combinación lineal de los vectores característicos: v = α1 φ1 + α2 φ2 + α3 φ3 + … + αn φn
(3.8)
Por ejemplo, con los vectores característicos antes obtenidos: 1 3 1 1 1 = 2 − 2 2 1 − 1 3.1.4 Ortogonalidad de los Vectores Vectores Característicos
Si las matrices A y B son Hermitianas (o simplemente simétricas) y definidas positivas, los valores característicos de A φ = λ B φ son todos todos reales y positivos. positivos. Para probar esto basta considerar: * * φ s Αφ r = λ r φ s Bφ r
(3.9a)
φ *r Αφ s = λ s φ *r Bφ s
(3.9b)
El superíndice * denota aquí conjugada traspuesta. La conjugada transpuesta de la segunda de estas expresiones es (recuérdese que λ s es un escalar): φ *s Αφ r = λ*s φ *s Bφ r
(3.9c)
y al ser A y B Hermitianas (es decir A* = A y B* = B), restando (3.9c) de (3.9a) se obtiene:
(λ
r
− λ*s ) φ *s Bφ r = 0
(3.9d)
Si r=s, al ser B una matriz definida positiva se tendría φ *r Bφ r > 0 . Por lo tanto, siendo λ r = λ s , se tendría λ r − λ*r = 0 lo que implica que todos los λ son números reales. Si H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
3-3
además A es definida positiva, es decir φ *r Aφ r > 0 , se concluye que los valores característicos son todos positivos. Por otro lado, si λ r ≠ λ s se tiene que λ r − λ*s ≠ 0 y en consecuencia (3.9d) implica que: φ *s B φ r = br δ rs (es decir, cero si r ≠ s )
(3.10a)
y además, observando las expresiones precedentes: φ *s A φ r = a r δ rs
(3.10b)
Las propiedades de ortogonalidad expresadas en (3.10) son la base para la descomposición modal utilizada al resolver sistemas de ecuaciones diferenciales en aplicaciones tales como el análisis sísm ico lineal. Refiriéndose nuevamente al ejemplo inicial: 2
(1 1) − 1
− 1 1
=2
2 1
2
− 1 1 =6 2 − 1
(1 − 1) − 1 2
(1 1) − 1
− 1 1 =0 2 − 1
3.1.5 Normalización de los Vectores Característicos
Como se mencionó anteriormente los vectores característicos se definen por la proporción de sus elementos, pudiéndose escalar o "normalizar" en forma arbitraria. Es frecuente escalarlos de modo que: φ *s B φ r = δ rs
(3.11a)
Se dice entonces que los vectores están normalizados respecto a la matriz B. En tal caso se tiene también: φ *s A φ r = λ r δ rs
(3.11b)
3.1.6 Cociente de Rayleigh
Si se conoce un vector característico φi , el correspondiente valor λ i puede determinarse con el cociente de Rayleigh: ρ(φ i ) =
φT i Αφ i B φT φ i i
(3.12)
Esta expresión puede aplicarse también con aproximaciones a los vectores propios. Si x es una aproximación a un vector característico con un error de orden ε, el cociente de Rayleigh, ρ( x ), aproxima el correspondiente valor característico con un error de orden ε 2. 3.1.7 Teorema de Gershgorin
Supóngase que λ i es un valor característico de la matriz A y que φ i es el correspondiente vector, con componentes v1 v 2 v3 L :
H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
3-4
A φi = λ i φi
(3.13a)
La componente de mayor valor absoluto en φ i es v s . Dividiendo la ecuación s en (3.13a) entre v s e intercambiando ambos miembros: v1 v v + a s 2 2 + L + a ss + L + a sn n v s v s v s
λ i = a s1
(3.13b)
y por lo tanto: λ i − a ss = a s1 + a s 2 + L + 0 + L + a sn
(3.13c)
En consecuencia, cada valor característico λ i está dentro de por lo menos uno de los círculos con centro en a ss y radio igual a la suma de los valores absolutos de la correspondiente fila s. Por ejemplo, considerando la matriz: 2
− 1
A =
4
− 1
que es definida positiva, puede asegurarse que sus valores característicos (que son números reales) están dentro de los intervalos (1,3) y (3,5). Efectivamente, en este caso λ =3± 2 . 3.1.8 Formas polinómicas
Supóngase que se conocen los valores y vectores característicos de una matriz, A: A φ = λ φ
(3.14a)
¿Cuáles son los valores característicos de la matriz A2 = A A?
(AA) φ = A (A φ) = A (λ φ) = λ (A φ)= λ2 φ Este resultado puede extenderse para la matriz Ak (siendo k un exponente). Los vectores característicos son los mismos que los de la matriz A, mientras que los correspondientes valores característicos son λk : A φ = λ φ k
k
(3.14b)
Esto es incluso válido para exponentes negativos. Por ejemplo, multiplicando ambos miembros de (3.15a) por λ-1A-1 se obtiene: -1
-1
A φ = λ φ
(3.14c)
Por otro lado, combinando linealmente expresiones de la forma (3.14b) y teniendo en cuenta que A0 = I (así como λ0 = 1): (c0 I+ c1 A+ c2 A2 + c3 A3 +...) φ = (c0 + c1 λ+ c2 λ2 + c3 λ3 +...) φ
(3.14d)
Por ejemplo, si: 2
A =
− 1
− 1
2
tiene valores característicos 1 y 3, la matriz: 5
A = AA = 2
− 4
− 4
5
H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
3-5
tiene valores característicos 1 y 9 (es decir, los cuadrados de 1 y 3). Los vectores característicos son los mismos para ambas matrices.
3.2 Métodos de Iteración con Vectores Los métodos que se presentan en esta sección son los más eficientes cuando sólo se requieren un valor característico y su vector asociado, o en todo caso cuando el número de valores y vectores característicos por determinar es pequeño. 3.2.1 Iteración Directa
En la iteración "directa" se considera un vector inicial vectores corregidos, x k , mediante:
x 0 y
se obtiene una secuencia de
B x j +1 = A x j
x j +1 =
x j +1 r j +1
(3.15a)
(3.15b)
donde r j +1 es un escalar que normaliza el vector utilizado en la iteración. Lo habitual es tomar r j +1 como el elemento de máximo valor absoluto en x j +1 , lo que significa escalar el vector de aproximación de modo que la mayor componente sea igual a 1.. Este proceso converge al vector característico φ n , asociado al valor característico de mayor módulo, λ n . En efecto, la aproximación inicial x0 puede escribirse como: x0 = α1 φ1 + α2 φ2 + α3 φ3 + … + αn -1 φn -1 + αn φn
(3.16a)
Recuérdese que los n vectores característicos son linealmente independientes y constituyen una base completa en el espacio de dimensión n. Entonces (suponiendo que B no es singular): A x0 =
∑ α Aφ = ∑ α λ B φ i
i
i
i
i
x1 = B −1 Ax 0 = (α1 λ1) φ1 + (α2 λ2) φ2 + …+ (αn λn ) φn
(3.16b) (3.16c)
y por lo tanto: x1 =
1 r 1
∑ (α λ ) φ i
i
i
(3.16d)
n
Se observa que, si las componentes de x0 eran αi, aquellas de x1 resultan proporcionales a αi λi. Repitiendo pasos análogos a los indicados en (3.18), puede comprobarse que la aproximación xk puede expresarse como combinación lineal de los vectores característicos con coeficientes proporcionales a α i λ k i (en este caso k es un exponente). En consecuencia, si λ n ≥ λ n −1 ≥ λ n −2 ≥ L ≥ λ 1 , las componentes según φn crecen más rápidamente que las otras y se tiene que: Lim x k = φn
(3.17a)
Lim r k = λ n
(3.17b)
k → ∞
k →∞
Esto es válido aún cuando αn = 0 puesto que, por lo menos al tratar con grandes matrices, los errores de redondeo (debidos a la aritmética imperfecta del computador) introducen siempre una componente según φn . La convergencia es muy rápida si H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
3-6
o si x0 es aproximadamente paralelo a φn (es decir, si la componente αn es importante en relación a las demás). En cambio, si los últimos valores característicos son similares la convergencia es en general muy lenta. Por otro lado, no se tienen dificultades para el caso (más académico que práctico) en que λ n = λ n −1 : en tal caso el proceso converge a un vector característico que resulta ser la proyección de x0 en el subespacio definido por los vectores φn y φn -1. λ n >> λ n −1
Considérese por ejemplo el problema A φ = λ B φ con las matrices: 5 − 2 0 1 3 − 1 A = − 2 B = 0 0 0 − 1 1 Aún cuando en este caso se tienen matrices aplica a matrices cuadradas cualesquiera.
0
0
2
0
0
3
simétricas, el procedimiento descrito se
En este caso se obtienen: k
x k
Ax k
x k +1
r k+1
0
1.00000
5.00000
5.00000
5.00000
0.00000
-2.00000
-1.00000
0.00000
0.00000
0.00000
1.00000
5.40000
5.40000
-0.20000
-2.60000
-1.30000
0.00000
0.20000
0.06667
1.00000
5.48148
5.48148
-0.24074
-2.73457
-1.36728
0.01235
0.25309
0.08436
1.00000
5.49887
5.49887
-0.24944
-2.76370
-1.38185
0.01539
0.26483
0.08828
1.00000
5.50259
5.50259
-0.25130
-2.76994
-1.38497
0.01605
0.26735
0.08912
1.00000
5.50339
5.50339
-0.25169
-2.77128
-1.38564
0.01620
0.26789
0.08930
1.00000
5.50356
5.50356
-0.25178
-2.77156
-1.38578
0.01623
0.26801
0.08934
1
2
3
4
5
6
ρ (xk+1)
5.481481 5.40000 5.502594 5.48148 5.503559 5.49887 5.503603 5.50259 5.503605 5.50339 5.503605 5.50356 5.503605
1.00000 El procedimiento converge al valor característico: φ3 = − 0.25180 0.01623 que corresponde al valor característico de mayor módulo, λ3 = 5.503605.
H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
3-7
El valor de r es aproximadamente λn, pero el cociente de Rayleigh, ρ(x)" proporciona siempre una aproximación mejor. 3.2.2 Iteración Inversa
El proceso de iteración directa antes descrito converge al vector característico asociado al valor característico de mayor módulo. Éste puede ser útil al considerar el condicionamiento de las matrices de coeficientes en grandes sistemas de ecuaciones, o al analizar la estabilidad numérica de ciertos métodos para integrar sistemas de ecuaciones diferenciales, pero por lo general tiene poca importancia en la respuesta del sistema estudiado. Para determinar la respuesta de un sistema se requieren más bien los valores característicos de menor módulo y sus vectores asociados. Para determinar el vector característico asociado al valor propio de menor módulo (el modo fundamental) puede usarse una "iteración inversa": A x j +1 = B x j
x j +1 =
x j +1 r j +1
(3.18a)
(3.18b)
En este caso si: x0 = α1 φ1 + α2 φ2 + α3 φ3 + … + αn-1 φn-1 + αn φn
(3.19a)
la aproximación xk puede expresarse como combinación lineal de los vectores característicos con coeficientes proporcionales a α i λk i (nuevamente, k es aquí un exponente): xk =
α1 λk 1
α2
φ1 +
λk 2
φ2 +
α3 λk 3
φ3 + … +
α n −1 λk n −1
φn -1 +
αn λk n
φn
(3.19b)
En consecuencia, si λ1 ≤ λ 2 ≤ λ 3 ≤ ⋅ ⋅ ⋅ ≤ λ n al emplear la iteración inversa se tiene que: Lim x k = φ1
(3.20a)
k → ∞
Lim r k = k →∞
1
λ1
(3.20b)
Los comentarios anteriores relativos a la convergencia de la iteración directa son también válidos. En este caso la velocidad de convergencia depende de la razón λ2 / λ1. Para las matrices del caso anterior y considerando, por ejemplo, el vector inicial: 0 x 0 = 1 2
se obtiene el vector asociado al valor característico de menor módulo, es decir, λ1. Nótese que r es ahora una aproximación de 1 / λ1, mientras que en la iteración directa lo era de λn. También en este caso se observa que el cociente de Rayleigh es siempre una mejor aproximación al valor característico.
H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
3-8
k
x k
Bx k
x k +1
r k +1
0
0.00000
0.00000
2.66667
12.66667
1.00000
2.00000
6.66667
2.00000
6.00000
12.66667
0.21053
0.21053
1.42105
0.52632
1.05263
3.44737
1.00000
3.00000
6.44737
0.22041
0.22041
1.42993
0.53469
1.06939
3.46463
1.00000
3.00000
6.46463
0.22119
0.22119
1.43102
0.53594
1.07187
3.46696
1.00000
3.00000
6.46696
0.22128
0.22128
1.43116
0.53610
1.07221
3.46727
1.00000
3.00000
6.46727
0.22129
0.22129
1.43118
0.53613
1.07225
3.46731
1.00000
3.00000
6.46731
0.22129
0.22129
1.43118
0.53613
1.07226
3.46731
1.00000
3.00000
6.46731
1
2
3
4
5
6
ρ (xk +1)
0.154734 6.44737 0.154625 6.46463 0.154624 6.46696 0.154624 6.46727 0.154624 6.46731 0.154624 6.46731 0.154624
En muchas aplicaciones B es diagonal y A no lo es, por lo que la iteración directa es más simple. Sin embargo, un paso típico de la iteración inversa requiere aproximadamente el mismo número de operaciones que un paso de iteración directa. Supóngase que se tienen matrices de orden n y que A es una matriz de alta densidad (es decir, con pocos coeficientes no significativos). El número de operaciones requeridas para efectuar el producto Ax es de orden n2. Aquí se cuenta como una "operación" la combinación de una multiplicación o división con una suma o resta. También se ha supuesto que n es grande, por lo que n2 es mucho mayor que n. La división de Ax entre los coeficientes (de la diagonal principal) de B requiere un número de operaciones de orden n, que puede despreciarse. Es interesante observar que si previamente se realizó (una sola vez) la factorización A = LU, la solución del sistema de ecuaciones Ax = b requiere también un número de operaciones de orden n2, mientras que el producto Bx demanda sólo n operaciones. Por otro lado, si la matriz A es de baja densidad y tiene un ancho de semibanda promedio m, tanto un producto de la forma Ax como la solución de las ecuaciones Ax = b requieren aproximadamente mn operaciones. 3.2.3 Traslación
La velocidad de convergencia de la iteración inversa depende de las razones 1 / λi . Si λ 2 ≈ λ 1 la convergencia es lenta; siendo en cambio muy rápida si λ 1 << λ 2 . La convergencia puede acelerarse mediante una "traslación" µ ≈ λ 1 : H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
3-9
A φ = λ B φ
(3.21a)
(A - µB) φ = (λ−µ) B φ
(3.21b)
Nótese que el nuevo sistema (3.21b) tiene los mismos vectores característicos que el sistema original (3.21a) y valores característicos λ i - µ. Desde el punto de vista del polinomio característico, se ha trasladado el origen: 150 ) µ µ ( p
100
50
0 0
1
2
3
4
5
6
7 µ
-50
Si µ ≈ λ 1 puede lograrse que: λ 1 − µ >> λ 2 − µ
1
y por tanto:
λ1 − µ
<<
1
λ2 − µ
con lo que la convergencia mejora en forma apreciable.
,
Para el ejemplo anterior, efectuando una traslación µ = 0.154 se tiene: 4.846 A − 0.154 B = − 2 0 y por iteración inversa:
−2
0
2.692
−1
−1
0.538
k
x k
Bx k
x k +1
r k+1
0
0.00000
0.00000
692.29
3129.01
1.00000
2.00000
1677.41
2.00000
6.00000
3129.01
0.22125
0.22125
354.79
0.53608
1.07216
859.55
1.00000
3.00000
1603.26
0.22130
0.22130
354.80
0.53613
1.07226
859.57
1.00000
3.00000
1603.29
1
2
ρ (xk+1)
0.000624 1603.26 0.000624 1603.29 0.000624
0.22129 Se obtienen: λ1 = 0.154 + 0.000624 = 0.154624 y φ1 = 0.53613 . 1.00000
El siguiente algoritmo usa el cociente de Rayleigh para efectuar la traslación. Iniciando el proceso con y 0 = B x 0 y µ0 = 0:
(A − µ k B ) x k +1 = y k
H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
3 - 10
y k +1 = B x k +1 T
µ k +1 = µ k +
(
x k +1 y k x T k +1 y k +1
T
y k +1 = x k +1 y k +1
)−
1 2
(3.22)
y k +1
En relación con las expresiones precedentes: Lim y k = B φ1
(3.23a)
Lim µ k = λ 1
(3.23b)
k →∞
k →∞
La convergencia es cúbica. 3.2.4 Determinación de Otros Vectores Característicos
En los párrafos precedentes se ha visto cómo mediante iteración directa o inversa pueden obtenerse φn o φ1 respectivamente. Podría determinarse un valor característico intermedio y su vector asociado por iteración inversa con una traslación adecuada; sin embargo, esto requeriría un procedimiento previo para definir la traslación. En los que sigue se describe la determinación de sucesivos vectores característicos aprovechando las condiciones de ortogonalidad para el caso en que las matrices A y B son simétricas. La idea básica consiste en iterar con vectores ortogonales a los previamente obtenidos. Desafortunadamente, el proceso acumula los errores de los vectores previos y cada nuevo vector se determina siempre con menos precisión que el anterior. En la práctica se observa que se pierde una cifra significativa por cada nuevo vector; por tanto, no es factible determinar por este método más de unos 10 vectores característicos. En algunas aplicaciones esto puede no ser suficiente. A partir de un vector arbitrario: v = α1 φ1 + α2 φ2 + α3 φ3 + ... + αn φn
(3.24a)
puede obtenerse un vector ortogonal a los vectores característicos ya conocidos haciendo uso de las relaciones de ortogonalidad: φi T B v = α1 φi T B φ1 + α2 φi T B φ2 + ... + αn φi T B φn φi T B v = αi φi T B φi
(3.24b)
αi = ( φi T B v ) / ( φi T B φi )
(3.24c)
es decir: Luego es suficiente restar de v los αi φi para obtener un vector que (salvo por la imprecisión en la aritmética no tiene componentes según los vectores característicos previamente hallados. Para el ejemplo antes tratado, suponiendo que se haya obtenido el primer vector característico: 0.221295029 φ1 = 0.536128843 1.000000000 H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
3 - 11
0 y considerando v = 1 se obtiene α1 = 0.295889928, de donde: 0 − 0.06548 x0 = v - α1 φ1 = 0.84136 − 0.29589 es un vector ortogonal a φ1. Si se hace una iteración inversa con x0 se obtiene (suponiendo que se operara con una aritmética infinitamente precisa) el vector característico φ2: k
x k
Bx k
x k +1
r k +1
0
-0.06548
-0.06548
0.24319
0.64072
0.84136
1.68272
0.64072
-0.29589
-0.88767
-0.24696
0.37956
0.37956
0.40775
1.00000
2.00000
0.82960
-0.38544
-1.15631
-0.32671
0.49150
0.49150
0.43668
1.00000
2.00000
0.84594
-0.39382
-1.18147
-0.33553
0.51620
0.51620
0.44210
1.00000
2.00000
0.84715
-0.39663
-1.18990
-0.34275
0.52187
0.52187
0.43603
1.00000
2.00000
0.82913
-0.40460
-1.21379
-0.38466
1
2
3
4
ρ (xk +1)
1.20534 0.82960 1.17649 0.84594 1.17517 0.84715 1.17504 0.82913 1.17113
Es importante hacer notar que, como consecuencia de los errores de redondeo se introducen en las aproximaciones x j componentes según los vectores característicos originalmente eliminados.. En los resultados precedentes se tienen las siguientes componentes según φ1: k
α1
0
-1.565 x 10
-6
1
-1.580 x 10
-5
2
-0.000123
3
-0.000941
4
-0.007188
5
-0.056063
Como estas componentes tienden a crecer más rápidamente que la propia solución, es necesario eliminarlas cada 4 ó 5 pasos, utilizando el mismo proceso inicial:
H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
3 - 12
j −1
x j = v -
∑α
i
(3.25)
φi
i =1
Para el caso del ejemplo: 0.52588 0.221295029 0.53829 x 0 = 1.00000 − (− 0.056093) 0.536128843 = 1.03006 − 0.46393 1.000000000 − 0.40787 y luego de escalar este vector: k
x k
Bx k
5
0.52258
6
7
x
k+1
r k+1
0.52258
0.44489
0.85094
1.00000
2.00000
0.85094
-0.39597
-1.18790
-0.33696
0.52282
0.52282
0.44496
1.00000
2.00000
0.85098
-0.39599
-1.18796
-0.33698
ρ (xk+1)
1.17511 0.85098 1.17511
0.52288 1.00000 -0.39599
0.52288 = 1.17511 y φ2 = 1.00000 se obtienen: λ 2 = 0.85098 − 0.39599 1
3.2.5 Deflación
~
~
Otra alternativa es hacer una deflación , obteniendo un nuevo sistema A φ = λ B φ , de orden menor, con los mismos valores característicos del problema original, excepto los previamente determinados. En lo que sigue se aplica esta idea a un problema de la forma clásica H φ = λ φ . Considérese una matriz ortogonal, P , cuya última columna es igual al vector característico φ1 previamente determinado: P = ( p1
p2
L
p n −1
φ1 )
(3.26)
Al hacer el cambio de variable φ = P z se obtiene H P z = λ P z y premultiplicando por T T P : ( P H P ) z = λ z . Sin embargo, al ser la última columna de P igual a φ1 y suponiendo que ese vector haya sido normalizado de modo que φ1T φ1 = 1 se tiene: ~ H P HP = 0 T
0
λ 1
(3.27)
Esta matriz tiene los mismos valores característicos que la matriz original, H . Lo mismo ~ se puede decir de H , excepto por λ1 . Hay múltiples posibilidades para formar P . En el proceso propuesto por Rutishauser se hacen las operaciones equivalentes a trabajar con P = J1 J 2 L J n −1 , donde:
H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
3 - 13
1 J k =
O
c k
s k
− s k
c k O
columna
k
1
fila k fila k + 1
(3.28a)
columna k + 1
Nótese que P T H P = J T n−1 L J T 2 ( J 1T H J 1 ) J 2 L J n −1 , pudiéndose evaluar fácilmente los sucesivos productos, ya que en cada caso sólo se alteran dos filas y dos columnas. Los coeficientes ck y sk se determinan a partir de las componentes x1 x2 x3 L del vector característico previamente hallado, φ1 . Definiendo: qk 2 = x12 + x22 + x32 + L xk 2
(3.28b)
se tiene: s k = c k =
q k q k +1 x k +1
(3.28c)
q k +1
Para el ejemplo considerado anteriormente, sería necesario primero convertir el problema a la forma clásica. Al ser B diagonal: H=B
−1
2
AB
−1
2
5 = − 1.414214 0
− 1.414214 1.5
− 0.408248
− 0.408248 0.333333 0
Y para H φ = λ φ pueden obtenerse (por iteración inversa): Para H φ = λ φ pueden obtenerse (por iteración inversa): 0.11625 φ1 = 0.39827 0.90987
Luego: q1 = 0.11625
s1 = 0.28019
c1 = 0.95994
q 2 = 0.41489
s 2 = 0.41489
c 2 = 0.90987
q3 = 1.00000
Con el propósito de observar que, efectivamente, la última columna de P es igual a φ1 se está evaluando aquí la referida matriz: 0.95994 P = − 0.28019 0
0.28019 0.95994 0
0 1
0 0 1 0
0 0.90987 0.414898
− 0.41489 0.90987 0
es decir:
H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
3 - 14
0.95994 P = − 0.28019 0
0.25494
0.11625
0.87342
0.39827
0.90987
− 0.41489
de donde: 5.48594 T P H P = − 0.27561 0
− 0.27561 1.19271 0
0 0.15462 0
Nótese el λ 1 en la esquina inferior derecha. ~ 5.48594 H = − 0.27561
Los valores característicos de
− 0.27561
son λ 2 = 1.1751 y λ 3 = 5.5036 , es decir, iguales a los
1.19271
restantes valores característicos del problema original. Los correspondientes vectores resultan: 0.0638 0.9980 z2 = y z 3 = 0.9980 − 0.0638
de donde: φ k = B
El factor B original.
−1
2
−1
2
z k 0
P
se requiere para obtener los vectores del problema general en su forma
3.3 Métodos de Transformación Los métodos de este grupo son eficientes sólo si se requieren todos o una alta proporción de los valores y vectores característicos. La idea básica de estos procesos consiste en hacer un cambio de variables: φ = P z
para transformar
(3.29a) A φ = λ B φ
en:
( P-1A P ) z = λ ( P-1B P ) z
(3.29b)
Este sistema tiene los mismos valores característicos que el sistema original y vectores propios relacionados por (3.29a). Si las transformaciones son tales que las nuevas matrices tienen valores y vectores característicos fáciles de determinar, se ha resuelto indirectamente el problema original. 3.3.1 Método de Jacobi
El método de Jacobi (1846) puede considerarse como prototipo de los métodos de transformación. En este procedimiento se transforma el problema original a uno de la forma: a1
a2
O
z1 b1 b2 z 2 =λ M a n z n
O
z1 z 2 M bn z n
H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
(3.30)
3 - 15
que tiene como vectores característicos las columnas de la matriz identidad y como valores característicos los λ i = a i bi . Los valores característicos del sistema original son los mismos. P es en este caso una matriz ortogonal: -1
T
P = P
(3.31)
cuyas columnas son la propia solución buscada. proceso iterativo que se describe a continuación.
Ésta se determina mediante un
En la forma que aquí se presenta, este método se aplica a problemas de la forma clásica, A φ = λ φ, siendo A una matriz simétrica (real). Más adelante se consideran las modificaciones requeridas para problemas de la forma general. Empezando con A(0) = A y llamando φ(0) a los vectores característicos del problema original, el paso k del proceso se define como: φ(k ) = Pk φ(k +1) (k +1)
A(k ) ( Pk φ
(3.32a)
) = λ ( Pk φ(k +1) )
y si P es una matriz ortogonal, premultiplicando por PT se obtiene: ( Pk TA(k ) Pk ) φ(k +1) = λ φ(k +1) Lo que equivale a considerar un problema similar al original: (
)
(
)
(
)
A k+1 φ k +1 = λ φ k +1
(3.32b)
T
(3.32c)
Siendo: A(k +1) = Pk A(k ) Pk
Nótese que se mantiene la simetría de la matriz. Los valores característicos de esta nueva matriz son los mismos de la matriz original; los correspondientes vectores se relacionan por expresiones de la forma (3.32a). En el método de Jacobi las matrices Pk corresponden a una rotación plana: col i
1 Pk =
col j
O
cos θ k
− sen θ k
sen θ k
cos θ k O
1
(3.33)
fila i fila j
El objetivo de un paso es hacer cero un coeficiente aij = a ji. fácilmente que:
(
a ij +1 = a ji +1 = a jj − a ii (k
)
(k
)
(k)
(k)
) cos θ
k )
k
sen θ k + aij(
(cos
2
Puede verificarse
θ k − sen 2 θ k ) = 0
(3.34a)
y por tanto: tg 2θ k =
2 a ij( k ) k )
a ii(
− a jj(k )
0 ≤ θ k ≤
π 4
(3.34b)
Sólo los elementos de dos filas y de dos columnas (i , j ) se alteran en cada paso. Además, como se mantiene la simetría de la matriz A sólo deben calcularse los
H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
3 - 16
coeficientes de la submatriz triangular superior (o inferior) de A(k+1). Con la notación c = cos θ k ; s = sen θ k : aii +1 = a ii c 2 + 2a ij cs + a jj s 2 (k
)
(k)
(k)
(k)
(k +1 ) a jj = a ii( k ) s 2 -2aij( k ) cs + a jj( k ) c 2 (k +1 )
aij
)
(k +1 ) a jr
(3.35a)
= a ji(k +1 ) = (a jj(k) − a ii(k) ) cs + a ij( k ) (c 2 − s 2 ) = 0
( ) a ir +1 = a ri +1 = a ir c + a jr s (k
(k
=
)
(k +1 ) a rj
(k)
=−
k
(k) a ir s
+
( k ) a jr c
(3.35b)
En un cierto paso se hacen cero los elementos aij y a ji. Sin embargo, las sucesivas rotaciones reintroducen valores significativos en estas posiciones, por lo que es necesario repetir el proceso en varios "ciclos" para todos los elementos de fuera de la diagonal principal. El proceso es convergente. Si en un ciclo dado los cocientes
[a ] ( k ) ij
γ ij =
k a ii( )
2
k a jj( )
(3.36)
son de orden ε, éstos se reducen a orden ε2 en el siguiente ciclo. El número de ciclos completos necesarios para que la matriz A sea suficientemente aproximada a una matriz diagonal depende del orden de la matriz. Para matrices de orden 50 ó 60 pueden ser necesarios 8 a 10 ciclos. Cada ciclo demanda O(2n3) operaciones. Desde un punto de vista teórico sería más eficiente hacer cero los elementos aij en orden decreciente de los γ ij , definidos por (3.36), pero las comparaciones necesarias son relativamente lentas. Por eso se prefiere seguir un orden fijo en la selección de los elementos y efectuar las rotaciones sólo si γ ij es mayor que una tolerancia, variable en función del número de ciclo, m (por ejemplo 10-2m). La convergencia del proceso se puede verificar con una medida similar. Para determinar los vectores característicos es suficiente efectuar el producto de las matrices Pk ya que: φ(k ) = Pk φ(k +1)
(3.37a)
y por lo tanto: φ = φ(0) = P1 P2 P3 ... Pm
(3.37b)
Para ilustrar el método de Jacobi considérese el problema A φ = λ φ con:
A ( 0)
2 - 3 = 1 0
-3
1
0
6
-3
-3
6
1 − 3
1
−3
4
En el primer paso se hacen cero los coeficientes a12 y a21. precedentes: ( 0)
a11 = 2
(0)
a 22 = 6
(0)
En las expresiones
(0 )
a12 = a 21 = −3
H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
3 - 17
tg (2θ ) =
2 (− 3) 2−6
= 1.5 ⇒ cos θ = 0.881675 sen θ = 0.471858
0.881675 0.471858 P1 = 0 0
A
(1)
-0.471858
0
0
0.881675
0
0
0
1
0
0
0.394449 0 T (0) = P1 A P1 = -0.533899 0.471858
0
1 0
-0.533899
0.471858
7.60555
-3.11688
0.881675
-3.11688
6
-3
0.881675
-3
4
Luego se hacen cero los coeficientes a31: (1)
a11 = 0.394449
cos θ = 0.995574
(1)
a33 = 6
(1)
sen θ = 0.093978
0.995574 0 P2 = 0.0939783 0
A ( 2)
(1)
a13 = a 31 = −0.533899
0
- 0.0939783
0
1
0
0
0
0.995574
0
0
0.344051 -0.292919 = P2T A (1) P2 = 0 0.187835
0
1
-0.292919
0
0.187835
7.60555
-3.10309
0.881675
-3.10309
6.0504
-3.03107
0.881675
-3.03107
4
Nótese que se tienen nuevamente valores significativos en las posiciones otro lado: 0.877773 0.469770 P1 P2 = 0.0939783 0
- 0.471858
- 0.0828582
0
0.881675
- 0.0443444
0
0
0.995574
0
0
1
12 y 21.
Por
0
Procediendo en forma similar: (2)
a 22 = 7.60555
cos θ = 0.788374
A
( 3)
( 2)
a33 = 6.0504
( 2)
(2 )
a 23 = a 32 = −3.10309
sen θ = 0.615196
0.344051 - 0.23093 = P3T A (2) P3 = - 0.180203 0.187835
- 0.23093
- 0.180203
0.187835
10.027
0
0
3.62895
2.55979 - 1.84721
2.55979
- 1.84721
4
H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
3 - 18
0.877773 0.46977 P1 P2 P3 = 0.0939783 0 ( 3)
- 0.321026
- 0.355609
0
0.72237
0.507443
0
- 0.612474
0.784885
0
0
( 3)
a11 = 0.344051 a 44 = 4
cos θ = 0.99869
A ( 4)
( 4)
-0.0854342
10.027
0
0
3.62895
2.54462
-1.85401
2.54462 -1.85401 4.00963 0
0.0449207
0.722370
0.507443
0.0240408
- 0.612474
0.784885
0
0
( 4)
0.00480941
0.99869
( 4)
a 24 = a 42 = 2.54462
sen θ = −0.343849
0.334426 -0.339576 = P5T A (4) P5 = -0.0854342 0.124345
a33 = 3.62895
-0.361627
- 0.355609
( 4)
0.876622 0.469154 P1 P2 P3 P4 P5 = 0.0938551 - 0.0511758 ( 5)
( 3)
- 0.321026
a 44 = 4.00963
cos θ = 0.939025
1
a14 = a 41 = 0.187835
0.334426 -0.361627 = P4T A (3) P4 = -0.0854342 0
a 22 = 10.027
0
sen θ = 0.0511758
0.876622 0.469154 P1 P2 P3 P4 = 0.0938551 - 0.0511758
A ( 5)
( 3)
-0.339576
-0.0854342
0.124345
10.9588
-0.6375
0
-0.6375
3.62895
0
-1.74096
-1.74096 3.07785
- 0.286006
- 0.355609
0.152566
0.686590
0.507443
- 0.225811
- 0.573474
0.784885
0.343398
0
0.937795
(5)
a 44 = 3.07785
( 5)
0.215115
( 5)
a34 = a 43 = −1.74096
cos θ = 0.760371 sen θ = 0.649489
A (6)
0.334426 -0.339576 = P6T A (5) P6 = -0.145722 0.0390598
-0.339576
-0.145722
0.0390598
10.9588
-0.484737
-0.414049
-0.484737
5.11603
0
-0.414049
0
0.876622 0.469154 Φ ≈ P1 P2 P3 P4 P5 P6 = 0.0938551 - 0.0511758
1.59076
- 0.286006
- 0.369485
- 0.114957
0.686590
0.532507
- 0.573474
0.457089
0.157878 0.673341
0.343398
- 0.609087
0.713072
con lo que termina un primer “ciclo”. Análogamente, al terminar el segundo ciclo:
H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
3 - 19
A
(12 )
0.317649 - 0.007656 = - 0.0006033 - 0.0003418
- 0.007656
- 0.0006033
- 0.0003418
11.0268
0.0015155
0.0015155
5.08272
0.0000149 0
0.0000149
0
0.856314 0.505117 Φ ≈ P1 P2 L P12 = 0.0771368 - 0.0750522
1.57279
- 0.275908
- 0.421440
- 0.11397
0.619208
0.566580
0.201062
- 0.640155
0.399623
0.361467
- 0.584532
0.651578
0.722518
y al finalizar el tercer ciclo:
A (18)
0.317644 =
1.57279
11.0269 5.08272
No se muestran los coeficientes con valor absoluto menor que 10-6. Los coeficientes de la diagonal de A(18) son (aproximaciones a) los valores característicos de la matriz A. Nótese que no se obtienen en orden ascendente o descendente. Las columnas del producto P1 P2 L P18 son los correspondientes vectores, que se obtienen normalizados: T Φ Φ = I. Esto se comprueba fácilmente, ya que las matrices Pk son todas ortogonales. 0.856032 - 0.276628 - 0.421478 - 0.114202 0.505686 0.618991 0.566358 0.200923 Φ ≈ P1 P2 L P18 = 0.076907 - 0.640107 0.399777 0.651558 - 0.074671 0.361373 - 0.584614 0.722537 3.3.2 Caso de Matrices Hermitianas.
El método de Jacobi puede también emplearse para hallar los valores y vectores característicos de una matriz Hermitiana, H , cuyos coeficientes (en general complejos) tienen simetría conjugada. En este caso se hacen productos de la forma: H ( k +1) = U *k H ( k ) U k
(3.38)
en los que U k es una matriz unitaria, es decir, tal que U −k 1 = U *k (el superíndice * denota en este caso la conjugada traspuesta). Para hacer cero el coeficiente hij se utiliza: col i
1 U k =
col j
O
e
iθ
cos φ
− e i θ sen φ
sen φ
cos φ O
1
fila i
(3.39)
fila j
Suponiendo que: H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
3 - 20
hii(
k )
hij( k ) = b − ic
=a
( k )
h jj( k ) = d
h ji = b + ic
(3.40a)
Las partes real e imaginaria de los nuevos coeficientes (i, j ) resultan:
(d − a ) cos φ sen φ cos θ + b cos 2 φ − b sen 2 φ cos 2θ − c sen 2 φ sen 2θ = 0 (a − d ) cos φ sen φ sen θ + b cos 2 φ − c sen 2 φ sen 2θ − c sen 2 φ cos 2θ = 0 de donde: tan θ =
c b
tan 2φ =
2 (b cos θ + c sen θ )
(3.40b)
a − d
3.3.3 Método de Jacobi Generalizado.
Es posible modificar el método de Jacobi "clásico" antes descrito para resolver directamente el problema general A φ = λ B φ. En lo que sigue se supone que A y B son simétricas y que esta última es definida positiva (y posiblemente no diagonal). Debe anotarse que si B fuera diagonal sería más eficiente transformar el problema a la forma clásica. Un paso del proceso general se define por: T
A(k +1) = Pk A(k ) Pk
(3.41)
T
B(k +1) = Pk B(k ) Pk
donde Pk es una matriz similar a la utilizada para el proceso clásico: col i
1 Pk =
col j
O
1
α k
γ k
1 O
1
fila i
(3.42)
fila j
α y γ se determinan de: (k +1 )
= a ji(k +1 ) = α k aii(k) + (1 + α k γ k ) a ij( k ) + γ k a jj(k) = 0
(k +1 )
= b ji(k +1 ) = α k bii(k) + (1 + α k γ k ) bij( k ) + γ k b jj(k) = 0
aij bij
(3.43)
Estas dos ecuaciones son independientes, excepto en el caso en que a ii(k) bii(k)
=
aij(k) bij(k)
=
(k) a jj
b jj(k)
en el que puede considerase, por ejemplo:
(3.44a)
α k = 0
γ k =
aij(k) (k)
(3.44b)
a jj
Definiendo:
H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
3 - 21
k
k )
c1 = a jj( ) bij(
− b jj(k ) a ij( k )
c 2 = aii( k ) bij( k ) − bii(k ) a ij( k ) c3 =
1 2
(a
( k ) ii
b jj( k ) − bii( k ) a jj( k )
(3.45a)
)
d = c 3 + ( signo c 3 ) c 32 + c1c 2
se obtienen: γ k = −
c2
α k =
d
c1 d
(3.45b)
El radical en la expresión de d es siempre positivo si B es una matriz definida positiva. Puede observarse que si B fuera la matriz identidad se obtendrían: α k = − γ k = −tg θ . Los comentarios precedentes relativos a la convergencia son también aquí aplicables. El número de operaciones en cada ciclo es de O(3n3). El siguiente ejemplo ilustra los aspectos nuevos introducidos en esta sección. Se pretende determinar los valores y vectores característicos del sistema: A φ = λ B φ, donde: 1
2 1 B = − 1 1 1 2 Para estas pequeñas matrices con un paso es suficiente: A =
− 1
i=1, j=2 a11 = 1
a22 = 1
a12 = a21 = -1
b11 = 2
b22 = 2
b12 = b21 = 1
c1 = 3
c2 = 3
c3 = 0
d=3
α = -γ k =1
1 − 1 1 − 1 1 1 4 T = A(1) = P1 A(0) P1 = 1 1 1 1 1 1 − − 0 1 − 1 2 T B(1) = P1 B(0) P1 = 1 1 1
0
0
1 1
1 2 0 = 2 − 1 1 0 6
de donde: λ2 = 4/2 = 2 λ1 = 0/6 = 0
1
1 1
1 1 − 0
Φ = P1 diag (bi ½) =
0
2 1
6
0.7071
=
0.4082
− 0.7071 0.4082
La post multiplicación de P1 sólo es necesaria para escalar los vectores de modo que φ T i B φ j = δ ij . Al igual que en el procedimiento clásico los valores característicos (y los correspondientes vectores) no quedan necesariamente ordenados. 3.3.4 El Método QR
Este proceso se aplica al problema clásico A φ = λ φ, donde A no requiere ser simétrica, pudiendo tener valores característicos cero (o incluso negativos). En el caso H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
3 - 22
más general, para una matriz A cualquiera, el método QR es poco eficiente, ya que requiere O( 4 3 n3) operaciones por paso. Sin embargo, sólo se requieren O(4n2) operaciones por paso si A es de la forma Hessemberg: a11 a 21 0 A= 0 0 L
a12
a13
a14
a15
a 22
a 23
a 24
a 25
a 32
a 33
a 34
a 35
0
a 43
a 44
a 45
0
0
a 54
a 55
K
O
(3.46)
es decir si es casi triangular superior, excepto por una codiagonal inferior. Para el caso particular en que la matriz A es además simétrica (y por lo tanto tridiagonal): a1 b1 A=
b1 a2
b2
b2
a3
b3
b3
a4 O
O O
(3.47)
el método QR es aún más eficiente, requiriendo tan solo O(12n) por paso. En todo caso es siempre posible efectuar la transformación a la forma Hessemberg (tridiagonal si A y B son simétricas), requiriéndose un total de O( 5 3 n3) operaciones (una sola vez). Debe anotarse además que, a diferencia del método de Jacobi, el método QR mantiene la posible configuración banda de la matriz y permite efectuar traslaciones (análogas a las de una iteración inversa), tanto para acelerar la convergencia como para mejorar la precisión en los valores característicos de interés. El objetivo del proceso conocido como QR es la determinación de los valores característicos; conocidos estos, los correspondientes vectores pueden obtenerse por iteración inversa con traslación. Considerando A(0) = A, el paso básico del método QR consiste en hacer la descomposición: A ( k ) = Q k R k
(3.48a)
donde Q k es una matriz ortogonal (es decir, Q T k Q k = I ) y R k es una matriz triangular superior. Luego se efectúa el producto en orden cambiado: A ( k +1) = R k Q k
(3.48b)
Obsérvese que premultiplicando (3.48a) por Q T k se obtiene: ( k ) Q T = R k k A
(3.48c)
y por lo tanto: A(
k +1)
= R k Q k = Q T k A ( k ) Q k
H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
(3.48d)
3 - 23
Nótese que si A ( k ) es simétrica A ( k +1) también resulta simétrica. La expresión (3.48d) indica además que A ( k +1) es "similar" a A ( k ) : sus valores característicos son los mismos, los correspondientes vectores se relacionan por una transformación lineal: A ( k ) φ ( k ) = λ φ ( k )
(3.49a)
al efectuar el cambio de variables: φ ( k ) = Q k φ ( k +1)
(3.49b)
se obtiene: A
( k )
(Q
Q k φ
( k +1)
T ( k ) Q k k A
= λ Q k φ ( k +1)
)φ
( k +1)
= λ φ ( k +1)
(3.49c) (3.49d)
Ambas matrices tienen los mismos valores característicos (que en consecuencia son los de la matriz original) y vectores característicos relacionados por (3.49b). A medida que k crece A ( k ) converge a una matriz triangular superior (cuyos valores característicos son los elementos de la diagonal principal); para el caso simétrico A ( k ) converge a una matriz diagonal. Los valores característicos se obtienen en orden descendente; así la aproximación al valor característico de menor módulo se obtiene en la posición nn de la matriz. La convergencia del proceso es análoga a la de la iteración inversa. Cuando en pasos sucesivos se obtienen valores similares en el extremo inferior de la diagonal principal, puede afirmarse que se tiene una aproximación al primer valor característico. La convergencia puede acelerarse efectuando traslaciones: ( k ) µ k = a nn
(3.50a)
A ( k +1) = R k Q k − µ k I
(3.50b)
Nótese que los valores característicos de esta nueva matriz son iguales a los de la ( k ) matriz original menos la translación. Cuando se logra que a nn = 0 puede hacerse una traslación: µ k = a n( k −)1, n −1
(3.50c)
para mejorar la convergencia al segundo valor característico y análogamente se procede para los otros valores requeridos. Por regla general se requieren sólo 2 pasos por cada valor característico adicional. Al finalizar el proceso debe agregarse a los valores λ obtenidos la suma de las traslaciones µ k efectuadas. Los vectores característicos podrían obtenerse con el producto: φ ( 0 ) = Q 1Q 2 Q 3 L
(3.51)
pero este proceso es poco eficiente, siendo más conveniente obtener estos vectores por iteraciones inversas con traslaciones iguales a los valores característicos ya determinados. Esto permite también mejorar la precisión en los λ . La determinación de Q y R en un paso puede hacerse en diversas formas. El proceso más eficiente consiste en transformar A en una matriz triangular superior utilizando matrices de rotación plana (como en el método de Jacobi):
H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
3 - 24
(P
)
T T n , n −1 L P31
T P21 A=R
(3.52a)
y por lo tanto: Q = P21 P31 L Pn,n −1
(3.52b)
La matriz P ji , que permite hacer cero el coeficiente ji: col i
1 Pk =
col j
O
cos θ k
− sen θ k
sen θ k
cos θ k O
1
fila i
(3.53a)
fila j
se obtiene mediante: k )
cos θ =
sen θ = d =
aii(
d a ji( k ) d
(a ) ( k ) ji
2
(3.53b)
+ (a ii( k ) )
2
Sólo se requiere un ciclo de estas transformaciones para obtener iterar. 2 1 (0) Para un ejemplo del proceso considérese la matriz: A = 1 4 0 1
R. No es necesario
0
1 2
Esta es una matriz simétrica (lo cual no es un requisito para emplear el método QR) y, siendo tridiagonal, tiene la “forma Hessemberg”. Para transformar A en una matriz triangular superior R se hace primero cero el coeficiente a21: (0)
( 0)
d = 2.236068
cos θ = 0.894427
sen θ = 0.447214
a11 = 2
a 21 = 1
.894427 - .447214 0 P21 = .447214 .894427 0 0 1 0 T P21 A (0 )
2.236068 2.683281 0.447214 0 3.130494 .894427 = 1 2 0
Luego se hace cero a32, con lo que se obtiene una matriz triangular superior:
H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
3 - 25
(0)
a 22 = 3.130494
cos θ = 0.952579
( 0)
a32 = 1
d = 3.286335
sen θ = 0.304290
0 0 1 P32 = 0 .952579 - .304290 0 .304290 .952579
T T R 1 = P32 P21 A ( 0)
2.236068 2.683281 0.447214 = 0 3.286332 1.460532 0 1.632993 0
Y se completa el primer paso efectuando el producto: A (1)
0 3.200000 1.469694 = R 1Q1 = R 1 P21 P32 = 1.469694 3.244444 .496904 .496904 1.555556 0
Análogamente, en el segundo paso: .908739 - .417365 0 P21 = .417365 .908739 0 0 1 0 0 0 1 P32 = 0 .978097 - .208150 0 .208150 .978097
T T R 2 = P32 P21 A (1)
A (2)
L 3.521363 2.689686 0 2.387242 .765454 = 0 1.427490 0
0 4.322580 .996351 = R 2 Q 2 = R 2 P21 P32 = .996351 2.281193 .297132 .297132 1.396226 0
Y en el tercer paso: .974449 - .224610 0 P21 = .224610 .974449 0 0 1 0 0 0 1 P32 = 0 .989134 - .147017 0 .147017 .989134
T T R 3 = P32 P21 A (2)
4.435924 1.483271 .066739 0 2.021076 .491663 = 0 1.33850 0
H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
3 - 26
A
(3)
0 4.655737 .453953 = R 3 Q 3 = R 3 P21P32 = .453953 2.020318 .196780 .196780 1.323944 0
( 3) Suponiendo que el coeficiente a 33 = 1.323944 sea una buena aproximación al primer valor característico, se efectúa una traslación:
A (3)
0 3.331794 .453953 − 1.323944 I = .453953 .696375 .196780 .196780 0 0
obteniéndose en el cuarto paso: 0 3.405209 .088938 .678541 - .017396 A (4) = .088938 - .017396 - .055581 0 Se hace entonces una nueva traslación:
µ 4 = -.055581 ⇒ µ =
∑µ
k
= 1.268362
obteniéndose: 0 3.463559 .018800 (5) A = .018800 .731767 .000010 .000010 - .000413 0 y nuevamente:
µ 5 = -.000413 ⇒ µ =
∑µ
k
= 1.267949
obteniéndose: A (6)
3.464096 .003973 0 = .003973 .732056 0 0 0 0
Se observa ahora que el coeficiente a33 es menor que 10-6, lo que implica que λ1 es aproximadamente igual a la suma de las traslaciones previamente realizadas. Conviene luego hacer una traslación igual al resultado obtenido para a22 a fin de mejorar la precisión para el segundo valor característico: µ 6 = .732051 ⇒ µ =
A
(6)
∑µ
k
=2
0 2.732050 .003973 0 0 − 0.732051 I = .003973 0 - .732051 0
Y puede trabajarse con la submatriz de un orden menor: 2.732050 .003973 A (6) = .003973 0
H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
3 - 27
(6)
(6)
a11 = 2.732050
cos θ = 1
a 21 = 0.003973
d = 2.732051
sen θ = 0.000307
1.000000 - .000307 P21 = Q 7 = .000307 1.000000 2.732051
T R 7 = P21 A ( 6) =
0
.000840
− .000001
2.73205 0 0 0
A ( 7 ) =
Los coeficientes indicados como 0 son menores que 10-6. Los valores característicos de esta matriz son 0 y 2.732051. Para obtener aquellos de la matriz original deben sumarse las traslaciones: λ1 = -0.732051 + 2 = 1.267949 λ2 = 0 + 2 = 2 λ3 = 2.732051 + 2 = 4.732051 3.3.5. Transformación a la Forma Hessemberg
Si el método QR se aplicara a una matriz cualquiera sería en general poco eficiente, puesto que requiere O ( 4 3 n 3 ) operaciones por paso. Para reducir el número de operaciones a O( 4n 2 ) por paso debe previamente transformarse la matriz a la forma "Hessemberg" (es decir, una matriz que es casi triangular superior, teniendo además coeficientes significativos en la primera codiagonal inferior): h11 h12 h21 h22 0 h32 H= 0 0 M M 0 0
h14
L
h23
h24
L
h33
h34
L
h43
h44
L
M
M
O
0
0
hn,n−1
h13
h1n
h2 n h3n h4 n hnn
(3.54)
Si la matriz original fuera simétrica, la transformación a la forma Hessemberg, que puede hacerse conservando la simetría, produce una matriz tridiagonal. En tal caso el QR requiere apenas 12 n operaciones por paso. Cabe anotar que la forma Hessemberg (tridiagonal para el caso simétrico) no se pierde en los sucesivos pasos del método QR. La transformación a la forma Hessemberg sólo requiere hacerse una vez. Por lo tanto las O ( 5 3 n 3 ) que se gastan en la transformación están plenamente justificadas. Entre los procedimientos que se encuentran en la literatura para efectuar la transformación, se propone el cambio de variables φ = B φ , con lo que el problema original Α φ = λ φ se reescribiría como B −1 Α B φ = λ φ o bien H φ = λ φ . En este caso B −1 Α B = H o, lo que es lo mismo, Α B = B H . En el proceso original de Hessemberg se usa una matriz B de la forma: H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
3 - 28
0 0 1 0 0 0 0 1 0 b32 1 0 B= 0 b b43 1 42 M 0 bn 2 bn 3 bn 4
L 0
L 0
L 0 L 0
(3.55a)
1
L
con coeficientes arbitrarios en la primera columna (que por simplicidad se ha escrito como la primera columna de la matriz identidad). Los coeficientes de las sucesivas filas de H y columnas de B pueden entonces obtenerse con las expresiones: n
∑a
hir = a ir +
r
ik b kr
k = r +1
bi , r +1
−
∑b
ik h kr
i = 1, 2, L r + 1
(3.55b)
i = r + 2, L n
(3.55c)
k =1
n = a ir + a ik bkr − hr +1,r k = r +1
∑
1
r
∑ k =1
bik hkr
Este procedimiento podría fallar si en algún paso hr +1, r = 0 . El proceso podría recomenzarse con una primera columna de B diferente, lo que en general evitaría el error, aunque esto no puede garantizarse. Por otro lado, el procedimiento antes expuesto no mantiene la posible simetría de la matriz A . También puede hacerse la transformación a la forma Hessemberg por rotaciones planas (método de Givens) o reflexiones (Householder). El método de Householder utiliza matrices ortogonales y simétricas, de la forma: P = I − 2 w w T
(3.56)
donde w es un vector unitario: w T w = 1 . Es fácil probar que P = P T = P −1 . La matriz P refleja al espacio en el "plano" que pasa por el orígen y es ortogonal a w . Considérese un vector cualquiera v = α 0 w + α1u donde u T w = 0 . Entonces, T P v = ( I − 2 w w ) ( α 0 w + α 1u ) = − α 0 w + α 1u . Nótese que la componente según w ha cambiado de signo, es decir, el vector v ha sido reflejado en el plano ortogonal a w. La transformación de A en H mediante el método de Householder requiere n − 2 pasos ( n es aquí el orden del sistema) de la forma: A(
k +1)
= Pk A ( k ) Pk
(3.57a)
donde: Pk = I − θ k w k w T k
θ k =
2
w T k w k
(
( k )
(3.57b)
)
w k = v k + signo a k +1,k v k e k +1
siendo:
H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
3 - 29
0 M 0 ( k ) v k = a k +1, k a ( k ) k + 2,k M a ( k ) nk
e k +1
una matriz que contiene los coeficientes de la columna k de A ( k ) que están por debajo de la diagonal principal, y
0 M = 1 M 0
la columna k + 1 de la matriz identidad (de orden n ).
Para que el proceso sea más eficiente, debe observarse que al premultiplicar A , cuyas columnas son a1 a 2 a 3 L , por la matriz P , cada columna se modifica en forma independiente. Las columnas de A = PA resultan:
(
)
(
)
T a j = I − θ k w k w T k a j = a j − θ k w k a j w k
(3.58a)
Igualmente, al postmultiplicar A por P las filas se modifican en forma independiente. Llamando ahora a i a la fila i de la matriz A = PA , la correspondiente fila de A = A P resulta:
(
)
T ai = ai I − θ k w k w T k = a i − θ k ( a i w k ) w k
(3.58b)
Por ejemplo, considérese la matriz: 4 3 A= 2 1
3
2
1
4
3
2
3
4
2
3
=A 3 4
(1)
Transformación de la primera columna a la forma Hessemberg: 0 3 v1 = 2 1
v1 = 3.74166
0 0 0 3 1 6.74166 w 1 = v 1 + v1 e 2 = + 3.74166 = 2 0 2 1 0 1 θ1 w 1T A (1) = (1.00000
1.38619
1.23786
H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
θ1 = 0.03964
0.93096)
3 - 30
P1 A (1)
4.00000 - 3.74166 = 0 0
3
2
- 5.34522
- 5.34522
0.22762
1.52428
0.61381
1.76214
- 4.27618 1.13809 3.06904 1
1.00000 − 2.02190 (1) θ1 P1 A w 1 = 0.22681 0.42543
A ( 2)
4.00000 - 3.74166 = P1 A (1) P1 = 0 0
- 3.74166
0
8.28571
- 1.30143
- 1.30143
1.07067
- 2.25428
0.91128
- 2.25428 0.91128 2.64362 0
Transformación de la segunda columna a la forma Hessemberg: 0 0 v2 = 1 . 30143 − − 2.25428
v 2 = 2.60298
0 0 0 0 0 0 w 2 = v 2 + v 2 e3 = − 2.60298 = − 1.30143 1 − 3.90441 − 2.25428 0 − 2.25428
θ 2 w T 2 A ( 2 ) = (0
P2 A ( 2)
1.00000
4.00000 - 3.74166 = 0 0
θ 2 = 0.09840
- 0.93647 )
- 0.61346
- 3.74166
0
8.28571
- 1.30143
2.60298
- 1.32452
0
- 0.47162
- 2.25428 - 2.74510 0.53254 0
0 1 ( 2) θ 2 P2 A w 2 = 1.11774 0.06306
H = A ( 3)
4.00000 - 3.74166 = P2 A ( 2) P2 = 0 0
- 3.74166
0
8.28571
2.60298
2.60298
3.03959
0
- 0.22540
H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
0 - 0.22540 0.67470 0
3 - 31
3.4 Métodos Mixtos Los dos procesos que se describen en lo que sigue son adecuados para sistemas de orden grande en el caso en que se requieran muchos vectores característicos. 3.4.1 Iteración con la Determinante de (A - µ B)
Los valores propios de A φ = λ B φ son los ceros del polinomio característico p (λ ) = det (A − λ B ) = 0 . Por ejemplo, si: 4 2 A= 0 0
2
0
8
2
2
8
0
2
0
1 0 B= 0 0
0 2 4
0
0
0
2
0
0
0
2
0
0
0
1
Las raíces del polinomio: 4 − λ 2 p(λ ) = det 0 0
2
0
0
8 − 2λ
2
2
8 − 2λ
0 2
0
2
4 − λ
= 4 λ4 − 64 λ3 + 364 λ2 − 864 λ + 720 = 0
son los valores característicos. La determinación de los coeficientes del polinomio característico es factible (utilizando, por ejemplo, el método de Hessemberg). Una vez obtenidos los coeficientes del polinomio característico, se requiere determinar los valores de λ para los que p (λ ) = 0 . Sin embargo, éste es frecuentemente un problema mal condicionado: pequeños errores en los coeficientes causan grandes errores en las raíces. Por ello, los métodos en los que se hace una determinación explícita del polinomio característico sólo son adecuados para pequeñas matrices. Para matrices de orden elevado, pero con un ancho de banda comparativamente pequeño, pueden determinarse los valores característicos por iteración, evaluando la determinante de A − µ k B para una secuencia de valores µ k que se corrigen con procesos tales como el método de la secante. Así, dadas las aproximaciones µ k −1 y µ k a una raíz y habiéndose calculado p(µ k −1 ) = A − µ k −1 B y p(µ k ) = A − µ k B se obtiene una mejor aproximación, µ k +1 , mediante:
p(µ k ) ( ) ( ) p µ − p µ k k −1
µ k +1 = µ k − η
µ k − µ k −1
1≤ η ≤ 2
(3.59)
La evaluación de p (µ ) no requiere tener el polinomio p (λ ) en forma explícita. Si A − µ B se descompone en el producto de una matriz triangular inferior, L , con unos en la diagonal principal, por una matriz triangular superior, U , se tiene que: p(µ ) = det (A − µ B ) = det (LU ) = det (L ) det (U )
(3.60a)
donde: det (L ) = l11 l 22 l 33 l 44 L = 1 det (U ) = u11 u 22 u 33 u 44 L
H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
(3.60b)
3 - 32
y por lo tanto: p (µ ) = u11 u 22 u 33 u 44 L u nn La descomposición de A − µ B en factores triangulares LU requiere pocas operaciones si el ancho de banda es pequeño. Particularmente importante es el caso en el que las matrices A y B son simétricas y definidas positivas (todos los valores característicos son reales y positivos). En tal caso puede aplicarse la propiedad de Sturm: el número de coeficientes negativos en la diagonal principal de U al hacer la descomposición A − µ B = LU es igual al número de valores característicos menores que µ k . Esta propiedad, combinada con la iteración (3.59) u otra similar, permite obtener una primera aproximación a una raíz. Sin embargo, el proceso debe combinarse con iteraciones inversas usando el cociente de Rayleigh para refinar los valores obtenidos. Para las matrices A y B antes indicadas, con µ = 1.5 : 2.5 2 A − 1.5 B = 0 0
0 0 1 5 2 0 .800 1 0 = 2 5 2 0 .588 1 0 2 2.5 0 0 .523 2
0
0
0 2.5
0 0 0 0 1 0
2
0
0
3.4
2
0
0
3.824
0
0
2 1.454
p (1.5) = 2.5 ⋅ 3.4 ⋅ 3.824 ⋅ 1.454 = 47.25
Análogamente se obtienen:
λ1 λ2
λ3 λ4
Número de coeficientes negativos en la diagonal principal de U
µ k
p (µ k )
1.5
47.25
0
2.0
0
0
2.5
-8.75
1
3.0
0
1
3.5
11.25
2
4.0
16.00
2
4.5
11.25
2
5.0
0
2
5.5
-8.75
3
6.0
0
3
6.5
47.25
4
3.4.2 Iteración en Subespacio
El método tratado en la sección precedente es eficiente cuando las matrices tienen ancho de banda relativamente pequeño. Cuando el ancho de banda es grande es más adecuado un proceso de iteración en subespacio , como se describe en este acápite. Este método tiene por objeto determinar en forma simultánea los p vectores característicos asociados a los valores característicos de menor módulo. La idea básica es que es mucho más fácil iterar para obtener un subespacio que contenga a estos vectores que iterar para obtener cada uno de ellos por separado. Se trabaja con una colección de q vectores linealmente independientes ( q > p ). Los q vectores iniciales definen un subespacio que no necesariamente contiene a los p H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
3 - 33
vectores de interés. Si esos p vectores característicos si estuvieran contenidos en el subespacio, sería suficiente proyectar A φ = λ B φ para obtener el sistema A z = λ B z , de orden q << n , que sería fácil de resolver por métodos de transformación. Los valores característicos del problema proyectado serían los mismos del problema original, mientras que sus vectores característicos, z , corresponderían a las proyecciones de los vectores φ en el subespacio. No siendo éste el caso, se hacen iteraciones inversas para mejorar los q vectores con los que se trabaja, de modo que el subespacio por ellos definido sea más y más “paralelo” a los p vectores propios de interés. En lo que sigue, se supone que A y B son matrices simétricas. Siendo X k los q vectores de aproximación, en cada ciclo del proceso se realizan los pasos siguientes: a. Iteración inversa: AX k +1 = BX k
La matriz A debe factorizarse antes de iniciar las iteraciones. Los vectores X k +1 son más “paralelos” a los primeros p vectores característicos. b. Proyección de A y B en el subespacio definido por los vectores X k +1 : ( k +1)
= X T k +1 A X k +1
( k +1)
= X T k +1 B X k +1
A B
Las matrices A ( k +1) y B ( k +1) son cuadradas, simétricas, de orden q . c. Solución del problema de valores y vectores característicos proyectado: A ( k +1) Q k +1 = B ( k +1) Q k +1 Λ k +1
Λ k +1 es una matriz diagonal, cuyos coeficientes son los valores característicos del
problema proyectado. Si los X k +1 definen un subespacio que contiene a los p primeros vectores propios, los p menores valores en Λ k +1 son parte de la solución buscada. d. Determinación de nuevos vectores: X k +1 = X k +1 Q k +1
Como consecuencia de los pasos c y d: T
T
(
T
)
T
T
(
T
)
T
X k +1 A X k +1 = Q k +1 X k +1 A X k +1 Q k +1 = Q k +1 A ( T
X k +1 B X k +1 = Q k +1 X k +1 B X k +1 Q k +1 = Q k +1 B (
k +1)
k +1)
Q k +1 = Λ q Q k +1 = I q
es decir, los vectores X k +1 satisfacen las condiciones de ortogonalidad, lo que asegura que la iteración inversa no produce q vectores todos iguales a φ1 . Si en las X 0 hay componentes según todos los p vectores característicos de interés: Lim Λ k = diag λ1
λ2
L
λp
L
φ2
L
φp
L
k →∞
Lim X k = diag φ1 k →∞
Habiéndose obtenido en dos ciclos sucesivos los estimados λ( pk ) y λ( pk +1) para el mayor de los valores característicos requeridos, el cociente λ( pk +1) − λ( pk ) λ( pk +1) da una H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
3 - 34
medida adecuada del error relativo y es útil para verificar la convergencia. Adicionalmente, debe comprobarse que los valores y vectores obtenidos corresponden a los p menores valores característicos. Para ello puede usarse la propiedad de Sturm, factorizando A − µ B en LU con valores de µ ligeramente mayores a los λ calculados. Si A y B son simétricas, de orden n , ancho de semibanda m , y A es definida positiva, el número de operaciones iniciales requeridas es de O ( 12 nm 2 ) , esencialmente para la factorización de A . En cada ciclo de la iteración, considerando q << n , deben hacerse O ( nq ( 4m + 2q + 3)) operaciones. Esto puede reducirse a O ( nq ( 2m + 2q + 3)) cuando B es diagonal. Para el procedimiento como se ha descrito en los párrafos precedentes, se trabaja con q = min ( 2 p, p + 8) . Habitualmente unos 10 ciclos de iteración son suficientes para obtener 6 cifras significativas correctas en los p valores y vectores característicos. Las operaciones finales requieren O ( 12 nm 2 p ) operaciones adicionales. Aproximación Inicial
Para iniciar el proceso se requieren q vectores linealmente independientes, agrupados en X 0 . Si A y B fueran diagonales, los vectores característicos serían las columnas e k de la matriz identidad. Aún cuando A y B no sean diagonales, éste puede ser un buen criterio para construir la aproximación inicial X 0 . En particular, deberían escogerse las columnas cuyo índice k corresponde a los máximos bkk a kk . Con el propósito de introducir componentes según todos los vectores característicos, se acostumbra además considerar dos columnas con componentes arbitrarios (que podrían ser todos iguales a 1, o iguales a los bkk a kk ). En algunas aplicaciones es fácil obtener una buena aproximación al primer vector característico, por ejemplo, como solución de un sistema de ecuaciones de la forma A x 1 = b . Las sucesivas columnas x k para una excelente aproximación inicial pueden entonces obtenerse como vectores de Ritz, mediante un proceso recursivo que combina pasos de iteración inversa con ortogonalización: A y k = B x k −1
y T k B x j x x k = y k − x T jB x j j j =1 k −1
∑
Determinación de Grupos de Vectores Característicos Haciendo Traslaciones
Si se requieren muchos vectores característicos, el procedimiento estándar de iteración en subespacio puede hacerse más eficiente utilizando sucesivas traslaciones en combinación con procedimientos de eliminación de las componentes según los vectores ya conocidos. En este caso se trabaja con subespacios de dimensión q , con el propósito de determinar grupos de p ≈ q 2 vectores. Habitualmente q = máx ( 4, m ), siendo m el ancho (promedio) de semibanda. Para cada grupo de vectores, se realizan cómputos iniciales que incluyen: a. Determinación de la traslación (el proceso se inicia con µ = 0 ) H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
3 - 35
~
µ = 0.1 λ n + 0.9 λ n+1 λ n es el último valor característico para el que se ha logrado convergencia; ~ λ n +1 es la aproximación al siguiente valor característico.
b. Factorización: A − µB = L U c. Determinación de q vectores de aproximación inicial, X 0 . La iteración incluye los pasos siguientes: a. Eliminación de las componentes de X k según los vectores característicos previamente determinados (ver acápite 3.2.4). b. Iteración inversa: Yk +1 = B X k L U Z k +1 = Yk +1
c. Proyección de A − µB y B en el subespacio definido por los vectores Yk +1 : A
( k +1)
= Z T k +1 Yk +1
B ( k +1) = Z T k +1 B Z k +1
Las matrices A ( k +1) y B ( k +1) son cuadradas, simétricas, de orden q . d. Solución del problema de valores y vectores característicos proyectado:
(A
( k +1)
+ µ B ( k +1) )Q k +1 = B ( k +1) Q k +1 Λ k +1
e. Determinación de nuevos vectores: X k +1 = Z k +1 Q k +1
f. Verificación de la convergencia Como en el procedimiento estándar, debe verificarse que se tienen los valores característicos correctos utilizando la propiedad de Sturm. Ejemplo simple
Supóngase que se requieren dos vectores característicos de A φ = λ B φ , siendo: 2 −1 A= 0 0
−1
0
2
−1
−1
2
0
−1
0
0 − 1 2
0 B=
1 0
2
En este caso particular la iteración inversa produce en un solo paso el subespacio que incluye a los dos primeros vectores característicos, ya que dos de los valores característicos son infinitos. Para hacer más eficiente el proceso debe factorizarse primero la matriz A :
H. Scaletti - Métodos Numéricos: Valores y Vectores Característicos
3 - 36
