2.4.1 Método de Newton-Raphson - Ejemplo con Python

[ Newton-Raphson ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [gráf] [función]
..


1. Ejercicio

ReferenciaBurden 2.1 ejemplo 1 p38

La ecuación mostrada tiene una raíz en [1,2], ya que f(1)=-5 y f(2)=14.
Muestre los resultados parciales del algoritmo de Newton-Raphson con una tolerancia de 0.0001

f(x) = x^3 + 4x^2 -10 =0

Newton Raphson GIF animado
[ Newton-Raphson ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [gráf] [función]
..


2. Desarrollo Analítico

El método requiere  obtener la derivada f'(x) de la ecuación para el factor del denominador.

f(x) = x^3 + 4x^2 -10 f'(x) = 3x^2 + 8x x_{i+1} = x_i -\frac{f(x_i)}{f'(x_i)}

Para el desarrollo se inicia la búsqueda desde un punto en el intervalo [1,2], por ejemplo el extremo derecho, x1=2.

iteración 1

f(2) = (2)^3 + 4(2)^2 -10 = 14 f'(2) = 3(2)^2 + 8(2) = 28 x_{2} = 2 -\frac{14}{28} = 1.5 tramo = |2 -1.5| = 0.5

iteración 2

f(1.5) = (1.5)^3 + 4(1.5)^2 -10 = 2.375 f'(1.5) = 3(1.5)^2 + 8(1.5) = 18.75 x_{3} = 1.5 -\frac{2.375}{18.75} = 1.3733 tramo = |1.5 -1.3733| = 0.1267

iteración 3

f(1.3733) = (1.3733)^3 + 4(1.3733)^2 -10 = 0.1337 f'(1.3733) = 3(1.3733)^2 + 8(1.3733) = 16.6442 x_{4} = 1.3733 -\frac{0.1337}{16.6442} =1.3652 tramo = |1.3733 -1.3652| = 0.0081

La tabla resume los valores de las iteraciones

Método de Newton-Raphson
iteración xi xnuevo tramo
1 2 1.5 0.5
2 1.5 1.3733 0.1267
3 1.3733 1.3653 0.0081
4 ...

Observe que el error representado por el tramo se va reduciendo entre cada iteración. Se debe repetir las iteraciones hasta que el error sea menor al valor tolerado.

Las demás iteraciones se dejan como tarea

[ Newton-Raphson ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [gráf] [función]
..


3. Algoritmo con Python

El método de Newton-Raphson se implementa como algoritmo básico en Python

https://youtu.be/GjgzoafJKpU

Se muestra el resultado del algoritmo luego de que el tramo alcance un valor menor que tolera.

raiz en:  1.3652300139161466
con error de:  3.200095847999407e-05

A lo expuesto en el video, se añade el control de iteraciones "iteramax" por si se da el caso que el algoritmo no es convergente.

# Método de Newton-Raphson
import numpy as np

# INGRESO
fx  = lambda x: x**3 + 4*(x**2) - 10
dfx = lambda x: 3*(x**2) + 8*x

x0 = 2
tolera = 0.001
iteramax = 100

# PROCEDIMIENTO
itera = 0
tramo = abs(2*tolera)
xi = x0
while (tramo>=tolera and itera<iteramax):
    fi = fx(xi)
    dfi = dfx(xi)
    xnuevo = xi - fi/dfi
    tramo = abs(xnuevo-xi)
    xi = xnuevo
    itera = itera + 1

if itera>=iteramax:
    xi = np.nan

# SALIDA
print('raiz en: ', xi)
print('con error de: ',tramo)

[ Newton-Raphson ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [gráf] [función]
..


4. Gráfica en Python

La gráfica presentada para revisar f(x) se realiza con las instrucciones:

# GRAFICA
import matplotlib.pyplot as plt
a = 1
b = 2
muestras = 21

xj = np.linspace(a,b,muestras)
fj = fx(xj)
plt.plot(xj,fj, label='f(x)')
plt.plot(xi,0, 'o')
plt.axhline(0)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.grid()
plt.legend()
plt.show()

[ Newton-Raphson ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [gráf] [función]
..


5. Algoritmo en Python como Función

Se convierte el algoritmo a una función, con partes para ver la tabla, y se obtienen los siguientes resultados:

i ['xi', 'fi', 'dfi', 'xnuevo', 'tramo']
0 [ 2.  14.  28.   1.5  0.5]
1 [ 1.5     2.375  18.75    1.3733  0.1267]
2 [1.3733e+00 1.3435e-01 1.6645e+01 1.3653e+00 8.0713e-03]
3 [1.3653e+00 5.2846e-04 1.6514e+01 1.3652e+00 3.2001e-05]
raíz en:  1.3652300139161466

Instrucciones en Python

# Ejemplo 1 (Burden ejemplo 1 p.51/pdf.61)

import numpy as np
import matplotlib.pyplot as plt

def newton_raphson(fx,dfx,xi, tolera, iteramax=100,
                   vertabla=False, precision=4):
    '''fx y dfx en forma numérica lambda
    xi es el punto inicial de búsqueda
    '''
    itera=0
    tramo = abs(2*tolera)
    if vertabla==True:
        print('método de Newton-Raphson')
        print('i', ['xi','fi','dfi', 'xnuevo', 'tramo'])
        np.set_printoptions(precision)
    while (tramo>=tolera):
        fi = fx(xi)
        dfi = dfx(xi)
        xnuevo = xi - fi/dfi
        tramo = abs(xnuevo-xi)
        if vertabla==True:
            print(itera,np.array([xi,fi,dfi,xnuevo,tramo]))
        xi = xnuevo
        itera = itera + 1

    if itera>=iteramax:
        xi = np.nan
        print('itera: ',itera,
              'No converge,se alcanzó el máximo de iteraciones')

    return(xi)

# INGRESO
fx  = lambda x: x**3 + 4*(x**2) - 10
dfx = lambda x: 3*(x**2) + 8*x

x0 = 2
tolera = 0.001

# PROCEDIMIENTO
respuesta = newton_raphson(fx,dfx,x0, tolera, vertabla=True)
# SALIDA
print('raíz en: ', respuesta)

La gráfica se la puede añadir usando las instrucciones dadas en el algoritmo básico realizado al inicio par ala comprensión del método.

[ Newton-Raphson ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [gráf] [función]


6. Función en librería scipy.optimize.newton

El método de Newton-Raphson se encuentra implementado en la librería  Scipy, que también puede ser usado de la forma:

>>> import scipy as sp
>>> sp.optimize.newton(fx,x0, fprime=dfx, tol = tolera)
1.3652300139161466
>>> 

https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.newton.html


Tarea

Calcule la raíz de f(x) = e-x - x, empleando como valor inicial x0 = 0

  • Revisar las modificaciones si se quiere usar la forma simbólica de la función y obtener la derivada con Sympy.

[ Newton-Raphson ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [gráf] [función]

2.4 Método de Newton-Raphson - Concepto

[ Newton-Raphson ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ]
..


Método de Newton-Raphson

Referencia: Burden 2.3 p49, Chapra 6.2 p148, Rodríguez 3.3 p52

Se deduce a partir de la interpretación gráfica o por medio del uso de la serie de Taylor de dos términos.

Newton Raphson GIF animado

De la gráfica, se usa el triángulo formado por la recta tangente que pasa por f(xi), con pendiente f'(xi)  y el eje x.

Newton Raphson GIF_animado

f'(x_i) = \frac{f(x_i) - 0}{x_i - x_{i+1}}

El punto xi+1 es la intersección de la recta tangente con el eje x, que es más cercano a la raíz de f(x), valor que es usado para la próxima iteración.

Reordenando la ecuación de determina la fórmula para el siguiente punto:

x_{i+1} = x_i -\frac{f(x_i)}{f'(x_i)}

El error se determina como la diferencia entre los valores sucesivos encontrados |xi+1 - xi|

La gráfica animada muestra el proceso aplicado varias veces sobre f(x) para encontrar la raíz.


Recta Tangente para la gráfica

Si se quiere dibujar la recta tangente en el punto inicial, la ecuación de la recta se obtiene usando: el valor de la pendiente con la derivada f'(x) y cambiando la constante b por b0 para no confundirla con la b del intervalo [a,b]

y = mx + b y = f'(x) x + b_0

Es necesario disponer de un punto conocido de la recta (x0,f(x0)) y la pendiente en ese punto f'(x0), quedando solo el término de la constante como incógnita .

f(x_0) = f'(x_0) x_0 + b_0 f(x_0) - f'(x_0) x_0 = b_0

La función recta tangente se expresa para la gráfica como:

b_0 = f(x_0) - f'(x_0) x_0 rtg(x)= f'(x_0) x + b_0

Tarea

Use la serie de Taylor hasta la primera derivada para encontrar el siguiente punto de aproximación xi+1

[ Newton-Raphson ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ]

2.3.1 Método del Punto fijo - Ejemplo con Python

[ Punto fijo ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ gráfica ] [ función ]
..


1. Ejercicio

Referencia: Burden 2.2 p41, Chapra 6.1 p143, Rodríguez 3.2 p44

Encontrar la solución a la ecuación, usando el método del punto fijo, considerando tolerancia de 0.001

f(x): e^{-x} - x = 0

Punto Fijo 01a_GIF

[ Punto fijo ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ gráfica ] [ función ]

..


2. Desarrollo Analítico

Al igual que los métodos anteriores, es conveniente determinar el intervalo [a,b] donde es posible evaluar f(x). Se revisa si hay cambio de signo en el intervalo para buscar una raíz. Para el ejemplo, el rango de observación será [0,1], pues f(0)=1 es positivo y f(1)=-0.63 es negativo.Punto Fijo Balanza animado

Para el punto fijo, se reordena la ecuación para para tener una ecuación con la variable independiente separada. Se obtiene por un lado la recta identidad y=x, por otro se tiene la función g(x).

f(x): e^{-x} - x = 0 x = e^{-x}

g(x) = e^{-x}punto fijo 01b

Se buscará la intersección entre las dos expresiones .

Se puede iniciar la búsqueda por uno de los extremos del rango [a,b].

iteración 1

- iniciando desde el extremo izquierdo c=a,

c = 0

- se determina el valor nuevo de gc=g(c) ,

gc = g(0) = e^{-0} = 1

- se determina la diferencia de la aproximación o error = |gc-c|

tramo = |1-0| = 1

- se proyecta en la recta identidad, como el nuevo punto de evaluación
c = gc.

- se repite el proceso para el nuevo valor de c, hasta que el error sea menor al tolerado. En caso que el proceso no converge, se utiliza un contador de iteraciones máximo, para evitar tener un lazo infinito.

iteración 2

c = 1 gc = g(1) = e^{-1} = 0.3678 tramo =|0.3678 - 1|= 0.6322

iteración 3

c = 0.3678 gc = e^{-0.3678} = 0.6922 tramo =|0.6922-0.3678|= 0.3244

La tabla resume los valores de las iteraciones

Método del punto fijo
iteración c gc= g(c) |tramo|
1 0 1.0 1
2 1.0 0.3678 0.6322
3 0.3678 0.6922 0.3244
4 0.6922 ...
5

El proceso realizado en la tabla se muestra en la gráfica, para una función que converge.

Cuando el método no converge

Cuando la 'x' a separar en la expresión f(x) es otra, el resultado no siempre será convergente. Si para el ejercicio anterior se hubiese despejado la 'x' en el exponencial, la expresión g(x) se obtiene como:

f(x): e^{-x} - x = 0 e^{-x} = x \ln\Big(e^{-x}\Big) = \ln(x) -x= \ln(x) x= -\ln(x)

Al revisar la gráfica del algoritmo se observa que también se cumple que la raíz de f(x) coincide en el valor x con la intersección entre la identidad y g(x).

punto fijo No Converge

Sin embargo al realizar las iteraciones se tiene que:

iteración 1
c = 0.5
gc = -ln(0.5) = 0.6931
tramo = | 0.6931-0.5| =0.1931

iteración 2
c = 0.6931
gc = -ln(0.6931) = 0.3665
tramo = | 0.3665-0.6931| = 0.3266

iteración 3
 c = 0.3665
gc = -ln(0.3665) = 1.0037
tramo = | 1.0037-0.3665| = 0.6372

Se observa que el tramo o error aumenta en cada iteración, con lo que se concluye que de esta forma, la aplicación del método no es convergente.

[ Punto fijo ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ gráfica ] [ función ]
..


3. Algoritmo en Python

https://youtu.be/ID_5uXfOal8

El algoritmo en Python del video se adjunta para seguimiento. Se requiere la función gx, el punto inicial y la tolerancia, la variable de número de iteraciones máxima, iteramax, permite controlar la convergencia de la función con el método.

El resultado del algoritmo se presenta como:

raíz en:  0.5669089119214953
errado :  0.0006477254067881466
itera  :  14

Instrucciones en Python

# Algoritmo de punto fijo
# error = tolera
import numpy as np

# INGRESO
fx = lambda x: np.exp(-x) - x
gx = lambda x: np.exp(-x)

c = 0  # valor inicial
tolera = 0.001
iteramax = 40

# PROCEDIMIENTO
tramo = 2*tolera # al menos una iteracion
i = 0  # iteración inicial
while tramo>=tolera and i<iteramax:
    gc = gx(c)
    tramo = abs(gc-c)
    c = gc
    i = i + 1

if (i>=iteramax): # Valida convergencia
    c = np.nan

# SALIDA
print('raíz en: ', c)
print('errado : ', tramo)
print('itera  : ', i)

Para obtener la gráfica básica se determinan los puntos para cada función fx y gx. Como los valores de c y gc cambiaron durante el algoritmo, para la siguiente sección es necesario volver a escribirlos u obtener una copia de los valores en otra variable para reutilizarlos en ésta sección.

[ Punto fijo ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ gráfica ] [ función ]
..


4. Gráfica en Python

Se añaden las siguiente instrucciones al algoritmo anterior:

# GRAFICA
import matplotlib.pyplot as plt

a = 0  # intervalo de gráfica en [a,b]
b = 1
muestras = 21

# calcula los puntos para fx y gx
xi = np.linspace(a,b,muestras)
fi = fx(xi)
gi = gx(xi)
yi = xi

plt.plot(xi,fi, label='f(x)',
         linestyle='dashed')
plt.plot(xi,gi, label='g(x)')
plt.plot(xi,yi, label='y=x')

plt.axvline(c, color='magenta',
            linestyle='dotted')
plt.axhline(0)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('Punto Fijo')
plt.grid()
plt.legend()
plt.show()

[ Punto fijo ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ gráfica ] [ función ]
..


5. Algoritmo en Python como Función

El resultado del algoritmo como una función y mostrando la tabla de resultados por iteraciones es:

método del Punto Fijo
i ['xi', 'gi', 'tramo']
0 [0. 1. 1.]
1 [1.       0.367879 0.632121]
2 [0.367879 0.692201 0.324321]
3 [0.692201 0.500474 0.191727]
4 [0.500474 0.606244 0.10577 ]
5 [0.606244 0.545396 0.060848]
6 [0.545396 0.579612 0.034217]
7 [0.579612 0.560115 0.019497]
8 [0.560115 0.571143 0.011028]
9 [0.571143 0.564879 0.006264]
10 [0.564879 0.568429 0.003549]
11 [0.568429 0.566415 0.002014]
12 [0.566415 0.567557 0.001142]
13 [0.567557 0.566909 0.000648]
raíz en:  0.5669089119214953

Instrucciones en Python

# Algoritmo de punto fijo
# error = tolera
import numpy as np

def puntofijo(gx,c,tolera,iteramax=40,vertabla=True, precision=6):
    """
    g(x) se obtiene al despejar una x de f(x)
    máximo de iteraciones predeterminado: iteramax
    si no converge hasta iteramax iteraciones
    la respuesta es NaN (Not a Number)
    """
    itera = 0 # iteración inicial
    tramo = 2*tolera # al menos una iteracion
    if vertabla==True:
        print('método del Punto Fijo')
        print('i', ['xi','gi','tramo'])
        np.set_printoptions(precision)
    while (tramo>=tolera and itera<=iteramax):
        gc = gx(c)
        tramo = abs(gc-c)
        if vertabla==True:
            print(itera,np.array([c,gc,tramo]))
        c = gc
        itera = itera + 1
    respuesta = c
    # Valida respuesta
    if itera>=iteramax:
        respuesta = np.nan
        print('itera: ',itera,
              'No converge,se alcanzó el máximo de iteraciones')
    return(respuesta)

# PROGRAMA ----------------------
# INGRESO
fx = lambda x: np.exp(-x) - x
gx = lambda x: np.exp(-x)

c = 0  # valor inicial
tolera = 0.001
iteramax = 15

# PROCEDIMIENTO
respuesta = puntofijo(gx,c,tolera,iteramax,vertabla=True)

# SALIDA
print('raíz en: ', respuesta)

De añadirse la parte gráfica, el bloque requiere el intervalo de observación [a,b].

[ Punto fijo ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ gráfica ] [ función ]


6. Función en librería Scipy.optimize.fixed_point

El método del punto fijo se encuentra implementado en Scipy, que también puede ser usado de la forma:

>>> import numpy as np
>>> import scipy as sp
>>> gx = lambda x: np.exp(-x)
>>> c = 0
>>> sp.optimize.fixed_point(gx,c,xtol=0.001,maxiter=15)
array(0.5671432948307147)

el valor predeterminado de iteraciones si no se escribe es 500 y la tolerancia predeterminada es xtol=1e-08
https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fixed_point.html


Tarea

  • Revisar lo que sucede cuando el valor inicial a esta a la derecha de la raíz.
  • Validar que la función converge, revisando que |g'(x)|<1, o que la función g(x) tenga pendiente menor a pendiente de la recta identidad.

[ Punto fijo ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ gráfica ] [ función ]

2.3 Método del Punto fijo - Concepto

[ Punto fijo ][ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ]
..


1. Método del Punto fijo

Referencia: Burden 2.2 p41, Chapra 6.1 p143, Rodríguez 3.2 p44

El método del punto fijo es un método abierto, también llamado de iteración de un punto o sustitución sucesiva, que reordena la ecuación planteada para raíces de ecuacionesPunto Fijo Balanza animado

 f(x)= 0

separando una x al lado izquierdo de la igualdad.

x = g(x)

Esto permite buscar la intersección entre la recta Identidad y=x y la curva g(x), que es lo que quedó del lado derecho de la igualdad.
Una representación gráfica del proceso de muestra en la gráfica.

Punto Fijo 01a_GIF

Observe que la raíz de f(x) se encuentra en el mismo valor de x donde ocurre la intersección entre la recta identidad y=x , en color verde, y la función g(x), en color naranja. Como referencia, se usa la linea vertical en x=raíz en color morado.

El método consiste en establecer un punto inicial x0 para la búsqueda, a partir del cual se calcula el valor g(x0). En la siguiente iteración el nuevo valor para x es g(x0), que se refleja en la recta identidad y nuevamente se usa para calcular g(x).

El resultado iterativo se muestra en la figura animada, donde se observa que el resultado es convergente.

[ Punto fijo ][ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ]


2. Ejemplos

Se muestran algunos ejemplos para destacar lo indicado de forma gráfica.

Ejemplo 1

f(x):e^{-x} - x = 0

se reordena para tener:

x = e^{-x} g(x) = e^{-x}

Punto Fijo


Ejemplo 2

f(x): x^2 - 2x -3 = 0

se reordena para tener:

x = \frac {x^2 - 3}{2} g(x) = \frac {x^2 - 3}{2}

Punto Fijo


Ejemplo 3

f(x): \sin (x) = 0

puede ser complicado despejar x, por lo que se simplifica el proceso sumando x en ambos lados.

x = \sin (x) + x g(x) = \sin (x) + x

punto fijo 03

El método proporciona una fórmula para predecir un valor nuevo de x en función del valor anterior:

x_{i+1} = g(x_i)

con error aproximado calculado como:

\epsilon_a = \left| \frac{x_{i+1} - x_i}{x_{i+1}} \right| 100\%

[ Punto fijo ][ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ]


4. Tarea

Plantee como usar los siguientes conceptos:

  • ¿cuál sería el valor de tolerancia?
  • ¿parámetros de inicio?
  • compare con con otro método conocido
  • Revisar el resultado cuando no se cumple que |g'(x)|<1

[ Punto fijo ][ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ]

2.2.1 Método de la Falsa Posición - Ejemplo con Python

[ Falsa Posición ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ][gráfica]
..


1. Ejercicio

Referencia: Burden 2.1 ejemplo 1 p38

La ecuación mostrada tiene una raíz en [1,2], ya que f(1)=-5 y f(2)=14. Muestre los resultados parciales del algoritmo de la falsa posición con una tolerancia de 0.0001

f(x) = x^3 + 4x^2 -10 =0

Posicion Falsa animado GIF

[ Falsa Posición ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ][gráfica]
..


2. Desarrollo Analítico

Semejante a los métodos anteriores, el método posición falsa, falsa posición, regla falsa o regula falsi, usa un intervalo [a,b] para buscar la raíz.

Se divide el intervalo en dos partes al calcular el punto c que divide al intervalo siguiendo la ecuación:

c = b - f(b) \frac{a-b}{f(a)-f(b)}

iteración 1

a = 1 , b = 2 f(1) = (1)^3 + 4(1)^2 -10 = -5 f(2) = (2)^3 + 4(2)^2 -10 = 14 c = 2 - 14 \frac{1-2}{-5-14} = 1.2631 f(1.2631) = (1.2631)^3 + 4(1.2631)^2 -10 = -1.6031

el signo de f(c) es el mismo que f(a), se ajusta el lado izquierdo

tramo = |c-a| = | 1.2631 - 1| = 0.2631 a = c = 1.2631

iteración 2

a = 1.2631 , b = 2 f(1.2631) = -1.6031 f(2) = 14 c = 2 - 14 \frac{1.2631-2}{-1.6031-14} = 1.3388 f(1.3388) = (1.3388)^3 + 4(1.3388)^2 -10 = -0.4308

el signo de f(c) es el mismo que f(a), se ajusta el lado izquierdo

tramo = |c-a| = |1.3388 - 1.2631| = 0.0757 a = c = 1.3388

iteración 3

a = 1.3388 , b = 2 f(1.3388) = -0.4308 f(2) = 14 c = 2 - 14 \frac{1.3388-2}{-0.4308-14} = 1.3585 f(1.3585) = (1.3585)^3 + 4(1.3585)^2 -10 = -0.1107

el signo de f(c) es el mismo que f(a), se ajusta el lado izquierdo

tramo = |c-a| = |1.3585 - 1.3388| = 0.0197 a = c = 1.3585

valores que se resumen en la tabla

Método de posición falsa
a c b f(a) f(c) f(b) tramo
1 1.2631 2 -5 -1.6031 14 0.2631
1.2631 1.3388 2 -1.6031 -0.4308 14 0.0757
1.3388 1.3585 2 -0.4308 -0.1107 14 0.0197
1.3585 ... 2

se puede continuar con las iteraciones como tarea

[ Falsa Posición ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ][gráfica]
..


3. Algoritmo en Python

https://youtu.be/jYXmFepPtSk

Algoritmo básico del video:

# Algoritmo Posicion Falsa para raices
# busca en intervalo [a,b]
import numpy as np

# INGRESO
fx = lambda x: x**3 + 4*x**2 - 10

a = 1
b = 2
tolera = 0.001

# PROCEDIMIENTO
tramo = abs(b-a)
while not(tramo<=tolera):
    fa = fx(a)
    fb = fx(b)
    c = b - fb*(a-b)/(fa-fb)
    fc = fx(c)
    cambia = np.sign(fa)*np.sign(fc)
    if (cambia > 0):
        tramo = abs(c-a)
        a = c
    else:
        tramo = abs(b-c)
        b = c
raiz = c

# SALIDA
print(raiz)

Algoritmo aumentado para mostrar la tabla de cálculos

# Algoritmo Posicion Falsa para raices
# busca en intervalo [a,b]
# tolera = error

import numpy as np

# INGRESO
fx = lambda x: x**3 + 4*(x**2) -10

a = 1
b = 2
tolera = 0.0001

# PROCEDIMIENTO
tabla = []
tramo = abs(b-a)
fa = fx(a)
fb = fx(b)
while not(tramo<=tolera):
    c = b - fb*(a-b)/(fa-fb)
    fc = fx(c)
    unafila = [a,c,b,fa,fc,fb,tramo]
    cambio = np.sign(fa)*np.sign(fc)
    if cambio>0:
        tramo = abs(c-a)
        a = c
        fa = fc
    else:
        tramo = abs(b-c)
        b = c
        fb = fc
    unafila[6]=tramo
    tabla.append(unafila)
tabla = np.array(tabla)
ntabla = len(tabla)

# SALIDA
np.set_printoptions(precision=6)
print('método de la Posición Falsa ')
print('i', ['a','c','b'],[ 'f(a)', 'f(c)','f(b)'])
print('  ','tramo')
for i in range(0,ntabla,1):
    print(i, tabla[i,0:3],tabla[i,3:6])
    print('   ', tabla[i,6])

print('raiz:  ',c)

Observe el número de iteraciones realizadas, hasta presentar el valor de la raíz en 1.3652 con un error de 0.000079585 en la última fila de la tabla. Sin embargo, observe que la tabla solo muestra cálculos de filas completas, el último valor de c y error no se ingresó a la tabla, que se muestra como c y tramo, y es el más actualizado en los cálculos.

método de la Posición Falsa 
i ['a', 'c', 'b'] ['f(a)', 'f(c)', 'f(b)']
   tramo
0 [1.       1.263158 2. ] [-5.       -1.602274 14. ]
   0.26315789473684204
1 [1.263158 1.338828 2. ] [-1.602274 -0.430365 14. ]
   0.0756699440909967
2 [1.338828 1.358546 2. ] [-0.430365 -0.110009 14. ]
   0.019718502996940224
3 [1.358546 1.363547 2. ] [-0.110009 -0.027762 14. ]
   0.005001098217311428
4 [1.363547 1.364807 2. ] [-2.776209e-02 -6.983415e-03  1.400000e+01]
   0.0012595917846898175
5 [1.364807 1.365124 2. ] [-6.983415e-03 -1.755209e-03  1.400000e+01]
   0.0003166860575976038
6 [1.365124 1.365203 2. ] [-1.755209e-03 -4.410630e-04  1.400000e+01]
    7.958577822231305e-05
raíz en:  1.3652033036626001

[ Falsa Posición ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ][gráfica]
..


4. Función en Python

# Algoritmo de falsa posicion para raices
# Los valores de [a,b] son seleccionados
# desde la gráfica de la función
# error = tolera

import numpy as np

def posicionfalsa(fx,a,b,tolera,iteramax = 20,
                        vertabla=False, precision=6):
    '''fx en forma numérica lambda
    Los valores de [a,b] son seleccionados
    desde la gráfica de la función
    error = tolera
    '''
    fa = fx(a)
    fb = fx(b)
    tramo = np.abs(b-a)
    itera = 0
    cambia = np.sign(fa)*np.sign(fb)

    if cambia<0: # existe cambio de signo f(a) vs f(b)
        if vertabla==True:
            print('método de la Posición Falsa ')
            print('i', ['a','c','b'],[ 'f(a)', 'f(c)','f(b)'])
            print('  ','tramo')
            np.set_printoptions(precision)

        while (tramo >= tolera and itera<=iteramax):
            c = b - fb*(a-b)/(fa-fb)
            fc = fx(c)
            cambia = np.sign(fa)*np.sign(fc)
            if vertabla==True:
                print(itera,np.array([a,c,b]),
                      np.array([fa,fc,fb]))
            if (cambia > 0):
                tramo = np.abs(c-a)
                a = c
                fa = fc
            else:
                tramo = np.abs(b-c)
                b = c
                fb = fc

            if vertabla==True:
                print('  ',tramo)
            itera = itera + 1
        respuesta = c
        # Valida respuesta
        if (itera>=iteramax):
            respuesta = np.nan
    else: 
        print(' No existe cambio de signo entre f(a) y f(b)')
        print(' f(a) =',fa,',  f(b) =',fb) 
        respuesta=np.nan

    return(respuesta)

# PROGRAMA ----------------------
# INGRESO
fx  = lambda x: x**3 + 4*x**2 - 10
a = 1
b = 2
tolera = 0.0001

# PROCEDIMIENTO
respuesta = posicionfalsa(fx,a,b,tolera,vertabla=True)
# SALIDA
print('raíz en: ', respuesta)

[ Falsa Posición ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ][gráfica]
..


5. Gráfica en Python

La gráfica se puede obtener añadiendo las siguientes instrucciones:

# GRAFICA
import matplotlib.pyplot as plt
muestras = 21

xi = np.linspace(a,b,muestras)
fi = fx(xi)
plt.plot(xi,fi, label='f(x)')
plt.axhline(0)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('Falsa Posición')
plt.grid()
plt.legend()
plt.show()

[ Falsa Posición ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ][gráfica]

2.2 Método de la Falsa Posición - Concepto

[ Falsa posición ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ]
..


Método de la Falsa Posición

Referencia: Burden p56, Chapra 5.3 p131

El método de la posición falsa, falsa posición, regla falsa o regula falsi, considera dividir el intervalo cerrado [a,b] donde se encontraría una raíz de la función f(x) basado en la cercanía a cero que tenga f(a) o f(b).

Posicion Falsa animado GIF

El método une f(a) con f(b) con una línea recta, la intersección de la recta con el eje x representaría una mejor aproximación hacia la raíz.

Al reemplazar la curva de f(x) por una línea recta, se genera el nombre de "posición falsa" de la raíz. El método también se conoce como interpolación lineal.

A partir de la gráfica, usando triángulos semejantes, considerando que f(a) es negativo en el ejemplo, se estima que:

Falsa Posicion animado

\frac{0-f(a)}{c-a} = \frac{f(b)-0}{b-c} \frac{f(a)}{c-a} = -\frac{f(b)}{b-c} \frac{f(a)}{c-a} = \frac{f(b)}{c-b}

que al despejar c, se obtiene:

c = b - f(b) \frac{(a-b)}{f(a)-f(b)}

Calculado el valor de c, éste reemplaza a uno de los valores iniciales [a,b], cuyo valor evaluado tenga el mismo signo que f(c)

Nota: La forma de la expresión presentada para c, se usa para comparar con el método de la secante. Se obtiene sumando y restando b y reagrupando.

Control de iteraciones

Las correcciones del intervalo que se realizan en cada iteración tienen a ser más pequeñas, por lo que el control de iteraciones se realizan sobre la porción o tramo que se redujo el intervalo.

Si la reducción del intervalo es por la izquierda, tramo = c - a
Si la reducción del intervalo es por la derecha, tramo = b - c

[ Falsa posición ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ]

 

2.1.1 Método de la Bisección - Ejemplo con Python

[ Bisección ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ] [ gráfica ]
..


1. Ejercicio

Referencia: Burden 2.1 ejemplo 1 p38

La ecuación mostrada tiene una raíz en [1,2], ya que f(1)=-5 y f(2)=14 y existe cambio de signo. Muestre los resultados parciales del algoritmo de la bisección con una tolerancia de 0.0001

f(x) = x^3 + 4x^2 -10 =0

método de Biseccion animado

[ Bisección ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ] [ gráfica ]
..


2. Desarrollo Analítico

El desarrollo del ejercicio tradicionalmente realizado con lápiz, papel y calculadora, muestra el orden y detalle de las operaciones que se pueden traducir a un algoritmo en Python. El objetivo además de desarrollar la comprensión del método, permite en una evaluación observar si el estudiante conoce el método y usa apropiadamente los valores en cada iteración.

iteración 1

a = 1, b=2 c = \frac{a+b}{2} = \frac{1+2}{2} = 1.5 f(1) = (1)^3 + 4(1)^2 -10 = -5 f(1.5) = (1.5)^3 + 4(1.5)^2 -10= 2.37 f(2) = (2)^3 + 4(2)^2 -10 =14

cambio de signo a la izquierda

a = 1, b= c = 1.5 tramo = |1.5-1| =0.5

iteración 2

a = 1, b=1.5 c = \frac{1+1.5}{2} = 1.25 f(1) = -5 f(1.25) = (1.25)^3 + 4(1.25)^2 -10 = -1.794 f(1.5) = 2.37

cambio de signo a la derecha

a = c = 1.25, b=1.5 tramo = |1.5-1.25| = 0.25

iteración 3

continuar como tarea.

La tabla resume los valores de las iteraciones

tabla para Bisección
i a c b f(a) f(c) f(b) tramo
1 1 1.5 2 -5 2.37 14 0.5
2 1 1.25 1.5 -5 -1.79 2.37 0.25
3 1.25 ... 1.5

La misma tabla se puede realizar con un algoritmo para tener los resultados más rápido y observar el comportamiento del método.

Observe los resultados de f(c), principalmente en la iteración i=9 con tramo=0.00097 que representa el error de estimación del valor vs tolerancia.

i ['a', 'c', 'b'] ['f(a)', 'f(c)', 'f(b)']
   tramo
0 [1.  1.5 2. ] [-5.     2.375 14.   ]
   0.5
1 [1.   1.25 1.5 ] [-5.     -1.7969  2.375 ]
   0.25
2 [1.25  1.375 1.5  ] [-1.7969  0.1621  2.375 ]
   0.125
3 [1.25   1.3125 1.375 ] [-1.7969 -0.8484  0.1621]
   0.0625
4 [1.3125 1.3438 1.375 ] [-0.8484 -0.351   0.1621]
   0.03125
5 [1.3438 1.3594 1.375 ] [-0.351  -0.0964  0.1621]
   0.015625
6 [1.3594 1.3672 1.375 ] [-0.0964  0.0324  0.1621]
   0.0078125
7 [1.3594 1.3633 1.3672] [-0.0964 -0.0321  0.0324]
   0.00390625
8 [1.3633 1.3652 1.3672] [-3.2150e-02  7.2025e-05  3.2356e-02]
   0.001953125
9 [1.3633 1.3643 1.3652] [-3.2150e-02 -1.6047e-02  7.2025e-05]
   0.0009765625
raíz en:  1.3642578125
>>> 

Se realiza la gráfica los puntos [c,f(c)] de la tabla para observar el resultado, resaltando que los puntos al final se aglomeran alrededor de la solución o raíz de la ecuación.

método de la bisección

[ Bisección ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ] [ gráfica ]
..


3. Algoritmo en Python

El video presenta el desarrollo básico conceptual del algoritmo en Python para una comprensión del proceso paso a paso.

Instrucciones en Python del Algoritmo básico del video

# Algoritmo de Bisección
# [a,b] se escogen de la gráfica de la función
# error = tolera

import numpy as np
import matplotlib.pyplot as plt

# INGRESO
fx = lambda x: x**3 + 4*x**2 - 10 
a = 1
b = 2
tolera = 0.001

# PROCEDIMIENTO
tramo = b-a
while not(tramo<tolera):
    c = (a+b)/2
    fa = fx(a)
    fb = fx(b)
    fc = fx(c)
    cambia = np.sign(fa)*np.sign(fc)
    if cambia < 0: 
        a = a
        b = c
    if cambia > 0:
        a = c
        b = b
    tramo = b-a

# SALIDA
print('       raiz en: ', c)
print('error en tramo: ', tramo)

[ Bisección ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ] [ gráfica ]

..


4. Algoritmo en Python como función

El algoritmo presentado en el video se puede mejorar, por ejemplo simplificando los dos condicionales en uno.

            if (cambia<0):
                b = c
                fb = fc
            else:
                a = c
                fa = fc

Considere que en cada iteración, se evalúa la función en tres puntos, fa y fb se calculan en la iteración anterior. Se puede optimizar sustituyendo los valores de los extremos y solo evaluando el centro fc.

    fa = fx(a)
    fb = fx(b)
    tramo = np.abs(b-a)
    itera = 0
...
        while (tramo>=tolera and itera<=iteramax):
            c = (a+b)/2
            fc = fx(c)
            cambia = np.sign(fa)*np.sign(fc)

Para el desarrollo con lápiz y papel sería necesario observar los valores calculados en cada iteración, por lo que se añade una variable "vertabla" para activar las instrucciones que muestran los cálculos en cada iteración.

        if vertabla==True:
            print('método de Bisección')
            print('i', ['a','c','b'],[ 'f(a)', 'f(c)','f(b)'])
            print('  ','tramo')
            np.set_printoptions(precision)

Se incluye validar que el intervalo disponga de fa y fb con signos opuestos, antes de intentar usar el algoritmo

    if cambia<0: # existe cambio de signo f(a) vs f(b)

Finalmente se puede convertir el procedimiento en una función de Python.

Algoritmo en Python como función

# Algoritmo de Bisección
# [a,b] se escogen de la gráfica de la función
# error = tolera
import numpy as np

def biseccion(fx,a,b,tolera,iteramax = 50, vertabla=False, precision=4):
    '''
    Algoritmo de Bisección
    Los valores de [a,b] son seleccionados
    desde la gráfica de la función
    error = tolera
    '''
    fa = fx(a)
    fb = fx(b)
    tramo = np.abs(b-a)
    itera = 0
    cambia = np.sign(fa)*np.sign(fb)
    if cambia<0: # existe cambio de signo f(a) vs f(b)
        if vertabla==True:
            print('método de Bisección')
            print('i', ['a','c','b'],[ 'f(a)', 'f(c)','f(b)'])
            print('  ','tramo')
            np.set_printoptions(precision)
            
        while (tramo>=tolera and itera<=iteramax):
            c = (a+b)/2
            fc = fx(c)
            cambia = np.sign(fa)*np.sign(fc)
            if vertabla==True:
                print(itera,np.array([a,c,b]),
                      np.array([fa,fc,fb]))
            if (cambia<0):
                b = c
                fb = fc
            else:
                a = c
                fa = fc
            tramo = np.abs(b-a)
            if vertabla==True:
                print('  ',tramo)
            itera = itera + 1
        respuesta = c
        # Valida respuesta
        if (itera>=iteramax):
            respuesta = np.nan

    else: 
        print(' No existe cambio de signo entre f(a) y f(b)')
        print(' f(a) =',fa,',  f(b) =',fb) 
        respuesta=np.nan
    return(respuesta)

# PROGRAMA ----------------------
# INGRESO
fx  = lambda x: x**3 + 4*x**2 - 10
a = 1
b = 2
tolera = 0.001

# PROCEDIMIENTO
respuesta = biseccion(fx,a,b,tolera,vertabla=True)
# SALIDA
print('raíz en: ', respuesta)

[ Bisección ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ] [ gráfica ]
..


5. Gráfica en Python

La gráfica se puede obtener añadiendo las siguientes instrucciones:

# GRAFICA
import matplotlib.pyplot as plt
muestras = 21

xi = np.linspace(a,b,muestras)
fi = fx(xi)
plt.plot(xi,fi, label='f(x)')
plt.axhline(0)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('Biseccion')
plt.grid()
plt.legend()
plt.show()

[ Bisección ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ] [ gráfica ]

..


6. Función en librería Scipy.optimize.bisect

El método de la bisección se encuentra también implementado en las librería Scipy, que también puede ser usado de la forma:

>>> import scipy as sp
>>> fx = lambda x: x**3 + 4*x**2 - 10
>>> sp.optimize.bisect(fx,1,2,xtol=0.001)
1.3642578125

que es el resultado de la raíz para la última iteración del ejercicio. Lo que muestra que el algoritmo realizado tiene un valor más aproximado.

Sin embargo por didáctica y mejor comprensión de los métodos y su implementación en algoritmos que es parte del objetivo de aprendizaje, se continuará desarrollando la forma básica y detallada con Python.

Referencia: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.bisect.html

[ Bisección ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ] [ gráfica ]

2.1 Método de la Bisección - Concepto

[ Bisección ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ]

..


Método de la Bisección

Referencia: Burden 2.1 p36, Chapra 5.2 p124, Rodríguez 3.1 p36

El método más simple para buscar raíces de una ecuación, se basa en el teorema del valor intermedio, búsqueda binaria, partición de intervalos o de Bolzano.

método de Biseccion animado

En el intervalo donde existe un cruce por cero de la función f(x), el algoritmo busca la raíz al reducir el intervalo en la mitad (bisección), seleccionando el sub- intervalo donde se mantenga el cambio de signo de la función f(x).

Los pasos a seguir son los siguientes:

  • el intervalo [a,b] se divide siempre en la mitad c.
  • Si la función f(x) cambia de signo sobre un intervalo, se evalúa el valor de la función en el punto medio f(c).
  • La posición de la raíz se determina en el punto medio del sub-intervalo, izquierdo o derecho,  dentro del cual ocurre un "cambio de signo".
  • el proceso se repite hasta obtener una mejor aproximación

La gráfica muestra una animación del proceso, observe la forma en que progresivamente se acercan los puntos [a,b], donde se mantienen valores con signo diferente entre f(a) y f(b).

Para describir mejor el método, observamos la gráfica en una sola iteración.
Para la primera iteración se tiene que la función tiene un cambio de signo dentro del intervalo [a,b].biseccion 01

El intervalo se divide en la mitad, representado por el punto c, obteniendo el sub-intervalo izquierdo [a,c] o sub-intervalo derecho [c,b].

El sub-intervalo que contiene la función con un cambio de signo, se convierte en el nuevo intervalo por analizar en la siguiente iteración.

[ Bisección ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ]


Cota de Error

Referencia: Burden Teorema 2.1  p39.

El error del método de la bisección se estima como el ancho o tamaño del intervalo [a,b] de la última iteración realizada. Si el error es menor que la tolerancia del ejercicio, el algoritmo se detiene y se considera encontrada la raíz.

Suponga que f ∈ C[a,b] y f(a)*f(b)<0, f es una función en el intervalo [a,b] y que presenta un cambio de signo.

|p_n - p| \leq \frac{b-a}{2^n} \text{donde } n \geq 1

la desigualdad implica que pn converge a p con una razón de convergencia de orden:

O \Big(\frac{1}{2^n}\Big)

es decir:

p_n =p+O \Big( \frac{1}{2^n} \Big)

Con lo que se puede determinar el número de iteraciones necesarias para encontrar la raíz, tal como se muestra en el siguiente ejercicio.

[ Bisección ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ]


Cantidad de iteraciones - Ejercicio

Referencia: Burden ejemplo 2  p40

Determine la cantidad de iteraciones necesarias para resolver

f(x) = x^3 + 4x^2 -10 =0

con exactitud de 10 - 3 en el intervalo [1,2].

Desarrollo: Se busca encontrar un entero n que satisface la ecuación:

|p_n -p| \leq \frac{b-a}{2^{n}} 2^{-n}< 10^{-3}

usando logaritmos:

-n \log _{10}( 2) < -3 n > \frac{3}{\log _{10}( 2)} = 9.96

En consecuencia se requieren unas diez iteraciones para lograr la aproximación de 10-3. Verifique los resultados con los valores calculados.

[ Bisección ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ]