# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# + [markdown] id="view-in-github" colab_type="text"
# <a href="https://colab.research.google.com/github/alexmascension/ANMI/blob/main/notebook/T3.ipynb" target="_parent"><img src="https://colab.research.google.com/assets/colab-badge.svg" alt="Open In Colab"/></a>

# + [markdown] id="Vf3GioQwcepr"
# # Tema 3: Aproximación de autovalores

# + id="1hDZNaYFX936"
# !pip install -r https://raw.githubusercontent.com/alexmascension/ANMI/main/requirements.txt

# + id="dYo_JGn2TSxC"
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

# + id="l-V6kF8_X937"
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

# + [markdown] id="B3T7oD8E1bEY"
# ### 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$$
#

# + id="hWHzJKX-_fyl"
help(matriz_krylov)

# + id="YFMQq7yCwW6W"
# 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

# + id="ADu3VAd5wXZ0"
sk, a = sucesion_krylov(A, x)
sk

# + id="C0Ystsv66PVG"
# EJEMPLO 15
A = mat([[2, -1, 0], [-1, 2, -1], [0, -1, 2]])
x = mat([[-1, 0, 1]]).T

matriz_krylov(A, x)

# + [markdown] id="9W55409Jywkz"
# ### 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.

# + id="SiaUm0xL_fyn"
help(potencia_iterada)

# + id="GTIZmAOv2Uni"
A.eigenvals()

# + id="ZG--F3wJBtXr"
x = mat([[-2, 0, 1]]).T

matriz_krylov(A, x, 17)

# + id="moRETh9eBtax"
x = mat([[-1, 0, 0]]).T

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

# + id="52wGtGj_1A4t"
x = mat([[-1, 0, 0]]).T

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

# + id="3TExHlDA2aR9"
N(2+sqrt(2))  # Vemos que converge al mayor autovalor

# + id="5mRLmChE0Lsx"
np.array(potencia_iterada(A**-1, x, 300, devolver_ultimo=True), dtype=float)

# + id="fZeAlvoy0Qsj"
1/N(2-sqrt(2))  # Y lo mismo con el menor

# + id="aCtHNnOo32Y3"
# 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

# + id="DxNd9Bu25PY0"
Q * Q.T

# + id="4OTGpPI24XMP"
Q.eigenvals()

# + id="8JqObRQh4dVy"
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)

# + id="wQA0QfnG4u4b"
matriz_krylov(Q, x, 5)
N(matriz_krylov(N(Q), x, 30), 4)

# + id="vHbR8UGa4hgd"
potencia_iterada(N(Q), x, 100, devolver_ultimo=False)[:, -5:]  # No hay una convergencia

# + [markdown] id="cjhzEPbM5x7P"
# ### 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$.

# + id="mxXCIhve_fys"
help(metodo_autovals_QR)

# + id="pcOLWedR_AYp"
dict_QR = metodo_autovals_QR(A, n_iters=10)

# + id="5_VzEkvM_LkV"
N(dict_QR['A'][-2], 3)

# + id="gJ8ibIVd_HH7"
N(dict_QR['A'][-1], 3)

# + id="8dWLVfTX_lAz"
N(2- sqrt(2), 3), 2, N(2 + sqrt(2), 3), 

# + id="9AYS1YyAAAU9"
A = mat([[1, 1, 1], [0, 0, 1], [0, 1, 1]])
dict_QR = metodo_autovals_QR(A, n_iters=30, verbose=False)

# + id="IKEdAQ_CAAU_"
N(dict_QR['A'][-15], 3)

# + id="MtaeFkK0AAVA"
N(dict_QR['A'][-1], 3)

# + id="79oglbStAAVB"
[N(i, 3) for i in list(A.eigenvals().keys())]

# + [markdown] id="9IX9k4ZTSSiO"
# ## Ejercicios

# + [markdown] id="dxhzBKn7SSrD"
# ### 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.

# + id="qHvk_qhbSp4D"
A = Matrix([[0, -1, 1], [0, 1, -1], [-1, -1, 2]])
A

# + id="8_31S9zcSqNb"
x0 = Matrix([1, -1, 2])

# + id="MiAEZjhFTGsv"
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) 


# + id="jIGHZaZFSvz4"
m_krylov_x

# + id="zeLdxnm0Udkl"
m_krylov_w

# + id="-dR0ujFSU1bi"
[N(i) for i in (np.array(m_krylov_x[0, 1:]) / np.array(m_krylov_w[0, :-1])).tolist()[0]]

# + id="OwcyXCkwVlhT"
# 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 

# + id="5L6GS871VPUF"
[N(i) for i in (np.array(m_krylov_x[0, 1:]) / np.array(m_krylov_x[0, :-1])).tolist()[0]]

# + [markdown] id="zIdRJzsoV_Xc"
# ### Ejercicio 28
#
# Realizar una iteración con el método QR para el cálculo de autovalores de la matriz $A$

# + id="w059mn3BWMjj"
A = Matrix([[1, 1, 0], [2, 1, 0], [2, 0, 1]])

# + id="3ISZtAmqWYh1"
metodo_autovals_QR(A, n_iters=1, verbose=True)

# + id="JU4FJ9AkWplI"
dict_QR = metodo_autovals_QR(A, n_iters=15, verbose=False)
N(dict_QR['A'][-1])

# + id="s-9D_Cb_W8BJ"
N(dict_QR['A'][-2])

# + id="0GfeTbWJXEwu"
np.sort(solve(det(A - eye(3) * Symbol('lambda')), Symbol('lambda')))[::-1]

# + [markdown] id="gk5nafPlXZ8L"
# ### 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$.

# + id="IOeqbEldXuE9"
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)])

# + id="CoJNBxDIYD_L"
ss = simplify(matriz_krylov(A, x0, n_iters=8))
ss

# + [markdown] id="_l6W_2WVSSyi"
# 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.

# + id="8Gk6eFrWZF9l"
M = ss[:3, :3]
am = Matrix([Symbol('a0'), Symbol('a1'), Symbol('a2')])
rhs = ss[:3, 3]

# + id="baxYJ0ekZgv7"
sol = simplify(M.inv() * rhs)
sol

# + id="QUiHi5OFZ2lF"
l = Symbol('lambda')
simplify(sol[0] + sol[1] * l + sol[2] * l ** 2 - l ** 3)

# + [markdown] id="E7tNXMTwsQwe"
# ### 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.

# + id="1eqi22XlsgMw"
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

# + id="nhfb3jZksxHS"
v = Matrix([1, 1, 1])
v

# + id="Yl5fxN4Ps2c8"
pol, vecfinal = sucesion_krylov(A, v)
pol

# + id="IplEaEeYyCGE"
solve(pol, l)

# + id="PIymkNwntGSZ"
vecfinal

# + id="tW_NB3tNtJNO"
# Vamos a ir desglosando 
mk = matriz_krylov(A, v, n_iters=7)
mk

# + id="mZZ0JSKlyRxa"
simplify(Matrix(np.array(mk[:, 1:]) / np.array(mk[:, :-1])))

StarcoderdataJupyterScriptsDedupFiltered