# Métodos para valores y vectores propios

Uno de los problemas matemáticos que surgen más frecuentemente en física es, dada una matriz $\mathbf{A}$ de $n \times n$, encontrar sus valores y vectores propios (comúnmente mal llamados *eigenvalores* y *eigenvectores*). Un vector propio $\mathbf{v} \neq \mathbf{0}$ y su valor propio asociado $\lambda$ cumplen la propiedad de

$$
\mathbf{Av} = \lambda \mathbf{v}
$$

Los valores y vectores propios normalmente se obtienen para después generar una matriz diagonal $\mathbf{D}$ y una matriz ortogonal $\mathbf{P}$ (cuyas columnas son los vectores propios) que cumpla que $\mathbf{A} = \mathbf{PDP^{-1}}$. Es decir, estan relacionados con la **diagonalización** de una matriz

No todas las matrices son diagonalizables, es decir, no se puede construir una base de $n$ vectores propios de la matriz. Una condición suficiente y necesaria para que esto se pueda es la matriz $\mathbf{A}$ sea **simétrica**. Por simplicidad, supondremos que solo trabajamos con matrices simétricas

Obtener los valores y vectores propios de una matriz computacionalmente es un problema **muy difícil** pues involucra resolver un sistema de ecuaciones lineales en el que la matriz no es invertible. La mayoría de los métodos para encontrar todos los valores y vectores propios escapan al alcance de este curso debido a su complejidad.

Sin embargo, existe un método sencillo que nos permite obtener el valor propio más grande en valor absoluto y su vector propio asociado. Dicho método se puede extender (de manera no trivial) para obtener todos los valores y vectores propios dada una aproximación inicial suficientemente buena de ellos

## Método de potencias

Notemos que, la propiedad que define a los valores y vectores propios:

$$
\mathbf{Av} = \lambda \mathbf{v}
$$

Es muy similar a la propiedad que cumplen los puntos fijos $x^*$ de una función real $F(x)$

$$
F(x^*) = x^*
$$

Así, esto nos podría sugerir que, análogo a las iteraciones de punto fijo, la siguiente sucesión vectorial convergería a un vector propio:

$$
\mathbf{x}^{(t)} = \begin{cases}
\mathbf{x}^{(\text{inicial})} \neq \mathbf{0} & \text{si } t=1 \\
\mathbf{A}\mathbf{x}^{(t-1)} & \text{si } t >1
\end{cases}
$$

Sin embargo, esa sucesión puede **generar problemas**. Notemos que, si el tomamos como condición incial $\mathbf{x}^{(\text{inicial})} = \mathbf{v}$ un vector propio con valor propio asociado $\lambda$, entonces

$$
\mathbf{x}^{(t)} = \mathbf{A}^{(t-1)} \mathbf{v} = \lambda^{(t-1)} \mathbf{v}
$$

Para $\vert \lambda \vert > 1$, conforme $t \to \infty$ sucede que $\lambda^{(t-1)} \to \pm \infty$ y la sucesión vectorial divergirá. por otro lado, si $\vert \lambda \vert < 1$, conforme $t \to \infty$ sucede que $\lambda^{(t-1)} \to \pm 0$ y nuestro la sucesión converge a $\mathbf{0}$, que no puede ser vector propio por definición. Se puede probar que, aunque tomemos cualquier otro vector como condición inicial, eventualmente se cumplirá alguno de estos dos casos.

Una manera de resolver este problema es multiplicando el vector de la sucesión en cada paso de tal manera que su norma se haga más chica, es decir, haciendo un tipo de **normalización**. 

Sin embargo, dicha normalización se debe de hacer con cuidado. Así, proponemos la siguiente sucesión recursiva

$$
\mathbf{x}^{(t)} = \begin{cases}
\mathbf{x}^{(\text{inicial})} \neq \mathbf{0} & \text{si } t=1 \\
\frac{1}{c_t} \mathbf{A}\mathbf{x}^{(t-1)} & \text{si } t >1
\end{cases}
$$

con 

