<a href="https://colab.research.google.com/github/alexmascension/ANMI/blob/main/notebook/T3.ipynb" target="_parent"></a>

# Tema 3: Aproximación de autovalores


```python
!pip install -r https://raw.githubusercontent.com/alexmascension/ANMI/main/requirements.txt
```


```python
from sympy import *
from sympy.matrices import Matrix as mat
from sympy.matrices import randMatrix
from sympy import symbols
import sympy

import numpy as np

from scipy.linalg import orth
```


```python
from anmi.genericas import norma, print_verbose

from anmi.T2 import factorizacion_QR
from anmi.T3 import matriz_krylov, sucesion_krylov, potencia_iterada, metodo_autovals_QR
```

### Sucesiones de Krylov
Sea $A$ una matriz (aplicación lineal) y $x$ un vector. Si aplicamos la multiplicación de $A$ por $x$ de manera iterada obtenemos una serie de vectores $\{x, Ax, A^2x, A^3x, \cdots\}$. Si $x$ no es un autovector de $A$, entonces esa sucesión tendrá $n$ (dimensión de $A$) vectores independientes. Si $x$ es un autovector, con su autovalor $\lambda$, entonces la sucesión de vectores será, $\{x, \lambda x, \lambda^2x, \lambda^3x, \cdots\}$. Estas sucesiones de vectores se llaman *sucesiones de Krylov*.

Por otra parte, por el teorema de Cayley-Hamilton se tiene que $A^nx$ tiene que ser una combinación lineal de los siguientes elementos de la sucesión, es decir:
$$(-1)^nA^n + a_{n-1}A^{n-1} + \cdots + a_1A + a_0I = 0$$

Luego si tomamos $a = \begin{bmatrix}a_0\\a_1\\ \cdots \\ a_n\end{bmatrix}$ se tiene que
$$(x|Ax|\cdots|A^{n-1}x)a = (-1)^{n+1}A^nx$$

Y si resolvemos $a$, entonces se tienen los coeficientes del polinómio característico $p(\lambda) = a_0 + a_1\lambda + a_2\lambda^2 + \cdots + a_n\lambda^n$$



```python
help(matriz_krylov)
```


```python
# EJERCICIO 26
A = mat([[1, 1, 1], [0, 2, 2], [3, -1, 0]])
x = mat([[1, 0, 0]]).T

m_krylov = matriz_krylov(A, x)
m_krylov
```


```python
sk, a = sucesion_krylov(A, x)
sk
```


```python
# EJEMPLO 15
A = mat([[2, -1, 0], [-1, 2, -1], [0, -1, 2]])
x = mat([[-1, 0, 1]]).T

matriz_krylov(A, x)
```

### Método de la potencia iterada

En el método de la potencia iterada, se aplica la matriz de krylov hasta una potencia determinada, $k$. Entonces, se tiene que 
$$\lim_{k \to \infty} \frac{A^kw}{A^{k-1}w} = |\lambda_1|$$
Es decir, el mayor autovalor.

Además, si tomamos $ B= A^{-1} $, tenemos que
$$\lim_{k \to \infty} \frac{B^kw}{B^{k-1}w} = \frac{1}{|\lambda_n|}$$
Donde $\lambda_n$ es el menor autovalor.


```python
help(potencia_iterada)
```


```python
A.eigenvals()
```


```python
x = mat([[-2, 0, 1]]).T

matriz_krylov(A, x, 17)
```


```python
x = mat([[-1, 0, 0]]).T

np.array(potencia_iterada(A, x, 30, devolver_ultimo=False)[:, -3:], dtype=float)
```


```python
x = mat([[-1, 0, 0]]).T

np.array(potencia_iterada(A, x, 300, devolver_ultimo=True), dtype=float)
```


```python
N(2+sqrt(2))  # Vemos que converge al mayor autovalor
```


```python
np.array(potencia_iterada(A**-1, x, 300, devolver_ultimo=True), dtype=float)
```


```python
1/N(2-sqrt(2))  # Y lo mismo con el menor
```


```python
# Si tomamos una matriz ortogonal, el metodo de la potencia no tiene validez 
# porque se requiere que haya autovalores dominantes, y en este caso los 
# autovalores tienen módulo 1.

dict_QR = factorizacion_QR(A)
Q = dict_QR['Q']

Q
```


```python
Q * Q.T
```


```python
Q.eigenvals()
```


```python
N(-1/2 + 3*sqrt(70)/70 + 3*sqrt(14)/28 + sqrt(5)/5 + sqrt(70)*I*sqrt(6*sqrt(14) + 20*sqrt(5) + 73)/140)
```


```python
matriz_krylov(Q, x, 5)
N(matriz_krylov(N(Q), x, 30), 4)
```


```python
potencia_iterada(N(Q), x, 100, devolver_ultimo=False)[:, -5:]  # No hay una convergencia
```

### Método QR

El método QR consiste en aplicar los siguientes pasos:

$$A^{(1)} = A$$

De ahí sacamos que 
$$A^{(1)}  = Q^{(1)}R^{(1)}$$

De ahí construimos:
$$A^{(2)}  = R^{(1)}Q^{(1)}$$

Y se cumple que $A^{(1)}$ y $A^{(2)}$ son semejantes, luego tienen los mismos 
autovalores.

Con ello se reitera el proceso, y se cumple que las matrices equivalentes 
construidas convergen a una matriz triangular superior. Los la diagonal de $A^{(k)}$ converge a los autovalores de $A$.


```python
help(metodo_autovals_QR)
```


```python
dict_QR = metodo_autovals_QR(A, n_iters=10)
```


```python
N(dict_QR['A'][-2], 3)
```


