# Métodos de Newton-Raphson y Secante
## ¿Cuál es la ventaja del método de la Secante y Newton-Raphson en comparación a los vistos anteriormente?


```python
import numpy as np
import matplotlib.pyplot as plt
```


```python
def bisection(f, a, b, tol=1e-10):
    fa = f(a)
    fb = f(b)
    x = []
    
    # Just checking if the sign is not negative => not root  necessarily 
    if np.sign(f(a)*f(b)) >= 0:
        print('f(a)f(b)<0 not satisfied!')
        return None

    while (b-a) / 2 > tol:
        c = (a + b) / 2.
        x.append(c)
        fc = f(c)
        if fc == 0:
            break
        elif np.sign(fa*fc) < 0:
            b = c
            fb = fc
        else:
            a = c
            fa = fc
        
    return np.array(x)
```


```python
def fpi(g, x0, n, tol=1e-10):
    x = np.empty(n+1)
    x[0] = x0
    for i in range(n):
        x[i+1] = g(x[i])
        if np.abs(x[i+1] - x[i]) < tol:
            x = x[:i+2]
            break
        
    return x
```


```python
def newton_raphson(f, fp, x0, n, m=1, tol=1e-10):
    x = np.empty(n + 1)
    x[0] = x0
    for i in range(n):
        x[i+1] = x[i] - m * f(x[i]) / fp(x[i])
        if np.abs(x[i+1] - x[i]) < tol or f(x[i+1]) == 0:
            x = x[:i+2]
            break
    return x
```


```python
def secant(f, x0, x1, n, tol=1e-10):
    x = np.empty(n + 2)
    x[0] = x0
    x[1] = x1
    for i in range(1, n + 1):
        x[i+1] = x[i] - (f(x[i]) * (x[i] - x[i-1])) / (f(x[i]) - f(x[i-1]))
        if np.abs(x[i+1] - x[i]) < tol:
            x = x[:i+2]
            break
    return x
```

# Ejemplo

Encontrar una raíz de la ecuación:
\begin{equation}
    f(x) = x^2-x-1
\end{equation}

Sabemos que sus raíces son:
\begin{equation}
    x = \frac{1\pm\sqrt{5}}{2}
\end{equation}


```python
f = lambda x: x ** 2 - x - 1
```

Utilizamos los distintos algoritmos.

## Newton-Raphson


```python
x0 = 1
n = 15
fp = lambda x: 2 * x - 1
```


```python
x_n = newton_raphson(f, fp, x0, n)
```


```python
print(x_n.shape, x_n[-1], f(x_n[-1]))
```

    (7,) 1.618033988749895 0.0


## Secante


```python
x0 = 1
x1 = 1.5
n = 30
```


```python
x_s = secant(f, x0, x1, n)
```


```python
print(x_s.shape, x_s[-1], f(x_s[-1]))
```

    (8,) 1.6180339887498947 -2.220446049250313e-16


## Bisección


```python
a = 1
b = 2
```


```python
x_b = bisection(f, a, b, tol=1e-16)
```


```python
print(x_b.shape, x_b[-1], f(x_b[-1]))
```

    (49,) 1.618033988749895 0.0


## IPF


```python
# Funciones para x=g(x)
g1 = lambda x: 1 + 1 / x
g2 = lambda x: 1 / (x - 1)
g3 = lambda x: 1/2 *(x + 1/x + 1)
```


```python
x0 = 1
n = 50
```


```python
x_f = fpi(g1, x0, n, tol=1e-10)
```


```python
print(x_f.shape, x_f[-1], f(x_f[-1]))
```

    (27,) 1.618033988738303 -2.591993286671368e-11


# Análisis

Error conociendo la raíz:

\begin{equation}
    e = |x - r|
\end{equation}


```python
error = lambda x, r: np.abs(r - x)
```

Aproximación de la tasa de convergencia para el exponente $a$.

