# NUME

Instructor: Keny **ORDAZ**
    
# Systems of Linear Algebraic Equations: Iterative Methods

$\newcommand{\vb}[1]{\boldsymbol{#1}}\newcommand{\RR}[1]{\mathbb{R}^{#1}}
\newcommand{\vt}[1]{\vb{#1}^\top}$

## Introduction


Iterative methods vs Direct methods of solution.

$$\vb{A} \vb{x} = \vb{b}$$

| Direct methods | Iterative (indirect) methods |
|:---------------|:------------------|
| Finite number of operations | Unknown number of iterations, slower |
| Best for dense and small $\vb{A}$ | Best for sparse and large $\vb{A}$|
| Always provide a solution | Not always converge to the solution |


- Iterative methods start with an initial guess of the solution $\vb{x}$ and then repeatedly improve the solution until the change in $\vb{x}$ becomes negligible.

- Because the required number of iterations can be large, the indirect methods are, in general, slower than their direct counterparts.

- Iterative methods do have the following two advantages that make them attractive for certain problems:

    1. It is feasible to store only the nonzero elements of the coefficient matrix. This makes it possible to deal with very large matrices that are sparse, but not necessarily banded. In many problems, there is no need to store the coefficient matrix at all.

    2. Iterative procedures are self-correcting, meaning that roundoff errors (or even arithmetic mistakes) in one iterative cycle are corrected in subsequent cycles.

- It can be shown that _convergence_ is **guaranteed only if the coefficient matrix is diagonally dominant**. The initial guess for $\vb{x}$ plays no role in determining whether convergence takes place: if the procedure converges for one starting vector, it would do so for any starting vector. The initial guess affects only the number of iterations that are required for convergence.

$\newcommand{\x}{\vb{x}}
\newcommand{\xT}{\vt{x}}
\newcommand{\b}{\vb{b}}
\newcommand{\A}{\vb{A}}
\newcommand{\bT}{\vt{b}}
\newcommand{\AT}{\vt{A}}$


```python
import numpy as np
import scipy
import sympy as sp
import matplotlib.pyplot as plt
```


```python
!ls book/*.py
```

    book/LUdecomp.py       book/example8_8.py     book/newtonRaphson2.py
    book/LUdecomp3.py      book/gaussElimin.py    book/plotPoly.py
    book/LUdecomp5.py      book/gaussNodes.py     book/polyFit.py
    book/LUpivot.py        book/gaussPivot.py     book/polyRoots.py
    book/bisection.py      book/gaussQuad.py      book/powell.py
    book/bulStoer.py       book/gaussQuad2.py     book/printSoln.py
    book/choleski.py       book/gaussSeidel.py    book/rational.py
    book/conjGrad.py       book/gerschgorin.py    book/ridder.py
    book/cubicSpline.py    book/goldSearch.py     book/romberg.py
    book/downhill.py       book/householder.py    book/rootsearch.py
    book/eigenvals3.py     book/inversePower.py   book/run_kut4.py
    book/error.py          book/inversePower3.py  book/run_kut5.py
    book/euler.py          book/inversePower5.py  book/sortJacobi.py
    book/evalPoly.py       book/jacobi.py         book/stdForm.py
    book/example8_1.py     book/lamRange.py       book/sturmSeq.py
    book/example8_3.py     book/linInterp.py      book/swap.py
    book/example8_4.py     book/midpoint.py       book/trapezoid.py
    book/example8_5.py     book/neville.py        book/triangleQuad.py
    book/example8_6.py     book/newtonPoly.py
    book/example8_7.py     book/newtonRaphson.py


$\A \x = \lambda \x$


```python
def build_householder_projector(v):
    #P = I - 2 * v v^T / v^T v
    return np.eye(len(v)) - 2 * np.outer(v, v) / np.dot(v, v)

v = np.array([1, 2, 3, 4, 5, 6])

P = build_householder_projector(v)
P
```




    array([[ 0.97802198, -0.04395604, -0.06593407, -0.08791209, -0.10989011,
            -0.13186813],
           [-0.04395604,  0.91208791, -0.13186813, -0.17582418, -0.21978022,
            -0.26373626],
           [-0.06593407, -0.13186813,  0.8021978 , -0.26373626, -0.32967033,
            -0.3956044 ],
           [-0.08791209, -0.17582418, -0.26373626,  0.64835165, -0.43956044,
            -0.52747253],
           [-0.10989011, -0.21978022, -0.32967033, -0.43956044,  0.45054945,
            -0.65934066],
           [-0.13186813, -0.26373626, -0.3956044 , -0.52747253, -0.65934066,
             0.20879121]])




```python
u = np.array([4, 2, 3, 0, 0, 0])
w = np.array([1, 0, 0, 0, 0, 3])
y = np.array([ 0, 0, 0, 1, 0, 0])

yp = P @ y
wp = P @ w
up = P @ u
up, wp, yp
```




    (array([ 3.62637363,  1.25274725,  1.87912088, -1.49450549, -1.86813187,
            -2.24175824]),
     array([ 0.58241758, -0.83516484, -1.25274725, -1.67032967, -2.08791209,
             0.49450549]),
     array([-0.08791209, -0.17582418, -0.26373626,  0.64835165, -0.43956044,
            -0.52747253]))



$$\begin{aligned}P &= I - 2 \frac{ v v^T }{ v^T v}\end{aligned}$$

Aplicamos el proyector:

$$\begin{aligned}w' &= P w \\
 &= \left[I - 2 \frac{ v v^T }{ v^T v}\right] w\\
 &= w - \left[2 \frac{ v v^T }{ v^T v}\right] w\\
  &= w - \frac{2  }{ v^T v} v v^T w \\
  &= w - 2\frac{v^T w  }{ v^T v} v  \\
  &= w - 2\frac{v \cdot w  }{ v \cdot v} v \\
  &= w - 2\frac{v \cdot w  }{ \|v\| \|v\|} v \\
  &= w - 2(\hat{v} \cdot w)  \hat{v}
\end{aligned}$$


donde:

- $w'$ es la imagen de $w$ bajo la proyección (transformación) $P$
- $\hat{v} = \frac{v}{ \|v\|}$ es el vector unitario (normalizado de $v$).
- $v \cdot v = \| v \|^2$



```python
def hh_projector(w, v):
    # P = I - 2 * v v^T / v^T v
    # w' = w - 2 * (np.dot(v, w)) * v
    return w - (2 * np.dot(v, w) / np.dot(v, v)) * v 
```


```python
def hh_projector_normalized(w, v):
    # P = I - 2 * v v^T / v^T v
    # w' = w - 2 * (np.dot(v, w)) * v
    return w - 2 * (np.dot(v, w)) * v

"""
w - v ==> n ops
w . v ==> n-1 sums + n prods = 2 n - 1
2 op
a . v ==> n ops

2(2n - 1) + 2 + n + n = 4 n - 2 + 2 + 2n = 6n

n = 3 --> fu
"""
None
```


```python
ui = hh_projector(u, v)
wi = hh_projector(w, v)
yi = hh_projector(y, v)

# n = 3
# Forma matricial P @ w --> n*m sumas + n*m prods => n^2 + n^2 = 2 * n^2 ops = 2 * (3)^2 = 2 * 9 = 18 ops
# Forma funcional (operaciones vectoriales): 6 n; con n=3 --> 6(3) = 18 ops

# n = 6
# Forma matricial P @ w --> 2 * n^2 ops = 2 * (6)^2 = 2 * 36 = 72 ops
# Forma funcional (operaciones vectoriales): 6 n; con n=6 --> 6(6) = 36 ops

## O(n) vs O(n^2), 6 n  vs 2 n^2


ui, up, wi, wp, yi, yp
ui - up, wi - wp, yi - yp
```




    (array([ 4.44089210e-16,  0.00000000e+00, -4.44089210e-16,  2.22044605e-16,
             0.00000000e+00,  0.00000000e+00]),
     array([ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
             0.00000000e+00, -1.11022302e-16]),
     array([ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
            -5.55111512e-17,  0.00000000e+00]))




```python
for n in [2, 3, 6, 20, 50, 200]:
    print(f'{n:5} {6*n:8} {2*n**2:10}')
```

        2       12          8
        3       18         18
        6       36         72
       20      120        800
       50      300       5000
      200     1200      80000


Con $\vb{v}$ normalizado, i.e. $\hat{\vb{v}}$:


```python
vnorm = v / np.linalg.norm(v)

projec_v = lambda w : w - 2 * (np.dot(vnorm, w)) * vnorm
```


```python
uin = projec_v(u)
win = projec_v(w)
yin = projec_v(y)
uin - up, win - wp, yin - yp, uin == ui, win == wi, yin == yi, np.allclose(win, wi)
```




    (array([ 4.44089210e-16,  0.00000000e+00, -4.44089210e-16,  2.22044605e-16,
            -2.22044605e-16,  0.00000000e+00]),
     array([ 1.11022302e-16,  1.11022302e-16,  0.00000000e+00,  2.22044605e-16,
             0.00000000e+00, -1.11022302e-16]),
     array([ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
            -5.55111512e-17,  0.00000000e+00]),
     array([ True,  True,  True,  True, False,  True]),
     array([False, False,  True, False,  True,  True]),
     array([ True,  True,  True,  True,  True,  True]),
     True)



# Sparse Matrix


```python
from scipy.sparse import csr_matrix

# create a 2-D representation of the matrix
A = np.array([[1, 0, 0, 0, 0, 0], [0, 0, 2, 0, 0, 1],\
 [0, 0, 0, 2, 0, 0]])
print("Dense matrix representation: \n", A)

# convert to sparse matrix representation 
S = csr_matrix(A)
print("Sparse matrix: \n",S)

# convert back to 2-D representation of the matrix
B = S.todense()
print("Dense matrix: \n", B)
```

    Dense matrix representation: 
     [[1 0 0 0 0 0]
     [0 0 2 0 0 1]
     [0 0 0 2 0 0]]
    Sparse matrix: 
       (0, 0)	1
      (1, 2)	2
      (1, 5)	1
      (2, 3)	2
    Dense matrix: 
     [[1 0 0 0 0 0]
     [0 0 2 0 0 1]
     [0 0 0 2 0 0]]



```python

```


```python

```


```python
A = np.zeros((3, 6))
A[0, 0] = A[1, 5] = 1
A[1, 2] = A[2, 3] = 2
A
#A = np.array([[1 ,0 , 0, 0 , 0....]])

x = np.ones(6)
x[-1] = -5
# GAXPY  y = y + A x
A, x.reshape((-1,1)), A @ x # 3 * 6 prod y 3 * 5 sum = 33 ops =  m * (2n - 1) = 2 m*n - m
```




    (array([[1., 0., 0., 0., 0., 0.],
            [0., 0., 2., 0., 0., 1.],
            [0., 0., 0., 2., 0., 0.]]),
     array([[ 1.],
            [ 1.],
            [ 1.],
            [ 1.],
            [ 1.],
            [-5.]]),
     array([ 1., -3.,  2.]))




```python
# b0 = a00 * x0
# b1 = a12 * x2 + a15 * x5
# b2 = a23 * x3 = b2

# SPARSE: prods: 4, sum: 1 --> 5 ops
# DENSE: 2 * 3 * 6 = 36 ops

def function_example_sparse(x):
    # 5 ops vs  2mn - m = 36 - 6 = 30 ops matricial ? (6 p + 5s)*3 = 33
    assert len(x) == 6
    b = np.zeros(3)
    b[0] = S[0, 0] * x[0]
    b[1] = S[1, 2] * x[2] + S[1, 5] * x[5]
    b[2] = S[2, 3] * x[3]
    return b
    

```


```python
b_i = function_example_sparse(x)
b_i, b_i - A @ x
```




    (array([ 1., -3.,  2.]), array([0., 0., 0.]))



# Cases


```python
from scipy.sparse import coo_matrix, diags

a0 = np.array([[-2., 5, 9], [7, 1, 1], [-3, 7, -1]])
b0 = np.array([1, 6, -26.])

a1 = np.array([[12, -2., 3, 1],
               [-2, 15, 6, -3],
               [1, 6, 20, -4],
               [0, -3, 2, 9]])
b1 = np.array([0, 0, 20., 0])

a2 = np.array([[4, -1, 0, 0],
               [-1, 4, -1, 0],
               [0, -1, 4, -1],
               [0, 0, -1, 3]])
b2 = np.array([15.] + [10]*3)
w2 = 1.1
initial2 = b2 / np.diag(a2)

"""
a3 =
2, -1, 0
-1, 2, -1
0, -1, 1
b3 = 1, 1, 1

a4 =
3, 0, -1
0, 4, -2
-1, -2, 5
b4 = [4, 10, -10]
"""

row  = np.array([0, 5])
col  = np.array([0, 5])
data = np.array([-1., -1])
a5 = coo_matrix((data, (row, col)), shape=(6, 6)).toarray()
a5 += diags([1, 1, -2, 4, -2, 1, 1], [-5, -2, -1, 0, 1, 2,5], shape=(6, 6)).toarray()
b5 = [10, -8, 10, 10, -8, 10]


# ex 2.17

def sparse_exam_exer(a00, a01, a0n, bn, n):
    a = diags([a0n, a01, a00, a01, a0n], [-n+1, -1, 0, 1, n-1], shape=(n, n)).toarray()
    b = np.zeros(n)
    b[-1] = bn
    x = np.array([-n / 4 + i/2 for i in range(1, n + 1)])
    return a, b, x

n6 = 4
a6, b6, x6 = sparse_exam_exer(4, -1, 1, 100, n6)

ne17 = 20
ae17, be17, xe17 = sparse_exam_exer(2, -1, 1, 1, ne17)



"""
a1 = np.array([[]])
b2 = np.array([])

a1 = np.array([[]])
b2 = np.array([])

a1 = np.array([[]])
b2 = np.array([])
"""
initial2
a2, a0, a1, a5, a6, b6, ae17, be17, xe17, ae17 @ xe17 - be17
```




    (array([[ 4, -1,  0,  0],
            [-1,  4, -1,  0],
            [ 0, -1,  4, -1],
            [ 0,  0, -1,  3]]),
     array([[-2.,  5.,  9.],
            [ 7.,  1.,  1.],
            [-3.,  7., -1.]]),
     array([[12., -2.,  3.,  1.],
            [-2., 15.,  6., -3.],
            [ 1.,  6., 20., -4.],
            [ 0., -3.,  2.,  9.]]),
     array([[ 3., -2.,  1.,  0.,  0.,  1.],
            [-2.,  4., -2.,  1.,  0.,  0.],
            [ 1., -2.,  4., -2.,  1.,  0.],
            [ 0.,  1., -2.,  4., -2.,  1.],
            [ 0.,  0.,  1., -2.,  4., -2.],
            [ 1.,  0.,  0.,  1., -2.,  3.]]),
     array([[ 4., -1.,  0.,  1.],
            [-1.,  4., -1.,  0.],
            [ 0., -1.,  4., -1.],
            [ 1.,  0., -1.,  4.]]),
     array([  0.,   0.,   0., 100.]),
     array([[ 2., -1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
              0.,  0.,  0.,  0.,  0.,  0.,  1.],
            [-1.,  2., -1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
              0.,  0.,  0.,  0.,  0.,  0.,  0.],
            [ 0., -1.,  2., -1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
              0.,  0.,  0.,  0.,  0.,  0.,  0.],
            [ 0.,  0., -1.,  2., -1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
              0.,  0.,  0.,  0.,  0.,  0.,  0.],
            [ 0.,  0.,  0., -1.,  2., -1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
              0.,  0.,  0.,  0.,  0.,  0.,  0.],
            [ 0.,  0.,  0.,  0., -1.,  2., -1.,  0.,  0.,  0.,  0.,  0.,  0.,
              0.,  0.,  0.,  0.,  0.,  0.,  0.],
            [ 0.,  0.,  0.,  0.,  0., -1.,  2., -1.,  0.,  0.,  0.,  0.,  0.,
              0.,  0.,  0.,  0.,  0.,  0.,  0.],
            [ 0.,  0.,  0.,  0.,  0.,  0., -1.,  2., -1.,  0.,  0.,  0.,  0.,
              0.,  0.,  0.,  0.,  0.,  0.,  0.],
            [ 0.,  0.,  0.,  0.,  0.,  0.,  0., -1.,  2., -1.,  0.,  0.,  0.,
              0.,  0.,  0.,  0.,  0.,  0.,  0.],
            [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0., -1.,  2., -1.,  0.,  0.,
              0.,  0.,  0.,  0.,  0.,  0.,  0.],
            [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0., -1.,  2., -1.,  0.,
              0.,  0.,  0.,  0.,  0.,  0.,  0.],
            [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0., -1.,  2., -1.,
              0.,  0.,  0.,  0.,  0.,  0.,  0.],
            [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0., -1.,  2.,
             -1.,  0.,  0.,  0.,  0.,  0.,  0.],
            [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0., -1.,
              2., -1.,  0.,  0.,  0.,  0.,  0.],
            [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
             -1.,  2., -1.,  0.,  0.,  0.,  0.],
            [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
              0., -1.,  2., -1.,  0.,  0.,  0.],
            [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
              0.,  0., -1.,  2., -1.,  0.,  0.],
            [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
              0.,  0.,  0., -1.,  2., -1.,  0.],
            [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
              0.,  0.,  0.,  0., -1.,  2., -1.],
            [ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
              0.,  0.,  0.,  0.,  0., -1.,  2.]]),
     array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
            0., 0., 1.]),
     array([-4.5, -4. , -3.5, -3. , -2.5, -2. , -1.5, -1. , -0.5,  0. ,  0.5,
             1. ,  1.5,  2. ,  2.5,  3. ,  3.5,  4. ,  4.5,  5. ]),
     array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
            0., 0., 0.]))



