# NUME

Instructor: Keny **ORDAZ**
    
# Selected Topics

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


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

# Conjugate Gradient Visualization


```python
A = np.array([[4, -1], [-1, 3]])
b = np.array([9, -16])

np.linalg.solve(A, b)
```




    array([ 1., -5.])




```python


def conj_grad(Av, x, b, tol=1.0e-9, callback=None):
    ''' x, numIter = conj_grad(Av,x,b,tol=1.0e-9)
    Conjugate gradient method for solving [A]{x} = {b}.
    The matrix [A] should be sparse. User must supply
    the function Av(v) that returns the vector [A]{v}
    and the starting vector x.
''' 
    R = []
    n = len(b)
    r = b - Av(x)
    s = r.copy()
    R.append(x)
    for i in range(n):
        u = Av(s)
        alpha = np.dot(s,r)/np.dot(s,u)
        x = x + alpha*s
        r = b - Av(x)
        R.append(x)
        #print(x)
        if (np.linalg.norm(r)) < tol:
            break
        else:
            beta = -np.dot(r,u)/np.dot(s,u)
            s = r + beta*s
    return x, i, np.asarray(R)
```


```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, [0, 17], callback=cg_callback)
x_py, i_py, nb_iters
```

    xk: [ 6.83235818 -0.60646146]
    xk: [ 1. -5.]





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




```python
v_x, v_i, v_r = conj_grad(lambda x : A @ x, np.array([0., 17]), b)
```


```python
def visualization(A, b, v_r):
    N = 31
    x = np.linspace(-48, 48, N)
    y = np.linspace(-20, 20, N)

    X, Y = np.meshgrid(x, y)
    Z = np.zeros_like(Y)
    XY = np.array([X, Y])
    grad = np.zeros_like(XY)
    for i in range(N):
        for j in range(N):
            Z[i, j] = 0.5 * XY[:, i, j] @ A @ XY[:, i, j] - b @ XY[:, i, j]
            grad[:, i, j] = A @ XY[:, i, j] - b 
            
    U, V = grad
            
    fig, ax = plt.subplots()
    #cs = ax.contourf(X, Y, Z)
    CS = ax.contour(X, Y, Z)
    ax.quiver(X, Y, U, V)
    #P = np.array([[0, 17], [ 6.83235818, -0.60646146], [ 1., -5.]])
    
    #ax.plot(*P.T, 'r:o')
    ax.plot(*v_r.T, 'b:.')
    ax.set_aspect('equal')
    print(v_r.shape)

visualization(A, b, v_r)
```


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


    (3, 2)



```python
v_r[0] @ v_r[1], v_r[0] @ v_r[2]
v_r
```




    array([[ 2.60000000e+01, -6.70000000e+01],
           [-1.89358942e+01, -7.34825744e+00],
           [-8.88178420e-15,  7.10542736e-15]])




```python
P = np.array([[0, 17], [ 6.83235818, -0.60646146], [ 1., -5.]])
```


```python
x = P [-1]
x + v_r[1] + v_r[0]
```




    array([  8.06410583, -79.34825744])




```python

```


```python

```


```python

```

# Visualización


```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
# Repaso
```


```python
## LU
```


```python
A = np.array([[10, 4, 7.], [2, 5, 8], [3, 6, 10]])
A
```




    array([[10.,  4.,  7.],
           [ 2.,  5.,  8.],
           [ 3.,  6., 10.]])




```python
M = []
T = np.array([[1, 0, 0], [-1/5, 1, 0], [-3/10, 0, 1]])
M.append(T)
M[0] @ A
```




    array([[10. ,  4. ,  7. ],
           [ 0. ,  4.2,  6.6],
           [ 0. ,  4.8,  7.9]])




```python
T = np.array([[1, 0, 0], [0, 1, 0], [0, -8/7, 1]])
M.append(T)
M[1] @ M[0] @ A
```




    array([[10.        ,  4.        ,  7.        ],
           [ 0.        ,  4.2       ,  6.6       ],
           [ 0.        ,  0.        ,  0.35714286]])




```python
 M[1] M[0] A = (M[1] M[0] L) U = U
    => L = M[0]^{-1} M[1]^{-1}
```


```python
M
```




    [array([[ 1. ,  0. ,  0. ],
            [-0.2,  1. ,  0. ],
            [-0.3,  0. ,  1. ]]),
     array([[ 1.        ,  0.        ,  0.        ],
            [ 0.        ,  1.        ,  0.        ],
            [ 0.        , -1.14285714,  1.        ]])]




```python
[np.linalg.inv(M[i]) for i in range(2)]
```




    [array([[1. , 0. , 0. ],
            [0.2, 1. , 0. ],
            [0.3, 0. , 1. ]]),
     array([[ 1.        ,  0.        ,  0.        ],
            [-0.        ,  1.        , -0.        ],
            [ 0.        ,  1.14285714,  1.        ]])]




```python
np.any(np.isclose(, np.zeros(3)))
```




    False



1 0 1 == 0 0 0 ---> F V F ---> $!\exists A^{-1}$


```python
np.diag(M[1] @ M[0]).prod() == 0.0
```




    False




```python
np.isclose(np.diag(M[1] @ M[0]).prod(), 0.0)
```




    False




```python
a = 1e-300
f'{a:0.0f}', a == 0.0, (1. + a)/1.0 == 1.0
```




    ('0', False, True)




