**Métodos computacionales 2**

**Universidad de los Andes**


```python
import numpy as np
```

## Eliminación gaussiana con pivoteo parcial (repaso)



```python
def triangular_superior(A):
  B = A.copy().astype(float)
  n = B.shape[0]
  for i in range(n):
    indi_max = np.argmax(np.abs(B[i:, i]))
    if indi_max > 0:
      C = B[i].copy()
      B[i] = B[i + indi_max]
      B[i + indi_max] = C
    for j in range(i+1, n):
      B[j] = B[j] - (B[j, i]/B[i, i]) * B[i]
  return B

def matriz_diagonal(A_t):
  B = A_t.copy().astype(float)
  n = B.shape[0]
  for i in range(n-1, -1, -1):
    for j in range(i-1, -1, -1):
      B[j] = (B[j, i]/B[i, i]) * B[i] - B[j]
    B[i] = B[i] / B[i, i]
  return B
```

## Inversión de matrices usando eliminación gaussiana

Dada una matriz $A$ es posible encontrar su inversa haciendo eliminación gaussiana sobre la matriz $A$ y copiando cada operación sobre la matriz identidad.

Ejemplo: Encontrar la inversa de $A$

\begin{equation}A=   \left (
      \begin{array}{rrr}
          2 &  1 &  3  \\
         -1 &  2 &  4 \\
          0 &  1 &  3 
      \end{array}
   \right )\end{equation}

   \begin{equation}\left (
      \begin{array}{rrr:rrr}
          2 &  1 &  3 & 1 & 0 & 0 \\
         -1 &  2 &  4 & 0 & 1 & 0\\
          0 &  1 &  3 & 0 & 0 & 1
      \end{array}
   \right )\end{equation}


```python
A = np.array([
    [ 2, 1, 3],
    [-1, 2, 4],
    [ 0, 1, 3],
]).astype(float)

def matriz_inversa(A):
  M_identidad = np.identity(A.shape[0])
  B = np.concatenate((A, M_identidad), axis = 1)
  print(B)
  B_t = triangular_superior(B)
  D = matriz_diagonal(B_t)
  print(D)
  return D[:, A.shape[0]:]

inversa_A = matriz_inversa(A)
print('inversa de A:')
print(inversa_A)
print('comprobar inversa:')
print(np.round(inversa_A @ A, 14))
print(np.round(A @ inversa_A, 14))
```

    [[ 2.  1.  3.  1.  0.  0.]
     [-1.  2.  4.  0.  1.  0.]
     [ 0.  1.  3.  0.  0.  1.]]
    [[ 1.    0.    0.    0.5   0.   -0.5 ]
     [-0.    1.   -0.    0.75  1.5  -2.75]
     [ 0.    0.    1.   -0.25 -0.5   1.25]]
    inversa de A:
    [[ 0.5   0.   -0.5 ]
     [ 0.75  1.5  -2.75]
     [-0.25 -0.5   1.25]]
    comprobar inversa:
    [[1. 0. 0.]
     [0. 1. 0.]
     [0. 0. 1.]]
    [[ 1.  0. -0.]
     [ 0.  1. -0.]
     [ 0.  0.  1.]]


## Cálculo de determinantes

llamemos $B$ a la matriz que resulta al realizar eliminación gaussiana sobre la matriz $A$, obteniendo una matriz triangular superior.

Recordemos que 

* Intercambiar dos filas multiplica el determinante por -1
* Multiplicar una fila por un escalar distinto de cero multiplica el determinante por el mismo escalar
* Agregar a una fila un múltiplo escalar de otra no cambia el determinante.

Entonces el determinante de A es el producto de la diagonal de B con $(-1)^d$ donde $d$ es el numero de pivoteos.

$$\det(A) = (-1)^d*\prod\operatorname{diag}(B)$$



```python
def triangular_superior(A):
  B = A.copy().astype(float)
  n = B.shape[0]
  d = 0
  for i in range(n):
    indi_max = np.argmax(np.abs(B[i:, i]))
    if indi_max > 0:
      C = B[i].copy()
      B[i] = B[i + indi_max]
      B[i + indi_max] = C
      d += 1
    for j in range(i+1, n):
      B[j] = B[j] - (B[j, i]/B[i, i]) * B[i] 
  return B, d

def determinante(A):
   B, d = triangular_superior(A)
   prod = 1
   for i in range(len(B)):
        prod = prod * B[i, i]
   prod = prod * (-1)**d
   return prod
```