### Matriz diagonal estrictamente dominante

- Por renglones

$${\displaystyle |a_{i,i}|>\sum _{j=1,j\neq i}^{n}|a_{i,j}|,\forall i\in \mathbb {N} :1\leq i\leq n}$$

- Por columnas

$${\displaystyle |a_{j,j}|>\sum _{i=1,i\neq j}^{n}|a_{i,j}|,\forall j\in \mathbb {N} :1\leq j\leq n}$$

# Resolución de sistemas lineales con métodos iterativos


```python
!cat book/jacobi.py
```

    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    



```python
!cat book/gaussSeidel.py
```

    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    



```python

```


```python
import numpy as np

def solve_sor(A, b, w, x):
    x = np.asarray(x)
    n = b.shape[0]
    xold = x.copy()
    convergence = False
    while not convergence:
        for i in range(n):
            s = 0.0
            for j in range(n):
                if j != i:
                    s += A[i,j]*x[j]
            suma = np.sum(A[i,:] @ x) - A[i,i]*x[i]
            #print(f'{suma=} \t\t {s=}\t{A[i,:] @ x}')
            x[i] = x[i] + w * ( (b[i] - s) / A[i, i] - x[i])
            # print(f'\t{x=}')
        # print(f'{x=}')
        convergence = np.allclose(A @ x, b, 1e-14) # or np.allclose(xold, x, 1e-14)
        xold = x.copy()
    return  x

def solve_jacobi_no(A, b, w, x):
    x = np.asarray(x)
    L, D, U = -np.tril(A, -1), np.diagonal(A), -np.triu(A, 1)
    assert (A == -L + np.diag(D) + -U).all()
    n = b.shape[0]
    xold = x.copy()
    convergence = False 
    k = 0
    while not convergence:
        k += 1
        x = ((L + U )@x - b)
        print(f'1 {x=}')
        x = x / D
        print(f'2 {x=}')
        convergence = k > 10 or np.allclose(A @ x, b, 1e-14) # or np.allclose(xold, x, 1e-14)
        xold = x.copy()
    return  x

def solve_jacobi(A, b, w, x, max_iter=1000):
    x = np.asarray(x)
    L, D, U = -np.tril(A, -1), np.diagonal(A), -np.triu(A, 1)
    assert (A == -L + np.diag(D) + -U).all()
    n = b.shape[0]
    xold = x.copy()
    convergence = False 
    k = 0
    while not convergence:
        k += 1
        for i in range(n):
            s = 0.0
            for j in range(n):
                if j != i:
                    s += A[i,j]*xold[j]
            suma = np.sum(A[i,:] @ x) - A[i,i]*x[i]
            #print(f'{suma=} \t\t {s=}\t{A[i,:] @ x}')
            x[i] = xold[i] + ( (b[i] - s) / A[i, i] - xold[i])
            # print(f'\t{x=}')
        # print(f'{x=}')
        convergence = k > max_iter or np.allclose(A @ x, b, 1e-14) # or np.allclose(xold, x, 1e-14)
        xold = x.copy()
    return  x


A = np.array([
    [4, -1, -6, 0.],
    [-5, -4, 10, 8],
    [0, 9, 4, -2],
    [1, 0, -7, 5]
    ])
b = np.array([2., 21, -12, -6]) 

E2_A = np.array([[10, -1, 2,0], [-1, 11, -1, 3], [2, -1, 10, -1], [3, -1, 8, 15]])
E2_b = np.array([6., 25, -11, 15])

for solver in [solve_sor]: #, solve_jacobi]:
    x = solver(A, b, 0.5, [0.0, 0, 0, 0])
    print(x, A@x - b)

y = solve_jacobi(E2_A, E2_b, 0. , np.zeros_like(E2_b))
print(f'sol: {y}\n{b - A @ y}')

```

    [ 3. -2.  2.  1.] [-6.99743374e-09  4.48017090e-09 -1.49637991e-09 -1.09573328e-09]
    sol: [ 0.98082192  1.88066971 -0.96377473  1.44322679]
    [ -5.82526636  31.51872146 -22.18447489 -20.943379  ]


