# 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 (como lo enuncia el  **[teorema espectral](https://en.wikipedia.org/wiki/Spectral_theorem)**) 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 sin solución única  (i.e. en el que la matriz asociada 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 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 no cambia abruptamente hacia cero o infinito, 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}
\begin{align}
\mathbf{x}^{(\text{inicial})} \neq \mathbf{0} \quad & \text{ si } t=1 \\
& \\
\frac{1}{\Vert \mathbf{A}\mathbf{x}^{(t-1)} \Vert} \mathbf{A}\mathbf{x}^{(t-1)} \quad & \text{ si } t >1
\end{align}
\end{cases}
$$

Además de la sucesión vectorial, definimos una sucesión de reales $c_t$ como:

$$
c_{t} =  (\mathbf{x}^{(t-1)})^{T} \mathbf{A} \mathbf{x}^{(t-1)} = \mathbf{x}^{(t-1)} \cdot \mathbf{A} \mathbf{x}^{(t-1)} 
$$

Dadas ambas sucesiones,$\{c_t\}$ y $\{\mathbf{x}^{(t)}\}$, 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}$ y $\mathbf{v}$ su vector propio asociado. 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 que su función regrese un arreglo tal que el primer elemento sea la sucesión de las $c_t$ y el segundo que sea la de $\mathbf{x}^{(t)}$

### 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}
$$

Valores propios:
 
$$
(\lambda_1,\lambda_2)= (3,5)
$$

Vectores propios **normalizados**:

\begin{gather}
\mathbf{v}_1 =  (0,1) \\
\mathbf{v}_2 =  (1,0)
\end{gather}


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

Valores propios: 

$$
(\lambda_1,\lambda_2,\lambda_3)= (-3.18788,-0.88679, 7.07467)
$$

Vectores propios **normalizados**:

\begin{gather}
\mathbf{v}_1 =  (-0.255232 , -0.601302 ,  0.757161) \\
\mathbf{v}_2 =  ( 0.824038 , -0.544925 , -0.154979) \\
\mathbf{v}_3 =  (-0.505785 , -0.584374 , -0.634577)
\end{gather}

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

Valores propios:

$$
(\lambda_1,\lambda_2,\lambda_3,\lambda_4) = (-4.70726, -1.79468,  3.34440, 14.15754 )
$$

Vectores propios **normalizados**:

\begin{gather}
\mathbf{v}_1 = (0.319972 , -0.335648 , -0.536533 ,  0.705046) \\
\mathbf{v}_2 = ( 0.452603 ,  0.763066 , -0.429371 , -0.168882) \\
\mathbf{v}_3 = ( 0.802188 , -0.401573 ,  0.269615 , -0.350058) \\
\mathbf{v}_4 = (-0.221952 , -0.379222 , -0.674597 , -0.593167)
\end{gather}

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

**Sugerencia 1** si lo vas a hacer por primera vez, toma un número de iteraciones $T$ pequeño ($T \sim 10$) para que veas que tu función `potencias` si funciona bien y tengas una idea del tiempo que tarda en correr

**Sugerencia 2**: recuerden si $\mathbf{v}$ es vector propio, entonces, para todo $k \neq 0$, $k \mathbf{v}$ también es vector propio. El vector propio que obtengan como resultado del método de potencias, en esta versión **va a estar normalizado pero puede tener signo contrario a los mostrados como soluciones correctas**

## Otros valores propios: "desinflación" (deflation)

El método de potencias previamente mostrado solo nos genera el valor propio más grande en valor absoluto. Si quisiéramos extenderlo para obtener otros valores propios, necesitaríamos construir otra matriz en la cual el valor propio más grande en valor absoluto sea remplazado por un cero, pero todos los demás valores propio sean los mismos.

### Ejercicio 3

Sea $\mathbf{A}$ una matriz real y simétrica de $n \times n$ con valores y vectores propios $\lambda_1, \ldots, \lambda_n$ y $\mathbf{v}_1, \ldots, \mathbf{v}_n$ tales que $ |\lambda_1| \leq |\lambda_2| \leq  \ldots \leq |\lambda_n|$.