$c_t = \max_{1 \leq i \leq n}{ \{ |x_{i}^{(t-1)}|\}} = \max{\{|x_1^{(t-1)}|,|x_2^{(t-1)}|,\ldots, |x_n^{(t-1)}|\}}$

Dada esta sucesión, se puede garantizar que

$$
\lim_{t \to \infty} c_t = \lambda \quad \quad \lim_{t \to \infty} \mathbf{x}^{(t)} = \mathbf{v}
$$

Con $\lambda$ el valor propio más grande (en valor absoluto) de la matriz $\mathbf{A}$. A la construcción de dicha sucesión se le llama el **método de potencias** para valores y vectores propios,

### Ejercicio 1

Implementa una función `potencias(A,T,x_inicial)` con $\mathbf{A}$ una matriz simétrica de $n \times n$, $T$ un número natural y $x_{inicial}$ es un vector de longitud $n$. La función debe de construir los primeros $T$ términos de las sucesiones $c_t$ y $\mathbf{x}^{(t)}$ y debe de regresarlos.

**Sugerencia 1**: para regresar ambas sucesiones, puedes hacer `return [sucC,sucX]` para regresar un arreglo tal que el primer elemento sea la sucesión de las $c_t$ y el segundo que sea la de $x^{(t)}$

**Sugerencia 2**: la función `maximum(arr)` te regresa el elemento máximo del arreglo `arr` (siempr ey cuando `arr` es un arreglo que contenga números enteros o flontantes como elementos)

### Ejercicio 2

Utiliza la función `potencias` para encontrar el valor propio más grande en valor absoluto y su vector propio asociado para las siguientes matrices:

1. $$
\begin{pmatrix}
5 &  0 \\
0 &  3  \\
\end{pmatrix}
$$

Respuesta: $\lambda = 5$, $\mathbf{v} = (1,0)$

2. $$
\begin{pmatrix}
1 &  2 & 3 \\
2 &  1 & 4 \\
3 &  4 & 1 \\
\end{pmatrix}
$$

Respuesta: $\lambda =  7.07467$, $\mathbf{v} = (-0.50578,-0.58437,-0.63457)$

3. $$
\begin{pmatrix}
2 & 0 & 4 & 0 \\
0 & 1 & 3 & 5 \\
4 & 3 & 5 & 7 \\
0 & 5 & 7 & 3
\end{pmatrix}
$$

respuesta: $\lambda = 14.157537$, $\mathbf{v} = (-0.221951,-0.379221,-0.674596,-0.593167)$

Realiza un análisis del error absoluto como función del número de iteraciones.

## Funciones `eigvals` y `eigvecs` de la librería `LinearAlgebra`

En la librería `LinearAlgebra` de Julia, existen dos funciones `eigvals(A)` y `eigvecs(A)` que nos regresan los valores y los vectores propios, respectivamente, de una matriz cuadrada `A`. 

Estas funciones utilizan algoritmos muy sofisticados y rápidos para calcular los vectores propios.


```julia
using LinearAlgebra
```


```julia
# defino una matriz
A = [1.0 1.0 -1.0; 1.0 -1.0 1.0 ; -1.0 1.0 1.0]
display(A)
```


    3×3 Array{Float64,2}:
      1.0   1.0  -1.0
      1.0  -1.0   1.0
     -1.0   1.0   1.0



```julia
# obtengo sus valores propios
valPropios = eigvals(A)
display(valPropios)
```


    3-element Array{Float64,1}:
     -1.9999999999999996
      0.9999999999999997
      2.000000000000001



```julia
# obtengo los vectores propios
vecPropios = eigvecs(A)
# vecPropios será una matriz donde cada columna es un vector propio
# la columna i corresponde al vector propio asociado al iésimo elemento de `eigvals(A)`
display(vecPropios)
```


    3×3 Array{Float64,2}:
      0.408248  0.57735  -0.707107
     -0.816497  0.57735  -1.11022e-16
      0.408248  0.57735   0.707107



