2.4.1 Método del Punto fijo – Ejemplo con Python

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

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

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

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.

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).

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

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

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

  • iniciando desde el extremo izquierdo x=a,
  • se determina el valor de b = g(x) ,
  • se determina la diferencia de la aproximación o error = |b-a|, tolera=0.001
  • se proyecta en la recta identidad, como el nuevo punto de evaluación
    x=b.
  • se repite el proceso para el nuevo valor de x, 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 1

x_0 = 0 x_1 = g(0) = e^{-0} = 1 tramo =|1-0|= 1

iteración 2

x_1 = 1 x_2 = g(1) = e^{-1} = 0.3678 tramo =|0.3678 - 1|= 0.6322

iteración 3

x_2 = 0.3678 x_3 = 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 x xnuevo = g(x) |error|
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.


Algoritmo en Python

El algoritmo en python se presenta como una función para usarla fácilmente como un bloque en otros ejercicios. 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 convergéncia de la función con el método.

Se presenta el algoritmo mejorado, convirtiendo el procedimiento del método a una función Python.

# Algoritmo de punto fijo
# [a,b] intervalo de búsqueda
# error = tolera

import numpy as np

def puntofijo(gx,a,tolera, iteramax = 15):
    i = 1 # iteración
    b = gx(a)
    tramo = abs(b-a)
    while(tramo>=tolera and i<=iteramax ):
        a = b
        b = gx(a)
        tramo = abs(b-a)
        i = i + 1
    respuesta = b
    
    # Validar respuesta
    if (i>=iteramax ):
        respuesta = np.nan
    return(respuesta)

# PROGRAMA ---------

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

a = 0       # intervalo
b = 1
tolera = 0.001
iteramax = 15  # itera máximo
muestras = 51  # gráfico
tramos = 50

# PROCEDIMIENTO
respuesta = puntofijo(gx,a,tolera)

# SALIDA
print(respuesta)

la respuesta obtenida del problema es

0.566908911921

Para obtener la gráfica básica se determinan los puntos para cada función fx y gx. Se añaden las siguiente instrucciones al algoritmo anterior:

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

import matplotlib.pyplot as plt

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

if (respuesta!= np.nan):
    plt.axvline(respuesta)