# Conjugate Gradient Method

Find the vector $\vb{x}$ that minimizes the scalar function:
    
$$f(\vb{x}) = \frac{1}{2} \vt{x}\vb{A}\vb{x} - \vt{b}\vb{x}$$

where:

- $\vb{A}$ is symmetric ($\vb{A} = \vt{A}$) and positive definite ($\vt{x}\vb{A}\vb{x} > 0$)
- $f(\vb{x})$ is minimized when its gradient $\nabla f(\vb{x}) = \vb{A}\vb{x} - \vb{b} = \vb{0}$
- Minimization is equivalente to solving

$$\vb{A}\vb{x} = \vb{b}$$

## Conjugate Gradient Algorithm

- Choose $\vb{x}_0$, $\epsilon_\tau$ (error tolerance).
- $\vb{r}_0 \leftarrow \vb{b} - \A\x_0$
- $\vb{s}_0 \leftarrow \vb{r}_0$ (choose the direction of steepest descent), search direction vector
- For $i$ in $[0,1,2,\cdots]$ do:
    - $\alpha_k \leftarrow \frac{\vt{s}_k\vb{r}_k}{\vt{s}_k\A\vb{s}_k}$
    - $\x_{k+1} \leftarrow \x_k + \alpha_k \vb{s}_k$
    - $\vb{r}_{k+1} \leftarrow \vb{b} - \A\x_{k+1}$
    - if $\|\vb{r}_{k+1}\| \le \epsilon_\tau$ break
    - $\beta_k \leftarrow - \frac{\vt{r}_{k+1}\A\vb{s}_k}{\vt{s}_k\A\vb{s}_k}$
    - $\vb{s}_{k+1} \leftarrow \vb{r}_{k+1} + \beta_k \vb{s}_k$
    
    
