# Eigenvalue and eigenvectors calculation

$$
A\mathbf{x} = \lambda \mathbf{x}
$$

### Power method (vector iteration)
- find the largest eigenvalue $\lambda_{max}$
\begin{align}
\mathbf{q}_k & = \frac{\mathbf{z}_{k-1}}{\|\mathbf{z}_{k-1}\|_2}\\
\mathbf{z}_k & = A\mathbf{q}_{k}\\
\lambda_{max}^k & = \mathbf{q}^T_k \mathbf{z}_k
\end{align}


```python
%matplotlib inline
from numpy import *
from matplotlib.pyplot import *
import numpy.linalg
import scipy.linalg
```


```python
n = 9
h = 1./(n-1)
x=linspace(0,1,n)

# construct matrix
a = -ones((n-1,))
b = 2*ones((n,))
A = (diag(a, -1) + diag(b, 0) + diag(a, +1))
A /= h**2

#print A
```


```python
z0 = ones_like(x) # starting 

def PM(A, z0, tol=1e-5, nmax=500):
    q = z0/numpy.linalg.norm(z0,2)
    it = 0
    err = tol + 1.
    while (it < nmax and err > tol):
        z = dot(A,q)
        l = dot(q.T,z)
        err = numpy.linalg.norm(z-l*q,2)
        q = z/numpy.linalg.norm(z,2)
        
        it += 1
    print("error =", err, "iterations =", it)
    print("lambda_max =", l)
    return l,q

l,x = PM(A,z0)   
l_np, x_np = numpy.linalg.eig(A)

print("numpy")
print(l_np)
```

    error = 8.456086478475517e-06 iterations = 82
    lambda_max = 249.73523408577807
    numpy
    [249.73523409 231.55417528 203.23651229 167.55417528 128.
       6.26476591  24.44582472  88.44582472  52.76348771]


### Inverse power method
- find the eigenvalue $\lambda$ **closest** to $\mu$
\begin{align}
M & = A-\mu I\\
M & = LU \\
& \\
M\mathbf{x}_k &= \mathbf q_{k-1}\\
\mathbf{q}_k & = \frac{\mathbf{x}_k}{\|\mathbf{x}_k\|_2}\\
\mathbf{z}_k & = A\mathbf{q}_{k}\\
\lambda^k & = \mathbf{q}^T_k \mathbf{z}_k
\end{align}



```python
def IPM(A, x0, mu, tol=1e-5, nmax=500):
    M = A - mu*eye(len(A))
    P,L,U = scipy.linalg.lu(M)
    err = tol + 1.
    it = 0
    q = x0/numpy.linalg.norm(x0,2)
    while (it < nmax and err > tol):
        y = scipy.linalg.solve_triangular(L, dot(P.T,q), lower=True)
        x = scipy.linalg.solve_triangular(U, y)
        q = x/numpy.linalg.norm(x,2)
        z = dot(A,q)
        l = dot(q.T,z)
        err = numpy.linalg.norm(z-l*q,2)
        it += 1
    print("error =", err, "iterations =", it)
    print("lambda =", l)
    return l,q


l,x = IPM(A,z0,6.)
```

    error = 2.6310164587523326e-06 iterations = 3
    lambda = 6.2647659142204954


3bfb88ccb4493315028021ee13c22833ac48d7d3
ipynb
Jupyter Notebook
notebooks/GR_lab05_eigenvalues.ipynb
mapenzo-ph/numerical-analysis-2021-2022
952808d7fa7a9e718274592104be24882acef3ec
CC-BY-4.0
notebooks/GR_lab05_eigenvalues.ipynb
mapenzo-ph/numerical-analysis-2021-2022
952808d7fa7a9e718274592104be24882acef3ec
CC-BY-4.0
notebooks/GR_lab05_eigenvalues.ipynb
mapenzo-ph/numerical-analysis-2021-2022
952808d7fa7a9e718274592104be24882acef3ec
CC-BY-4.0
Qwen/Qwen-72B
1. YES
2. YES