```python
from decimal import Decimal
Decimal.from_float(1e-17)
```




    Decimal('1.00000000000000007154242405462192450852805618492324772617063644020163337700068950653076171875E-17')




```python
Decimal.from_float((1. + a)/1.0), Decimal.from_float(a * b)
```




    (Decimal('1'), Decimal('5.46699999999999963762320476234890520572662353515625'))




```python
b = 5.467e300
a * b
```




    5.467




```python
0.1 + 0.1 + 0.1 == 0.3
```




    False




```python
A
```




    array([[10.,  4.,  7.],
           [ 2.,  5.,  8.],
           [ 3.,  6., 10.]])




```python
M
```




    [array([[ 1. ,  0. ,  0. ],
            [-0.2,  1. ,  0. ],
            [-0.3,  0. ,  1. ]]),
     array([[ 1.        ,  0.        ,  0.        ],
            [ 0.        ,  1.        ,  0.        ],
            [ 0.        , -1.14285714,  1.        ]])]




```python
Ac = A.copy()
Ac
```




    array([[10.,  4.,  7.],
           [ 2.,  5.,  8.],
           [ 3.,  6., 10.]])




```python
Ac[[1, 2]] = Ac[[2, 1]]
Ac
```




    array([[10.,  4.,  7.],
           [ 3.,  6., 10.],
           [ 2.,  5.,  8.]])




```python
P = np.eye(3)
P[[1, 2]] = P[[2, 1]]
P
```




    array([[1., 0., 0.],
           [0., 0., 1.],
           [0., 1., 0.]])




```python
P @ A
```




    array([[10.,  4.,  7.],
           [ 3.,  6., 10.],
           [ 2.,  5.,  8.]])




```python
Ac[1,:] += (-Ac[1,0]/Ac[0,0]) * Ac[0,:]
```


```python
Ac
```




    array([[10. ,  4. ,  7. ],
           [ 0. ,  4.8,  7.9],
           [ 2. ,  5. ,  8. ]])




```python
# P[1,:] += (-A[1,0]/A[0,0]) * P[0,:]
P
```




    array([[1., 0., 0.],
           [0., 0., 1.],
           [0., 1., 0.]])




```python
P @ P.T == np.eye(3)
# Ortho: P^T == P^{-1}
# Sim: A = A^T
#
```




    array([[ True,  True,  True],
           [ True,  True,  True],
           [ True,  True,  True]])




```python
  A = P @ L @ U
```


```python
x = np.array([5, 1, 0])
```


```python
P @ x
```




    array([5., 0., 1.])




```python
A
```




    array([[10.,  4.,  7.],
           [ 2.,  5.,  8.],
           [ 3.,  6., 10.]])




```python
Mi = [np.linalg.inv(M[i]) for i in range(2)]
```


```python
Minv = Mi[0] @ Mi[1]
Minv
```




    array([[1.        , 0.        , 0.        ],
           [0.2       , 1.        , 0.        ],
           [0.3       , 1.14285714, 1.        ]])




```python
D = np.diag([1,0,3])
np.linalg.det(Minv)
```




    1.0




```python
# A es triangular (diagonal)
#det(A) = ∏aii

det(A): O(n!)
A = LU => det(A) = det(LU) = det(L) det(U) = 1 * det(U) = ∏Uii, O(n)
LU: O(n^3)
    
    O(n^3) + O(n) = O(n^3) det por LU
```


```python
import math
[(i, i**3, math.factorial(i)) for i in range(15)]
```




    [(0, 0, 1),
     (1, 1, 1),
     (2, 8, 2),
     (3, 27, 6),
     (4, 64, 24),
     (5, 125, 120),
     (6, 216, 720),
     (7, 343, 5040),
     (8, 512, 40320),
     (9, 729, 362880),
     (10, 1000, 3628800),
     (11, 1331, 39916800),
     (12, 1728, 479001600),
     (13, 2197, 6227020800),
     (14, 2744, 87178291200)]




```python
Minv
```




    array([[1.        , 0.        , 0.        ],
           [0.2       , 1.        , 0.        ],
           [0.3       , 1.14285714, 1.        ]])




```python
M[1] @ M[0]
```




    array([[ 1.        ,  0.        ,  0.        ],
           [-0.2       ,  1.        ,  0.        ],
           [-0.07142857, -1.14285714,  1.        ]])




```python
Minv[[1,2],:][:,[1,2]]
```




    array([[1.        , 0.        ],
           [1.14285714, 1.        ]])




```python
1 () 
-0.2 () + 1()
0.3() - 1.14...() + 1()
```


```python
adjoint ---- adjugate
```


```python
Mi[0] @ Mi[1] = Ai @ Bi = (B A)^{-1} => M[1] @ M[0]
```

7c263bfa6e00c95eaefdf2b6e5a2d1c2131efa79
ipynb
Jupyter Notebook
src/nume-s5-selectedtopics.ipynb
krontzo/nume.py
9d1e576fb3474333a8e2cf4f26f4236ee4f9deea
MIT
src/nume-s5-selectedtopics.ipynb
krontzo/nume.py
9d1e576fb3474333a8e2cf4f26f4236ee4f9deea
MIT
src/nume-s5-selectedtopics.ipynb
krontzo/nume.py
9d1e576fb3474333a8e2cf4f26f4236ee4f9deea
MIT
Qwen/Qwen-72B
1. YES
2. YES