```python
A = np.random.rand(3,3)

print('Determinante con nuestro método:',determinante(A))
print('Determinante con numpy         :',np.linalg.det(A))
```

    Determinante con nuestro método: 0.0821413592337298
    Determinante con numpy         : 0.08214135923372984


## Valores propios de una matriz

Dada una matriz $A$ diagonalizable de tamaño $n \times n$ sus valores propios se corresponden con las raíces del polinomio característico de $A$ que es

$$P(z) = \det(zI - A) = 0$$

Buscar valores propios de forma numérica por este método se considera completamente inadecuado, el solo hecho de encontrar el polinómio caracteristico engendra inconvenientes, sobre todo para matrices grandes. Por lo cual se han desarrollado técnicas, de las cuales solo veremos introducciones básicas que nos permitirán encontrar los valores propios de algunas matrices, pero no para cualquier matriz cuadrada general, dadas las limitaciónes de tiempo de la clase.   

## Método de las potencias

Permite conocer el valor propio mayor (en valor absoluto) de una matriz $A$ diagonalizable de tamaño $n \times n$. El método parte de algún vector $v_0$ que puede ser aleatorio y calcula iterativamente

$$v_{k+1} = \frac{Av_k}{\|Av_k\|}$$

Entonces $v_k$ converge normalmente al autovector de mayor autovalor.

Se asume que $A$ tiene un valor propio que es estrictamente mayor en magnitud que sus otros valores propios y el vector inicial $v_0$ tiene un componente distinto de cero en la dirección de un autovector asociado con el autovalor dominante.
La aproximación al valor propio dominante $\mu_k$ se puede calcular a partir del vector propio como
$$\mu_{k} = \frac{v_{k}^{\intercal}Av_{k}}{v_{k}^{\intercal}v_{k}}$$



**Ejemplo:** Encontrar el vector propio de mayor valor propio de:

\begin{equation}A = 
\left[\begin{matrix}
1 & 2 & 3\\
1 & 2 & 1\\
3 & 2 & 1\\
\end{matrix}\right]\end{equation}



```python
def metodo_de_potencias(A, v0, iteraciones:int):
  Av0 = A @ v0
  v = Av0 / (np.sqrt(sum(Av0 * Av0)))
  if iteraciones == 1: return v
  return metodo_de_potencias(A, v, iteraciones-1)

def valor_propio(A, v_prop):
  return (v_prop.T @ A @ v_prop) / (v_prop @ v_prop)

```


```python
A = np.array([
              [ 1, 2, 3],
              [ 1, 2, 1],
              [ 3, 2, 1]
])
v0 = np.array([1, 0, 0])
iteraciones = 20

v = metodo_de_potencias(A, v0, iteraciones)

v_teo = np.array([1, (1 + np.sqrt(5))*0.5 - 1, 1])
print('vector propio aprox:', v)
print('vector propio real:', v_teo/np.sqrt(v_teo @ v_teo), '\n')

mu = valor_propio(A, v)
print('valor propio dominante aprox:', mu)
print('valor propio dominante real: ', 3 + np.sqrt(5))
```

    vector propio aprox: [0.64793617 0.40044657 0.64793616]
    vector propio real: [0.64793616 0.40044657 0.64793616] 
    
    valor propio dominante aprox: 5.23606797749979
    valor propio dominante real:  5.23606797749979


## Algoritmo $QR$
Es un procedimiento para calcular los valores propios de una matriz. La idea básica es realizar una descomposición $QR$, escribiendo la matriz como un producto de una matriz ortogonal y una matriz triangular superior, multiplicar los factores en el orden inverso e iterar.

El algoritmo $QR$ consiste en  multiplicar los factores $Q$ y $R$ en el orden inverso para obtener $A_0$ que es el resultado de la primera iteración

$$A_0 = R_0 Q_0$$
En general
$$A_{k+1} = R_k Q_k$$
Se puede sustituir a $R_k$ haciendo

$$A_{k+1} = R_k Q_k = Q_k^{-1} Q_k R_k Q_k = Q_k^{-1} A_k Q_k $$

$$A_{k+1} = Q_k^{\mathsf{T}} A_k Q_k$$

Al iterar múltiples veces, los valores en la diagonal de $A_k$ se aproximan a los valores propios de $A$. En la forma como está planteado acá este algoritmo no puede encontrar valores propios complejos y tampoco funcionará para todos los casos con valores reales. La intención es simplemente dar algunas nociones que sirvan de introducción a este extenso tema.

## Descomposición $QR$