\begin{equation}
    tasa \approx \frac{|x_{i+1} - x_i|}{|x_i - x_{i-1}|^a}
\end{equation}


```python
# Utilizando el exponente correspondiente
conv_a = lambda x, a: error(x[1:-1], x[2:]) / error(x[:-2], x[1:-1]) ** a
```

Tasa conociendo $r$:

\begin{equation}
    tasa \approx \frac{|r - x_{i+1}|}{|r - x_i|^a}
\end{equation}


```python
conv_t = lambda x, r, e: error(x[1:], r) / error(x[:-1], r) ** e
```

## Convergencia

Aproximación del error 


```python
conv_b = conv_a(x_b, 0)
conv_f = conv_a(x_f, 0)
conv_n = conv_a(x_n, 0)
conv_s = conv_a(x_s, 0)
```


```python
plt.figure(figsize=(12, 6))
plt.plot(conv_b, 'ro', label="Bisección")
plt.plot(conv_f, 'yo', label="IPF")
plt.plot(conv_n, 'go', label="Newton-Raphson")
plt.plot(conv_s, 'bo', label="Secante")
plt.grid(True)
plt.legend()
plt.yscale('log')
plt.show()
```

## Aproximación tasa de convergencia

Calcular la tasa de convergencia según el método. Recordar:
- Lineal con tasa $S$ para Bisección e IPF
- Super-lineal para Secante.
- Cuadrática con tasa $M$ para Newton-Rapshon.


```python
conve_b = conv_a(x_b, 1) 
conve_f = conv_a(x_f, 1)
conve_n = conv_a(x_n, 2)
conve_s = conv_a(x_s, (1 + np.sqrt(5)) / 2)
```

Tabla para analizar los valores


```python
def table(approx, conv, conve, a):
    print("i \t x_{i} \t\t\t x_{i+1} \t\t |x_{i+1}-x_{i}| \t e_{i+1}/e_{i}^" + str(a))
    print("%d \t %.16f \t %.16f \t %.16f" % (0, approx[0], approx[1], error(approx[0], approx[1])))
    for i, (xip1, xi, c, ce) in enumerate(zip(approx[1:-1], approx[2:], conv, conve)):
        print("%d \t %.16f \t %.16f \t %.16f \t %.16f" % (i+1, xip1, xi, c, ce))
```

### Bisección