c.f. Householder projector:

$$ p(w) = w - 2\frac{v^T w  }{ v^T v} v$$

Los vectores residuales son mutuamente ortogonales:

$$\vb{r}_i \cdot \vb{r}_j = 0, i\ne j$$


```python
# asignación: [=] a = b
# igualdad:    [==] a==b 

# a <-- b
# a = b
```


```python
!cat book/conjGrad.py
```

    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    



```python
from book.conjGrad import conjGrad

A, b, x_sol = ae17.copy(), be17.copy(), xe17.copy()

def fAx(x):
    return A @ x

x_conjGrad, i_conjGrad = conjGrad(fAx, np.zeros_like(b), b)
```


```python
x_conjGrad, x_sol, i_conjGrad
```




    (array([-4.5, -4. , -3.5, -3. , -2.5, -2. , -1.5, -1. , -0.5,  0. ,  0.5,
             1. ,  1.5,  2. ,  2.5,  3. ,  3.5,  4. ,  4.5,  5. ]),
     array([-4.5, -4. , -3.5, -3. , -2.5, -2. , -1.5, -1. , -0.5,  0. ,  0.5,
             1. ,  1.5,  2. ,  2.5,  3. ,  3.5,  4. ,  4.5,  5. ]),
     9)




```python
x_conjGrad, i_conjGrad = conjGrad(fAx, np.ones_like(b), b, 1e-12)
x_conjGrad, x_sol, i_conjGrad
```




    (array([-4.50000000e+00, -4.00000000e+00, -3.50000000e+00, -3.00000000e+00,
            -2.50000000e+00, -2.00000000e+00, -1.50000000e+00, -1.00000000e+00,
            -5.00000000e-01, -2.22044605e-16,  5.00000000e-01,  1.00000000e+00,
             1.50000000e+00,  2.00000000e+00,  2.50000000e+00,  3.00000000e+00,
             3.50000000e+00,  4.00000000e+00,  4.50000000e+00,  5.00000000e+00]),
     array([-4.5, -4. , -3.5, -3. , -2.5, -2. , -1.5, -1. , -0.5,  0. ,  0.5,
             1. ,  1.5,  2. ,  2.5,  3. ,  3.5,  4. ,  4.5,  5. ]),
     9)