Consiste en una factorización de una matriz $A$ en un producto $A  =  QR$ de una matriz ortogonal $Q$ y una matriz triangular superior $R$. Existen varios métodos para calcular la descomposición $QR$, como por medio del proceso de Gram-Schmidt, transformaciones de Householder o rotaciones de Givens. Se explica a continuación el proceso de Gram-Schmidt. Podemos nombrar a las columnas de $A$ y $Q$ así:

\begin{equation}A = 
\left(\begin{matrix}
| & | &  & |\\
a_0 & a_1 & \cdots & a_{n-1}\\
| & | &  & |\\
\end{matrix}\right)\,\,\,\,\,\,\,\,\,\,Q = 
\left(\begin{matrix}
| & | &  & |\\
q_0 & q_1 & \cdots & q_{n-1}\\
| & | &  & |\\
\end{matrix}\right)\end{equation}

donde

\begin{align}
  q_0 = {} &\frac{a_0}{||a_0||} \\
  q_1 ={} & \frac{a_1^{\perp}}{||a_1^{\perp}||} \,\,\,\,\,\,\,\,\,\,\,\,\,\,   a_1^{\perp} = {} a_1 - \left \langle a_1, q_0  \right \rangle q_0 \\
  q_2 ={} & \frac{a_2^{\perp}}{||a_2^{\perp}||} \,\,\,\,\,\,\,\,\,\,\,\,\,\,   a_2^{\perp} = {} a_2 - \left \langle a_2, q_0  \right \rangle q_0 - \left \langle a_2, q_1  \right \rangle q_1\\
  & \vdots \\
  q_j ={} & \frac{a_j^{\perp}}{||a_j^{\perp}||} \,\,\,\,\,\,\,\,\,\,\,\,\,\,   a_j^{\perp} = {} a_{j} - \sum_{i = 0}^{j-1}\left \langle a_j, q_i  \right \rangle q_i
\end{align}



```python
#Teter en cuenta que este algoritmo solo encuetra valores propios para algunos casos
#El aogoritmo más general se escapa del alcance de este curso 
def calcula_Q(A):
  n = A.shape[0]
  Q = np.zeros(A.shape)
  Q[:, 0] = A[:, 0] / norm(A[:, 0])
  for j in range(1,n):
    Q[:,j] = calcula_qj(Q,A[:,j],j,n)
  return Q

def calcula_qj(Q,aj,j,n):
  suma = np.zeros(n)
  for i in range(j):
    suma += (aj @ Q[:, i])*Q[:, i]
  aj_p = aj - suma
  return aj_p / norm(aj_p)

def norm(v):
  return np.sqrt(v @ v)

def algoritmo_QR(A,iteraciones):
  Ak = A.copy().astype(float)
  for k in range(iteraciones):
    Qk = calcula_Q(Ak)
    Ak = Qk.T @ Ak @ Qk
  return np.diag(Ak)

```


```python
A = np.array([
              [ 1, 2, 3],
              [ 1, 2, 1],
              [ 3, 2, 1]
])

vals_cal = algoritmo_QR(A,15)
vals = np.linalg.eigvals(A)

print('valores propios con nuestro método:')
print(np.round(vals_cal,10),'\n')
print('valores propios con numpy:')
print(np.round(vals,10))

```

    valores propios con nuestro método:
    [ 5.23606798 -2.          0.76393202] 
    
    valores propios con numpy:
    [ 5.23606798 -2.          0.76393202]


**Referencias:**

**Universidad de los Andes**

**Profesor: Diego Alberto Castro Rodríguez**


```python

```

b4d98299e967c0e227139ba1b5dd50a09ce67532
ipynb
Jupyter Notebook
Notebooks/09 - inversion_matricial_y_calculo_de_autovalores.ipynb
Computational-Science-Uniandes/MetodosComputacionales1
bb74a6918265a54c414a3e2d2b2744b17cf0a4ac
MIT
2022-03-18T01:42:37.000Z
2022-03-18T01:42:37.000Z
Notebooks/09 - inversion_matricial_y_calculo_de_autovalores.ipynb
Computational-Science-Uniandes/MetodosComputacionales1
bb74a6918265a54c414a3e2d2b2744b17cf0a4ac
MIT
Notebooks/09 - inversion_matricial_y_calculo_de_autovalores.ipynb
Computational-Science-Uniandes/MetodosComputacionales1
bb74a6918265a54c414a3e2d2b2744b17cf0a4ac
MIT
2022-01-26T17:26:29.000Z
2022-03-05T23:23:44.000Z
Qwen/Qwen-72B
1. YES
2. YES