```python
table(x_b, conv_b, conve_b, a)
```

    i 	 x_{i} 			 x_{i+1} 		 |x_{i+1}-x_{i}| 	 e_{i+1}/e_{i}^1
    0 	 1.5000000000000000 	 1.7500000000000000 	 0.2500000000000000
    1 	 1.7500000000000000 	 1.6250000000000000 	 0.1250000000000000 	 0.5000000000000000
    2 	 1.6250000000000000 	 1.5625000000000000 	 0.0625000000000000 	 0.5000000000000000
    3 	 1.5625000000000000 	 1.5937500000000000 	 0.0312500000000000 	 0.5000000000000000
    4 	 1.5937500000000000 	 1.6093750000000000 	 0.0156250000000000 	 0.5000000000000000
    5 	 1.6093750000000000 	 1.6171875000000000 	 0.0078125000000000 	 0.5000000000000000
    6 	 1.6171875000000000 	 1.6210937500000000 	 0.0039062500000000 	 0.5000000000000000
    7 	 1.6210937500000000 	 1.6191406250000000 	 0.0019531250000000 	 0.5000000000000000
    8 	 1.6191406250000000 	 1.6181640625000000 	 0.0009765625000000 	 0.5000000000000000
    9 	 1.6181640625000000 	 1.6176757812500000 	 0.0004882812500000 	 0.5000000000000000
    10 	 1.6176757812500000 	 1.6179199218750000 	 0.0002441406250000 	 0.5000000000000000
    11 	 1.6179199218750000 	 1.6180419921875000 	 0.0001220703125000 	 0.5000000000000000
    12 	 1.6180419921875000 	 1.6179809570312500 	 0.0000610351562500 	 0.5000000000000000
    13 	 1.6179809570312500 	 1.6180114746093750 	 0.0000305175781250 	 0.5000000000000000
    14 	 1.6180114746093750 	 1.6180267333984375 	 0.0000152587890625 	 0.5000000000000000
    15 	 1.6180267333984375 	 1.6180343627929688 	 0.0000076293945312 	 0.5000000000000000
    16 	 1.6180343627929688 	 1.6180305480957031 	 0.0000038146972656 	 0.5000000000000000
    17 	 1.6180305480957031 	 1.6180324554443359 	 0.0000019073486328 	 0.5000000000000000
    18 	 1.6180324554443359 	 1.6180334091186523 	 0.0000009536743164 	 0.5000000000000000
    19 	 1.6180334091186523 	 1.6180338859558105 	 0.0000004768371582 	 0.5000000000000000
    20 	 1.6180338859558105 	 1.6180341243743896 	 0.0000002384185791 	 0.5000000000000000
    21 	 1.6180341243743896 	 1.6180340051651001 	 0.0000001192092896 	 0.5000000000000000
    22 	 1.6180340051651001 	 1.6180339455604553 	 0.0000000596046448 	 0.5000000000000000
    23 	 1.6180339455604553 	 1.6180339753627777 	 0.0000000298023224 	 0.5000000000000000
    24 	 1.6180339753627777 	 1.6180339902639389 	 0.0000000149011612 	 0.5000000000000000
    25 	 1.6180339902639389 	 1.6180339828133583 	 0.0000000074505806 	 0.5000000000000000
    26 	 1.6180339828133583 	 1.6180339865386486 	 0.0000000037252903 	 0.5000000000000000
    27 	 1.6180339865386486 	 1.6180339884012938 	 0.0000000018626451 	 0.5000000000000000
    28 	 1.6180339884012938 	 1.6180339893326163 	 0.0000000009313226 	 0.5000000000000000
    29 	 1.6180339893326163 	 1.6180339888669550 	 0.0000000004656613 	 0.5000000000000000
    30 	 1.6180339888669550 	 1.6180339886341244 	 0.0000000002328306 	 0.5000000000000000
    31 	 1.6180339886341244 	 1.6180339887505397 	 0.0000000001164153 	 0.5000000000000000
    32 	 1.6180339887505397 	 1.6180339886923321 	 0.0000000000582077 	 0.5000000000000000
    33 	 1.6180339886923321 	 1.6180339887214359 	 0.0000000000291038 	 0.5000000000000000
    34 	 1.6180339887214359 	 1.6180339887359878 	 0.0000000000145519 	 0.5000000000000000
    35 	 1.6180339887359878 	 1.6180339887432638 	 0.0000000000072760 	 0.5000000000000000
    36 	 1.6180339887432638 	 1.6180339887469017 	 0.0000000000036380 	 0.5000000000000000
    37 	 1.6180339887469017 	 1.6180339887487207 	 0.0000000000018190 	 0.5000000000000000
    38 	 1.6180339887487207 	 1.6180339887496302 	 0.0000000000009095 	 0.5000000000000000
    39 	 1.6180339887496302 	 1.6180339887500850 	 0.0000000000004547 	 0.5000000000000000
    40 	 1.6180339887500850 	 1.6180339887498576 	 0.0000000000002274 	 0.5000000000000000
    41 	 1.6180339887498576 	 1.6180339887499713 	 0.0000000000001137 	 0.5000000000000000
    42 	 1.6180339887499713 	 1.6180339887499144 	 0.0000000000000568 	 0.5000000000000000
    43 	 1.6180339887499144 	 1.6180339887498860 	 0.0000000000000284 	 0.5000000000000000
    44 	 1.6180339887498860 	 1.6180339887499002 	 0.0000000000000142 	 0.5000000000000000
    45 	 1.6180339887499002 	 1.6180339887498931 	 0.0000000000000071 	 0.5000000000000000
    46 	 1.6180339887498931 	 1.6180339887498967 	 0.0000000000000036 	 0.5000000000000000
    47 	 1.6180339887498967 	 1.6180339887498949 	 0.0000000000000018 	 0.5000000000000000