```python
import scipy.sparse.linalg

nb_iters = 0

def cg_callback(xk):
    global nb_iters
    nb_iters = nb_iters + 1
    print(f'xk: {xk}')

x_py, i_py = scipy.sparse.linalg.cg(A, b, np.zeros_like(b), callback=cg_callback)
x_py, i_py, nb_iters
```

    xk: [0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
     0.  0.5]
    xk: [-0.5  0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.
      0.   0.   0.   0.   0.5  1. ]
    xk: [-1.  -0.5  0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.
      0.   0.   0.   0.5  1.   1.5]
    xk: [-1.5 -1.  -0.5  0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.
      0.   0.   0.5  1.   1.5  2. ]
    xk: [-2.  -1.5 -1.  -0.5  0.   0.   0.   0.   0.   0.   0.   0.   0.   0.
      0.   0.5  1.   1.5  2.   2.5]
    xk: [-2.5 -2.  -1.5 -1.  -0.5  0.   0.   0.   0.   0.   0.   0.   0.   0.
      0.5  1.   1.5  2.   2.5  3. ]
    xk: [-3.  -2.5 -2.  -1.5 -1.  -0.5  0.   0.   0.   0.   0.   0.   0.   0.5
      1.   1.5  2.   2.5  3.   3.5]
    xk: [-3.5 -3.  -2.5 -2.  -1.5 -1.  -0.5  0.   0.   0.   0.   0.   0.5  1.
      1.5  2.   2.5  3.   3.5  4. ]
    xk: [-4.  -3.5 -3.  -2.5 -2.  -1.5 -1.  -0.5  0.   0.   0.   0.5  1.   1.5
      2.   2.5  3.   3.5  4.   4.5]
    xk: [-4.5 -4.  -3.5 -3.  -2.5 -2.  -1.5 -1.  -0.5  0.   0.5  1.   1.5  2.
      2.5  3.   3.5  4.   4.5  5. ]





    (array([-4.5, -4. , -3.5, -3. , -2.5, -2. , -1.5, -1. , -0.5,  0. ,  0.5,
             1. ,  1.5,  2. ,  2.5,  3. ,  3.5,  4. ,  4.5,  5. ]),
     0,
     10)