plt.axhline(0, color='k')
plt.title('Punto Fijo')
plt.legend()
plt.show()

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 scipy.optimize as opt
>>> opt.fixed_point(gx,a,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 raiz.
  • 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.
  • Realizar el mismo ejercicio para
f(x) = x^{2}-x-2 = 0

en el rango [-1,2]
Realice sus observaciones y recomendaciones.

2.3.1 Método de Newton-Raphson Ejemplo con Python

Referencia: Burden ejemplo 1 p51

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

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


Algoritmo con Python

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

De la sección anterior, luego de realizar tres iteraciones para la tabla, notamos la necesidad de usar un algoritmo para que realice los cálculos repetitivos y muestre la tabla o directamente el resultado.

 ['xi', 'xnuevo', 'tramo']
[[2.0000 1.5000 5.0000e-01]
 [1.5000 1.3733 1.2667e-01]
 [1.3733 1.3653 8.0713e-03]
 [1.3653 1.3652 3.2001e-05]]
raiz en:  1.3652300139161466
con error de:  3.200095847999407e-05

Al algoritmo básico se les añade lo necesario para mostrar la tabla con los valores de las iteraciones.

# Método de Newton-Raphson
# Ejemplo 1 (Burden ejemplo 1 p.51/pdf.61)

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

# PROCEDIMIENTO
tabla = []
tramo = abs(2*tolera)
xi = x0
while (tramo>=tolera):
    xnuevo = xi - fx(xi)/dfx(xi)
    tramo  = abs(xnuevo-xi)
    tabla.append([xi,xnuevo,tramo])
    xi = xnuevo

# convierte la lista a un arreglo.
tabla = np.array(tabla)
n = len(tabla)

# SALIDA
print(['xi', 'xnuevo', 'tramo'])
np.set_printoptions(precision = 4)
print(tabla)
print('raiz en: ', xi)
print('con error de: ',tramo)

scipy.optimize.newton

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

>>> import scipy.optimize as opt
>>> opt.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

  • convertir el algoritmo a una función
  • Revisar las modificaciones si se quiere usar la forma simbólica de la función.
  • incorpore la grafica básica 2D de la función f(x)

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

Referencia: Burden 9Ed ejemplo 1 p51

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

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

posicionfalsa01_GIF


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 raiz.

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 contrinuar con las iteraciones como tarea


Algoritmo en Python


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)
    tabla.append([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
        
tabla = np.array(tabla)
ntabla = len(tabla)

# SALIDA
np.set_printoptions(precision=4)
for i in range(0,ntabla,1):
    print('iteración:  ',i)
    print('[a,c,b]:    ', tabla[i,0:3])
    print('[fa,fc,fb]: ', tabla[i,3:6])
    print('[tramo]:    ', tabla[i,6])

print('raiz:  ',c)
print('error: ',tramo)

Observe el número de iteraciones realizadas, hasta presentar el valor de la raiz en 1.3652 con un error de 0.0003166860575976038 en la útlima 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.

iteración:   0
[a,c,b]:     [1.     1.2632 2.    ]
[fa,fc,fb]:  [-5.     -1.6023 14.    ]
[tramo]:     1.0
iteración:   1
[a,c,b]:     [1.2632 1.3388 2.    ]
[fa,fc,fb]:  [-1.6023 -0.4304 14.    ]
[tramo]:     0.26315789473684204
iteración:   2
[a,c,b]:     [1.3388 1.3585 2.    ]
[fa,fc,fb]:  [-0.4304 -0.11   14.    ]
[tramo]:     0.0756699440909967
iteración:   3
[a,c,b]:     [1.3585 1.3635 2.    ]
[fa,fc,fb]:  [-0.11   -0.0278 14.    ]
[tramo]:     0.019718502996940224
iteración:   4
[a,c,b]:     [1.3635 1.3648 2.    ]
[fa,fc,fb]:  [-2.7762e-02 -6.9834e-03  1.4000e+01]
[tramo]:     0.005001098217311428
iteración:   5
[a,c,b]:     [1.3648 1.3651 2.    ]
[fa,fc,fb]:  [-6.9834e-03 -1.7552e-03  1.4000e+01]
[tramo]:     0.0012595917846898175
iteración:   6
[a,c,b]:     [1.3651 1.3652 2.    ]
[fa,fc,fb]:  [-1.7552e-03 -4.4106e-04  1.4000e+01]
[tramo]:     0.0003166860575976038
raiz:   1.3652033036626001
error:  7.958577822231305e-05

Ejemplo 2

Referencia: Burden ejemplo 3 p.74/pdf.84

Con el método de la posición falsa, determine la cantidad de iteraciones necesarias para resolver:

f(x) = \cos (x) - x

Realice una tabla que muestre: [ i, a, b, c, f(c) ].
Seleccióne el rango [a,b] usando una gráfica y compruebe los resultados.
Compare con otros métodos para encontrar raíces.

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

Referencia: Burden 9Ed. ejemplo 1 p50

La ecuación mostrada tiene una raiz 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


Desarrollo Analítico

Como parte del desarrollo del ejercicio se presenta las iteraciones para el algoritmo, tradicionalmente realizadas con una calculadora.

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 tramo = 2-1 =1

cambio de signo a la izquierda

a = 1, b= c = 1.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 tramo = 1.5-1 = 0.5

cambio de signo a la derecha

a = c = 1.25, b=1.5

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 1
2 1 1.25 1.5 -5 -1.794 2.37 0.5
3 1.25 1.5

La misma tabla se puede relizar 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 iteracion i=9 con tramo=0.004, respecto a la obtenida en la última iteración.

[i, a,   c,    b,    f(a),   f(c), f(b),  tramo]
 1 1.000 1.500 2.000 -5.000  2.375 14.000 1.000 
 2 1.000 1.250 1.500 -5.000 -1.797  2.375 0.500 
 3 1.250 1.375 1.500 -1.797  0.162  2.375 0.250 
 4 1.250 1.312 1.375 -1.797 -0.848  0.162 0.125 
 5 1.312 1.344 1.375 -0.848 -0.351  0.162 0.062 
 6 1.344 1.359 1.375 -0.351 -0.096  0.162 0.031 
 7 1.359 1.367 1.375 -0.096  0.032  0.162 0.016 
 8 1.359 1.363 1.367 -0.096 -0.032  0.032 0.008 
 9 1.363 1.365 1.367 -0.032  0.000  0.032 0.004 
10 1.363 1.364 1.365 -0.032 -0.016  0.000 0.002 
11 1.364 1.365 1.365 -0.016 -0.008  0.000 0.001 
raiz:  1.36474609375
>>> 

Se grafican los puntos [c,f(c)]  de la tabla para obsevar el resultado de forma gráfica, resaltando que los puntos al final se aglomeran alrededor de la solución o raiz de la ecuación.

Escriba sus observaciones y preguntas sobre los resultados.


Algoritmo en Python

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

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)

Mejorando el algoritmo

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

Considere que siempre se evalúa la función en tres puntos y se puede optimizar sustituyendo los valores de los extremos  y solo evaluando el centro.

El resultado se mejorar usando los valores finales del último intervalo [a,b] y obteniendo una nueva midad ‘c’ al final.

Para mostrar la tabla de la sección anterior se agregan los resultados parciales a una tabla que permita mostrar el resultado al final

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

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

# PROCEDIMIENTO
tabla = []
tramo = b-a

fa = fx(a)
fb = fx(b)
i = 1
while (tramo>tolera):
    c = (a+b)/2
    fc = fx(c)
    tabla.append([i,a,c,b,fa,fc,fb,tramo])
    i = i + 1
                 
    cambia = np.sign(fa)*np.sign(fc)
    if (cambia<0):
        b = c
        fb = fc
    else:
        a=c
        fa = fc
    tramo = b-a
c = (a+b)/2
fc = fx(c)
tabla.append([i,a,c,b,fa,fc,fb,tramo])
tabla = np.array(tabla)

raiz = c

# SALIDA
np.set_printoptions(precision = 4)
print('[ i, a, c, b, f(a), f(c), f(b), tramo]')
# print(tabla)

# Tabla con formato
n=len(tabla)
for i in range(0,n,1):
    unafila = tabla[i]
    formato = '{:.0f}'+' '+(len(unafila)-1)*'{:.3f} '
    unafila = formato.format(*unafila)
    print(unafila)
    
print('raiz: ',raiz)

El ultimo complemento al algoritmo consiste en realizar la gráfica, seleccionando solo las columnas correspondientes a [c,f(c)].

Se ordenan los datos en forma ascendente antes de graficarlos usando solo sus índices con np.argsort(xi)

# Algoritmo de Bisección
# GRAFICA
import matplotlib.pyplot as plt

xi = tabla[:,2]
yi = tabla[:,5]

# ordena los puntos para la grafica
orden = np.argsort(xi)
xi = xi[orden]
yi = yi[orden]

plt.plot(xi,yi)
plt.plot(xi,yi,'o')
plt.axhline(0, color="black")

plt.xlabel('x')
plt.ylabel('y')
plt.title('Bisección en f(x)')
plt.grid()
plt.show()

Scipy.optimize.bisect

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

>>> import scipy.optimize as opt
>>> opt.bisect(fx,1,2,xtol=0.001)
1.3642578125

que es el valor de la variable ‘a’ de la tabla 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 en Python.

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