### IPF


```python
table(x_f, conv_f, conve_f, 1)
```

    i 	 x_{i} 			 x_{i+1} 		 |x_{i+1}-x_{i}| 	 e_{i+1}/e_{i}^1
    0 	 1.0000000000000000 	 2.0000000000000000 	 1.0000000000000000
    1 	 2.0000000000000000 	 1.5000000000000000 	 0.5000000000000000 	 0.5000000000000000
    2 	 1.5000000000000000 	 1.6666666666666665 	 0.1666666666666665 	 0.3333333333333330
    3 	 1.6666666666666665 	 1.6000000000000001 	 0.0666666666666664 	 0.3999999999999989
    4 	 1.6000000000000001 	 1.6250000000000000 	 0.0249999999999999 	 0.3750000000000000
    5 	 1.6250000000000000 	 1.6153846153846154 	 0.0096153846153846 	 0.3846153846153846
    6 	 1.6153846153846154 	 1.6190476190476191 	 0.0036630036630036 	 0.3809523809523809
    7 	 1.6190476190476191 	 1.6176470588235294 	 0.0014005602240896 	 0.3823529411764706
    8 	 1.6176470588235294 	 1.6181818181818182 	 0.0005347593582887 	 0.3818181818181530
    9 	 1.6181818181818182 	 1.6179775280898876 	 0.0002042900919306 	 0.3820224719101963
    10 	 1.6179775280898876 	 1.6180555555555556 	 0.0000780274656680 	 0.3819444444447312
    11 	 1.6180555555555556 	 1.6180257510729614 	 0.0000298044825942 	 0.3819742489263303
    12 	 1.6180257510729614 	 1.6180371352785146 	 0.0000113842055531 	 0.3819628647190943
    13 	 1.6180371352785146 	 1.6180327868852458 	 0.0000043483932688 	 0.3819672131268725
    14 	 1.6180327868852458 	 1.6180344478216819 	 0.0000016609364362 	 0.3819655522155682
    15 	 1.6180344478216819 	 1.6180338134001251 	 0.0000006344215568 	 0.3819661866826649
    16 	 1.6180338134001251 	 1.6180340557275543 	 0.0000002423274292 	 0.3819659445990088
    17 	 1.6180340557275543 	 1.6180339631667064 	 0.0000000925608479 	 0.3819660373066749
    18 	 1.6180339631667064 	 1.6180339985218035 	 0.0000000353550971 	 0.3819660027299530
    19 	 1.6180339985218035 	 1.6180339850173580 	 0.0000000135044456 	 0.3819660156562402
    20 	 1.6180339850173580 	 1.6180339901755971 	 0.0000000051582392 	 0.3819660085655353
    21 	 1.6180339901755971 	 1.6180339882053252 	 0.0000000019702719 	 0.3819659928497893
    22 	 1.6180339882053252 	 1.6180339889579018 	 0.0000000007525767 	 0.3819658851324578
    23 	 1.6180339889579018 	 1.6180339886704433 	 0.0000000002874585 	 0.3819657369182581
    24 	 1.6180339886704433 	 1.6180339887802426 	 0.0000000001097993 	 0.3819656758314164
    25 	 1.6180339887802426 	 1.6180339887383031 	 0.0000000000419396 	 0.3819657345316001


### Secante