```julia
# checando que si sean vectores propios
# iteramos sobre el número de vectores propios
for i in 1:3
    # multiplico A por cada vector propio
    res = A*vecPropios[:,i]
    print("El vector propio $i es: ")
    println(res)
    print("El valor propio $i es: ")
    println(valPropios[i])
    println()
    # itero sobre cada coordenada para dividirla y verla
    for j in 1:3
        println("La entrada $j del vector propio $i despues de multiplicarlo por A")
        println("dividida por la misma entrada antes de multiplicar es:")
        println(res[j]/vecPropios[j,i])
    end
    println()
    println()
end
```

    El vector propio 1 es: [-0.816496580927726, 1.6329931618554514, -0.816496580927726]
    El valor propio 1 es: -1.9999999999999996
    
    La entrada 1 del vector propio 1 despues de multiplicarlo por A
    dividida por la misma entrada antes de multiplicar es:
    -2.0000000000000018
    La entrada 2 del vector propio 1 despues de multiplicarlo por A
    dividida por la misma entrada antes de multiplicar es:
    -1.9999999999999991
    La entrada 3 del vector propio 1 despues de multiplicarlo por A
    dividida por la misma entrada antes de multiplicar es:
    -2.0000000000000018
    
    
    El vector propio 2 es: [0.5773502691896252, 0.5773502691896267, 0.5773502691896254]
    El valor propio 2 es: 0.9999999999999997
    
    La entrada 1 del vector propio 2 despues de multiplicarlo por A
    dividida por la misma entrada antes de multiplicar es:
    0.9999999999999987
    La entrada 2 del vector propio 2 despues de multiplicarlo por A
    dividida por la misma entrada antes de multiplicar es:
    1.0000000000000024
    La entrada 3 del vector propio 2 despues de multiplicarlo por A
    dividida por la misma entrada antes de multiplicar es:
    0.9999999999999989
    
    
    El vector propio 3 es: [-1.4142135623730954, 0.0, 1.4142135623730951]
    El valor propio 3 es: 2.000000000000001
    
    La entrada 1 del vector propio 3 despues de multiplicarlo por A
    dividida por la misma entrada antes de multiplicar es:
    2.0
    La entrada 2 del vector propio 3 despues de multiplicarlo por A
    dividida por la misma entrada antes de multiplicar es:
    -0.0
    La entrada 3 del vector propio 3 despues de multiplicarlo por A
    dividida por la misma entrada antes de multiplicar es:
    2.0
    
    


### Ejercicio 3

Utiliza las funciones `eigvals` y `eigvecs` para obtener todos los valores propios y vectores propios de las matrices del ejercicio 2. Verifica que tus resultados sean consistentes con lo obtenido en dicho ejercicio

## Aplicación interesante: osciladores armónicos acoplados

Considera el sistema de dos carros de masas $m_1$ y $m_2$, acoplados entre si y a dos paredes por tres resortes Hookeanos (que siguen la ley de Hooke) con constantes $k_1$, $k_2$ y $k_3$. El siguiente diagrama representa al sistema:



Al hacer el diagrama de fuerzas y plantear la segunda ley de Newton, resolviendo la ecuación diferencial usando una exponencial compleja, llegamos a que el desplazamiento respecto a la posición de equilibrio de cada carrito, denotado por $x_1(t)$ y $x_2(t)$ es de la forma

$$
\begin{align}
x_1(t) &= \alpha_1 \cos{(\omega t - \delta_1)} \\
x_2(t) &= \alpha_2 \cos{(\omega t - \delta_2)}
\end{align}
$$

Supongamos que los carros son de masas iguales $m_1 = m_2 = m$ y que $\delta_1 = \delta_2 = 0$. Así, la frecuencia $\omega$ y el vector de amplitudes $\mathbf{A} = [\alpha_1, \alpha_2]$ se puede encontrar de la siguiente manera:

$\mathbf{A}$ es un vector propio, y $m \omega^2$  su valor propio asociado, de la **matriz de constante de resortes** $\mathbf{K}$, definida como