```python
N(dict_QR['A'][-1], 3)
```


```python
N(2- sqrt(2), 3), 2, N(2 + sqrt(2), 3), 
```


```python
A = mat([[1, 1, 1], [0, 0, 1], [0, 1, 1]])
dict_QR = metodo_autovals_QR(A, n_iters=30, verbose=False)
```


```python
N(dict_QR['A'][-15], 3)
```


```python
N(dict_QR['A'][-1], 3)
```


```python
[N(i, 3) for i in list(A.eigenvals().keys())]
```

## Ejercicios

### Ejercicio 27

Determinar las primera iteraciones generadas por el método de la potencia con normalización con norma infinito cuando se aplica a la matriz con vector inicial.


```python
A = Matrix([[0, -1, 1], [0, 1, -1], [-1, -1, 2]])
A
```


```python
x0 = Matrix([1, -1, 2])
```


```python
n_iters = 15
m_krylov_x, m_krylov_w = zeros(A.shape[0], n_iters), zeros(A.shape[0], n_iters)
m_krylov_x[:, 0] = x0
m_krylov_w[:, 0] = x0/max(x0)

for i in range(1, n_iters):
    kriv_i_x = A * m_krylov_w[:, i - 1]
    m_krylov_x[:, i] = kriv_i_x 
    m_krylov_w[:, i] = kriv_i_x / max(kriv_i_x) 

```


```python
m_krylov_x
```


```python
m_krylov_w
```


```python
[N(i) for i in (np.array(m_krylov_x[0, 1:]) / np.array(m_krylov_w[0, :-1])).tolist()[0]]
```


```python
# Probamos sin normalizar

n_iters = 15
m_krylov_x  = zeros(A.shape[0], n_iters)
m_krylov_x[:, 0] = x0

for i in range(1, n_iters):
    kriv_i_x = A * m_krylov_x[:, i - 1]
    m_krylov_x[:, i] = kriv_i_x 

m_krylov_x 
```


```python
[N(i) for i in (np.array(m_krylov_x[0, 1:]) / np.array(m_krylov_x[0, :-1])).tolist()[0]]
```

### Ejercicio 28

Realizar una iteración con el método QR para el cálculo de autovalores de la matriz $A$


```python
A = Matrix([[1, 1, 0], [2, 1, 0], [2, 0, 1]])
```


```python
metodo_autovals_QR(A, n_iters=1, verbose=True)
```


```python
dict_QR = metodo_autovals_QR(A, n_iters=15, verbose=False)
N(dict_QR['A'][-1])
```


```python
N(dict_QR['A'][-2])
```


```python
np.sort(solve(det(A - eye(3) * Symbol('lambda')), Symbol('lambda')))[::-1]
```

### Ejercicio 29
Sea $A$ con $\theta$ dado. Determinar el n-ésimo términa de la sucesión de Krylov asociada a $A$ y al vector $x$. Por medio de la sucesión de Krylov determia el polinomio característico de $A$.


```python
t = Symbol('theta')
A = Matrix([[1, 0, 0], [0, cos(t), -sin(t)], [0, sin(t), cos(t)]])
x0 = Matrix([1, cos(t), -sin(t)])
```


```python
ss = simplify(matriz_krylov(A, x0, n_iters=8))
ss
```

Así a priori tiene pintas de que $A^nx = (1, \cos(n-1)\theta, sin(n-1)\theta)^t$. 

Como $n=3$ solo necesitamos los primeros 4 términos de la sucesión de krylov para resolver el sistema y crear el polinomio característico.


```python
M = ss[:3, :3]
am = Matrix([Symbol('a0'), Symbol('a1'), Symbol('a2')])
rhs = ss[:3, 3]
```


```python
sol = simplify(M.inv() * rhs)
sol
```


```python
l = Symbol('lambda')
simplify(sol[0] + sol[1] * l + sol[2] * l ** 2 - l ** 3)
```

### Ejercicio 31

Si $\alpha$ es una constante, aproximar el autovalore de módulo máximo de la matriz $A$ usando el método de la potencia (sin normalización) como vector inicial $v$, probando las diferentes componentes.


```python
a = Symbol('alpha')
n = Symbol('n')
A = Matrix([[1, 0, 0], [0, -2*cos(a)**2, 2*cos(a)*sin(a)], [0, 2*cos(a)*sin(a), -2*sin(a)**2]])
A
```


```python
v = Matrix([1, 1, 1])
v
```


```python
pol, vecfinal = sucesion_krylov(A, v)
pol
```


```python
solve(pol, l)
```


```python
vecfinal
```


```python
# Vamos a ir desglosando 
mk = matriz_krylov(A, v, n_iters=7)
mk
```


```python
simplify(Matrix(np.array(mk[:, 1:]) / np.array(mk[:, :-1])))
```

a85e3f72b386d4ca40805d3e7aa1e12ccead16ec
ipynb
Jupyter Notebook
notebook/T3.ipynb
alexmascension/ANMI
9c51a497a5fa2650f1429f847c7f9df69271168b
BSD-3-Clause
2021-11-30T23:30:45.000Z
2021-11-30T23:30:45.000Z
notebook/T3.ipynb
alexmascension/ANMI
9c51a497a5fa2650f1429f847c7f9df69271168b
BSD-3-Clause
2021-04-11T20:39:15.000Z
2021-04-13T17:45:43.000Z
notebook/T3.ipynb
alexmascension/ANMI
9c51a497a5fa2650f1429f847c7f9df69271168b
BSD-3-Clause
2021-11-30T23:31:11.000Z
2021-11-30T23:31:11.000Z
Qwen/Qwen-72B
1. YES
2. YES