```python
table(x_s, conv_s, conve_s, (1+np.sqrt(5))/2)
```

    i 	 x_{i} 			 x_{i+1} 		 |x_{i+1}-x_{i}| 	 e_{i+1}/e_{i}^1.618033988749895
    0 	 1.0000000000000000 	 1.5000000000000000 	 0.5000000000000000
    1 	 1.5000000000000000 	 1.6666666666666667 	 0.1666666666666667 	 0.5115940846088300
    2 	 1.6666666666666667 	 1.6153846153846154 	 0.0512820512820513 	 0.9311948386620529
    3 	 1.6153846153846154 	 1.6179775280898876 	 0.0025929127052722 	 0.3170353894921610
    4 	 1.6179775280898876 	 1.6180340557275543 	 0.0000565276376667 	 0.8646491584172478
    5 	 1.6180340557275543 	 1.6180339887482036 	 0.0000000669793507 	 0.4999508090467787
    6 	 1.6180339887482036 	 1.6180339887498947 	 0.0000000000016911 	 0.6855491046247754


### Newton-Raphson


```python
table(x_n, conv_n, conve_n, 2)
```

    i 	 x_{i} 			 x_{i+1} 		 |x_{i+1}-x_{i}| 	 e_{i+1}/e_{i}^2
    0 	 1.0000000000000000 	 2.0000000000000000 	 1.0000000000000000
    1 	 2.0000000000000000 	 1.6666666666666667 	 0.3333333333333333 	 0.3333333333333333
    2 	 1.6666666666666667 	 1.6190476190476191 	 0.0476190476190477 	 0.4285714285714292
    3 	 1.6190476190476191 	 1.6180344478216819 	 0.0010131712259371 	 0.4468085106382766
    4 	 1.6180344478216819 	 1.6180339887499890 	 0.0000004590716929 	 0.4472134119832668
    5 	 1.6180339887499890 	 1.6180339887498949 	 0.0000000000000941 	 0.4467299333062025



```python
plt.plot(conve_b, 'ro', label="Bisección")
plt.plot(conve_f, 'yo', label="IPF")
plt.plot(conve_n, 'go', label="Newton-Raphson")
plt.plot(conve_s, 'bo', label="Secante")
plt.grid(True)
plt.legend()
plt.yscale('log')
plt.show()
```

Estos valores aproximan $S$ y $M$ según corresponda.

## Convergencia "teórica"

Utilizando el valor real de la raíz.


```python
r = (1 + np.sqrt(5)) / 2
convt_b = conv_t(x_b, r, 1)
convt_f = conv_t(x_f, r, 1)
convt_n = conv_t(x_n, r, 2)
convt_s = conv_t(x_s, r, (1 + np.sqrt(5)) / 2)
```


```python
plt.plot(convt_b, 'ro', label="Bisección")
plt.plot(convt_f, 'yo', label="IPF")
plt.plot(convt_n, 'go', label="Newton-Raphson")
plt.plot(convt_s, 'bo', label="Secante")
plt.grid(True)
plt.legend
plt.yscale('log')
plt.show()
```

Acá se obtienen los $S$ y $M$ 'reales'.

## Error


```python
r = (1 + np.sqrt(5)) / 2
err_b = error(x_b, r)
err_f = error(x_f, r)
err_n = error(x_n, r)
err_s = error(x_s, r)
```


```python
plt.figure(figsize=(12, 6))
plt.plot(err_b, 'ro', label="Bisección")
plt.plot(err_f, 'yo', label="IPF")
plt.plot(err_n, 'go', label="Newton-Raphson")
plt.plot(err_s, 'bo', label="Secante")
plt.grid(True)
plt.legend()
plt.xlabel("# Iteraciones")
plt.ylabel("Error")
plt.yscale('log')
plt.show()
```

## Conclusión

El método de Newton-Raphson y el de la secante converge a la solución utilizando un menor número de iteraciones.

## 'Problema' de Newton-Raphson

Raíces con multiplicidad $m$.


```python
f2 = lambda x: x ** 2
```


```python
x0 = 1
n = 20
f2p = lambda x: 2 * x
```