$$
\mathbf{K} = \begin{pmatrix}
k_1 + k_2 & -k_2 \\
-k_2 & k_2 + k_3 \\
\end{pmatrix}
$$

Notemos que, por ser simétrica, la matriz $\mathbf{K} $ tendrá dos  vectores propios distintos, asociados a valores propios distintos. Cada uno de ellos define dos movimientos distintos de los carros acoplados que corresponden a sus **modos normales**.

### Ejercicio 4

Utilizando las funciones `eigvals` y `eigvecs`, escribe una función `carrosAcoplados2(k1,k2,k3,m)` que tome como argumento números flotantes correspondientes a las constantes de resorte $k_1$, $k_2$, $k_3$ y la masa $m$ y te regrese las dos frecuencias $\omega_1, \omega_2$ y sus vectores de amplitudes correspondientes $\mathbf{A}_1, \mathbf{A}_2$.

Revisa que para el caso $k_1=k_2=k_3=k$, se cumpla

$$
\omega_1 = \sqrt{\frac{k}{m}} \quad , \quad \mathbf{A}_1 = \frac{1}{\sqrt{2}} (1,1)
$$

$$
\omega_2 = \sqrt{\frac{3k}{m}} \quad , \quad \mathbf{A}_2 = \frac{1}{\sqrt{2}} (1,-1)
$$

**Sugerencia**: para comprobar que se cumplen esas relaciones, deberán darle valores arbitrarios a $k$ y $m$

### Ejercicio 5

Revisa la parte del tutorial de Plots sobre animaciones. Después de revisarlo y aprenderlo bien, dale valores arbitrarios a $l_1$, $l_2$ y $l_3$, las longitudes en equilibrio de los resortes 1, 2 y 3, respectivamente, y realiza dos animaciones del movimiento de los carritos, una para cada modo normal.

### Ejercicio 6 (solo para los que ya cursaron Mecánica Analítica)

Generaliza el sistema anterior para $n$ carritos de la misma masa acoplados mediante $n+1$ resortes con constantes $k_1,\ldots, k_{n+1}$ y escribe una función `carrosAcoplados(ks,m)` que toma como argumento un vector `ks` de longitud $n+1$, que contiene las constantes de los resortes, y la masa $m$ de los carritos y regresa las $n$ frecuencias $\omega_1, \ldots \omega_n$ y los $n$ vectores de amplitud $\mathbf{A}_1, \ldots, \mathbf{A}_n$ correspondientes a los modos normales del sistema.

###


```julia
A = [-4 2 -2 ; 2 -7 4 ; -2 4 -7]
for i in 1:3
    for j in i+1:3
        A[i,j] = A[j,i]
    end
end
display(A)
```


    3×3 Array{Int64,2}:
     -4   2  -2
      2  -7   4
     -2   4  -7



```julia
using LinearAlgebra
```


```julia
res = eigvecs(A)
```




    3×3 Array{Float64,2}:
      0.333333  0.894427  -0.298142
     -0.666667  0.447214   0.596285
      0.666667  0.0        0.745356




```julia
transpose(res)*A*res
```




    3×3 Array{Float64,2}:
     -12.0          0.0          -1.77636e-15
       0.0         -3.0           2.22045e-16
      -1.9984e-15   4.44089e-16  -3.0




```julia

```

d1017ce666b6493bc25ece0a63e731dff841478c
ipynb
Jupyter Notebook
files/fiscomp_2020-4/material/clase13.ipynb
sayeg84/sayeg84.github.io
18f2e36dd7252603fad8f7093dc5aa00fc721be4
MIT
files/fiscomp_2020-4/material/clase13.ipynb
sayeg84/sayeg84.github.io
18f2e36dd7252603fad8f7093dc5aa00fc721be4
MIT
files/fiscomp_2020-4/material/clase13.ipynb
sayeg84/sayeg84.github.io
18f2e36dd7252603fad8f7093dc5aa00fc721be4
MIT
Qwen/Qwen-72B
1. YES
2. YES