Utiliza el **[teorema espectral](https://en.wikipedia.org/wiki/Spectral_theorem)** para demostrar que la matriz $\mathbf{B}_k (\mathbf{A})$, definida como

$$
\mathbf{B}_k (\mathbf{A}) = \mathbf{A} - \frac{\lambda_k}{\Vert \mathbf{v}_k \Vert^2}\mathbf{v}_k^T \mathbf{v}_k
$$

Tiene valores propios $0,\lambda_1,\ldots,\lambda_{k-1},\lambda_{k+1},\ldots,\lambda_n$ y sus correspondientes vectores propios asociados $\mathbf{v}_k,\mathbf{v}_1,\ldots,\mathbf{v}_{k-1}, \mathbf{v}_{k+1},\ldots,\mathbf{v}_n$. Es decir, la matriz $\mathbf{B}_k$ tiene los mismos valores propios y vectores propios de $\mathbf{A}$ pero $\lambda_k$ se ha convertido en $0$

**Sugerencia 1** recuerda que
$$
\mathbf{v}^T \mathbf{w} = \begin{pmatrix}
v_1 w_1 & v_1 w_2 & \ldots & v_1 w_n \\
v_2 w_1 & v_2 w_2 & \ldots & v_2 w_n \\
\vdots & \vdots & \ddots & \vdots \\
v_n w_1 & v_n w_2 & \ldots & v_n w_n \\
\end{pmatrix}
$$

Y que  $(\mathbf{v}^T \mathbf{w}) \mathbf{u} = (\mathbf{w} \cdot \mathbf{u}) \mathbf{v}^T  $ 

**Sugerencia 2**: el teorema espectral dice que una matriz real y simétrica tiene una base ortonormal de vectores propios.

Así, podemos proponer el siguiente procedimiento para encontrar todos los valores y vectores propios: utilizamos el método de potencias sobre $\mathbf{A}$ para encontrar el valor propio más grande en magnitud  $\lambda_n$  y el vector propio asociado $\mathbf{v}_n$. 

Desupués, volvemos a aplicar el método de potencias pero ahora a $\mathbf{B}_n(\mathbf{A})$, el cuál nos regresara el valor propio más grande en valor absoluto(y su vector propio asociado ) de $\mathbf{B}_n(\mathbf{A})$: $\lambda_{n-1}$, segundo valor propio más grande de $\mathbf{A}$. 

Posteriormente, si ahora le volvemos a aplicar el método de potencias a $\mathbf{B}_n(\mathbf{B}_n (\mathbf{A})) $, obtendremos así $\lambda_{n-2}$ y su vector propio asociado. De forma general, si aplicamos el método de potencias a

$$
C_i = \underbrace{\mathbf{B}_n \circ \ldots \circ \mathbf{B}_{n}}_{\text{i veces}}(\mathbf{A})
$$

Obtendremos el valor propio $\lambda_{n-i}$ de $\mathbf{A}$. A esta técnica se le conoce como la *desinflación* (deflation, en inglés) de la matriz $\mathbf{A}$

### Ejercicio 4

Implementa una función `matrizB(A,w,s)` con $\mathbf{A}$ una matriz simétrica de $n \times n$ y $\mathbf{w}$ un vector de longitud $n$ y $s$ un número y te regrese la matriz $\mathbf{B}$ definida como

$$
\mathbf{B} = \mathbf{A} - \frac{s}{\Vert\mathbf{w}\Vert^2} \mathbf{w}^T \mathbf{w}
$$

### Ejercicio 5

Implementa una función `valoresPropios(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$. 

Tu función debe de regresar dos arreglos: uno con todos los valores propios $\lambda_i$ de la matriz $\mathbf{A}$ y otro con los vectores propios asociados $\mathbf{v_i}$ aplicando el método de potencias con $T$ iteraciones y vector inicial $x_{inicial}$ con la técnica de desinflación. Usa como base las funciones `matrizB` y `potencias`.

### Ejercicio 6

Aplica tu función `valoresPropios` a las matrices del ejercicio 2 para obtener sus valores y vectores propios. Compara con los resultados analíticos ahí mostrados para garantizar que tu método funciona bien.

**Sugerencia**: recuerda que el método de potencias regresa los valores propios normalizados pero no necesariamente con la misma orientación, por lo que la comparación con los resultados analíticos tiene que hacerse tomando esto en cuenta.

### Ejercicio 7

Obtén la matriz diagonal $\mathbf{D}$ para las matrices $\mathbf{A}$ del ejercicio 2. tal que

$\mathbf{A} = \mathbf{PDP^{-1}}$

**Sugerencia** Recuerda que las columnas de $\mathbf{P}$ son los vectores propios normalizados.

## La forma fácil: 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 (no necesariamente simétrica) `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 8

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 el ejercicio 2 y en el ejercicio 6

## 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 9

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 10

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 11 (solo intentar si ya se curso 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.

7eafa3a7979cbe9270ef89fc9a2c53dc77e86ee2
ipynb
Jupyter Notebook
files/fiscomp_2021-1/material/clase13.ipynb
sayeg84/sayeg84.github.io
18f2e36dd7252603fad8f7093dc5aa00fc721be4
MIT
files/fiscomp_2021-1/material/clase13.ipynb
sayeg84/sayeg84.github.io
18f2e36dd7252603fad8f7093dc5aa00fc721be4
MIT
files/fiscomp_2021-1/material/clase13.ipynb
sayeg84/sayeg84.github.io
18f2e36dd7252603fad8f7093dc5aa00fc721be4
MIT
Qwen/Qwen-72B
1. YES
2. YES