```python
x_n2 = newton_raphson(f2, f2p, x0, n)
```


```python
print(x_n2[-1], f2(x_n2[-1]))
```

    9.5367431640625e-07 9.094947017729282e-13



```python
x0 = 1
x1 = 1.1
n = 20
```


```python
x_s2 = secant(f2, x0, x1, n)
```


```python
print(x_s2[-1], f2(x_s2[-1]))
```

    5.982324949014275e-05 3.5788211795598648e-09



```python
r = 0
err_n2 = error(x_n2, r)
err_s2 = error(x_s2, r)
```


```python
plt.figure(figsize=(12, 6))
plt.plot(err_n2, 'go', label="Newton-Raphson")
plt.plot(err_s2, 'bo', label="Secante")
plt.grid(True)
plt.legend()
plt.xlabel("# Iteraciones")
plt.ylabel("Error")
plt.yscale('log')
plt.show()
```

El problema es que para este tipo de problemas se obtiene convergencia lineal. Se puede corregir con *Newton-Raphson* modificado.

## Utilizando $m$ para Newton-Raphson


```python
n = 5
x0 = 1.1
x_n2_m = newton_raphson(f2, f2p, x0, n, 2)
```


```python
x_n2_m
```




    array([1.1, 0. ])




```python
r = 0
err_n2_m = error(x_n2_m, r)
```


```python
plt.figure(figsize=(12, 6))
plt.plot(err_n2_m, 'go', label="Newton-Raphson Mejorado")
plt.plot(err_n2, 'gx', label="Newton-Raphson")
plt.plot(err_s2, 'bo', label="Secante")
plt.grid(True)
plt.legend()
plt.xlabel("# Iteraciones")
plt.ylabel("Error")
plt.show()
```

Se alcanzó la solución en una interación!

## Otro ejemplo

\begin{equation}
    f(x) = 0.5 x\cos(x)
\end{equation}


```python
f3 = lambda x: .5 * x * np.cos(x) 
```


```python
x = np.linspace(-10, 10, 100)
```


```python
plt.plot(x, f3(x))
plt.grid(True)
plt.show()
```


```python
x0 = -5
x1 = -4
n = 20
xx_s = secant(f3, x0, x1, n)
```


```python
print(xx_s.shape, xx_s[-1], f3(xx_s[-1]))
```

    (8,) -4.71238898038469 4.3282590608740267e-16



```python
x0 = 1
n = 20
f3p = lambda x: .5 * np.cos(x) - .5 * x * np.sin(x)
xx_n = newton_raphson(f3, f3p, x0, n)
```


```python
print(xx_n.shape, xx_n[-1], f3(xx_n[-1]))
```

    (8,) 1.5707963267948966 4.8091767343044744e-17


# Referencias

* Sauer, T. (2006). Numerical Analysis Pearson Addison Wesley.
* https://github.com/tclaudioe/Scientific-Computing/tree/master/SC1

2a0fe272bb9942ffdc2fb5a9654eb53ce60ce6c9
ipynb
Jupyter Notebook
material/03_raices_1D/newton_secant.ipynb
etra0/INF-285
189f8d66cf6997fc87b545378b2c5f1c7b908dbc
BSD-3-Clause
2020-04-24T01:25:24.000Z
2021-05-10T01:08:37.000Z
material/03_raices_1D/newton_secant.ipynb
etra0/INF-285
189f8d66cf6997fc87b545378b2c5f1c7b908dbc
BSD-3-Clause
2020-04-22T00:57:29.000Z
2020-04-23T23:59:28.000Z
material/03_raices_1D/newton_secant.ipynb
etra0/INF-285
189f8d66cf6997fc87b545378b2c5f1c7b908dbc
BSD-3-Clause
2020-04-20T01:15:09.000Z
2020-07-27T22:53:59.000Z
Qwen/Qwen-72B
1. YES
2. YES