```python
%matplotlib widget
```


```python
import ipywidgets as widgets
```


```python
def three(x, y, z):
    return (x, y, z)
 
_ = widgets.interact(
    three, 
    x=(0, 10, 1), 
    y=True, 
    z=['a', 'b', 'c']
)
```


    interactive(children=(IntSlider(value=5, description='x', max=10), Checkbox(value=True, description='y'), Drop…



```python
# set up plot
fig, ax = plt.subplots(figsize=(6, 4))
ax.set_ylim([-4, 4])
ax.grid(True)
 
# generate x values
x = np.linspace(0, 2 * np.pi, 100)
 
 
def my_sine(x, w, amp, phi):
    """
    Return a sine for x with angular frequeny w and amplitude amp.
    """
    return amp*np.sin(w * (x-phi))
 
 
@widgets.interact(w=(0, 10, 1), amp=(0, 4, .1), phi=(0, 2*np.pi+0.01, 0.01))
def update(w = 1.0, amp=1, phi=0):
    """Remove old lines from plot and plot new one"""
    [l.remove() for l in ax.lines]
    ax.plot(x, my_sine(x, w, amp, phi), color='C0')
```


    Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …



    interactive(children=(IntSlider(value=1, description='w', max=10), FloatSlider(value=1.0, description='amp', m…



```python

```

# NUME

Instructor: Keny **ORDAZ**
    
# Final Exam

$\newcommand{\vb}[1]{\boldsymbol{#1}}\newcommand{\RR}[1]{\mathbb{R}^{#1}}
\newcommand{\vt}[1]{\vb{#1}^\top}$

$\newcommand{\x}{\vb{x}}
\newcommand{\xT}{\vt{x}}
\newcommand{\b}{\vb{b}}
\newcommand{\A}{\vb{A}}
\newcommand{\bT}{\vt{b}}
\newcommand{\AT}{\vt{A}}$ 

## Instructions

- You have exactly 24 hours to submit your examination
- The results and the solving process are to be evaluated
- Clarity of code and good programming practices must be observed
- It must be presented as a **self-contained Jupyter Python Notebook** by e-mail.
- Take all the appropriate precautions to deliver the exam on time. 
- Out-of-time deliveries will not be graded.
- Python libraries accepted: `numpy`, `scipy`, `sympy`, `matplotlib`.


```python

```

0ab4439bec804c0d20e5dd63f41986903be3b04c
ipynb
Jupyter Notebook
src/nume-s4-iterativelinalg.ipynb
krontzo/nume.py
9d1e576fb3474333a8e2cf4f26f4236ee4f9deea
MIT
src/nume-s4-iterativelinalg.ipynb
krontzo/nume.py
9d1e576fb3474333a8e2cf4f26f4236ee4f9deea
MIT
src/nume-s4-iterativelinalg.ipynb
krontzo/nume.py
9d1e576fb3474333a8e2cf4f26f4236ee4f9deea
MIT
Qwen/Qwen-72B
1. YES
2. YES