Categoría: Unidades

Unidades de estudio

  • 2.7 Métodos de raíces con gráficos animados en Python

    GIF animado: [ Bisección ] [ Posicion Falsa ] [ Newton-Raphson ] [ Punto Fijo ] [ Secante ]

    Solo para fines didácticos, y como complemento para los ejercicios presentados en la unidad para raíces de ecuaciones, se presentan las instrucciones para las animaciones usadas en la presentación de los conceptos y ejercicios. Los algoritmos para animación NO son necesarios para realizar los ejercicios, que requieren una parte analítica con al menos tres iteraciones en papel y lápiz. Se lo adjunta como una herramienta didáctica de asistencia para las clases.

    La gráfica (graf_ani) se crea en una ventana (fig_ani), inicializando con la linea a partir f(x) y configurando los parámetros base para el gráfico.

    Se usan procedimientos para crear unatrama() para marca de iteración y en cada cambio se limpia la trama manteniendo la base con limpiatrama().

    En caso de requerir un archivo .gif animado se proporciona un nombre de archivo. Para crear el archivo se requiere de la librería 'pillow', a ser instalado.

    otros ejemplos de animación en el curso de Fundamentos de Programación:

    Movimiento circular – Una partícula, animación con matplotlib-Python

    GIF animado: [ Bisección ] [ Posicion Falsa ] [ Newton-Raphson ] [ Punto Fijo ] [ Secante ]
    ..


    Método de la Bisección con gráfico animado en Python

    método de Biseccion animado

    # Algoritmo de Bisecci n, Tabla
    # Los valores de [a,b] son aceptables
    # y seleccionados desde la gr fica de la funci n
    # error = tolera
    
    import numpy as np
    
    def biseccion_tabla(fx,a,b,tolera,iteramax = 20,
                        vertabla=False, precision=6):
        '''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)
        tabla=[]
        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)
                unafila = [a,c,b,fa,fc,fb] 
                if vertabla==True:
                    print(itera,[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)
                unafila.append(tramo)
                tabla.append(unafila)
                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
        tabla = np.array(tabla,dtype=float)
        return(respuesta,tabla)
    
    # PROGRAMA #######################
    
    # INGRESO
    fx = lambda x: x**3 + 4*x**2 - 10
    
    a = 1
    b = 2
    tolera = 0.001
    
    # PROCEDIMIENTO
    [raiz,tabla] = biseccion_tabla(fx,a,b,tolera,vertabla=True)
    
    # SALIDA
    print('ra z en: ', raiz)
    
    # GRAFICA
    import matplotlib.pyplot as plt
    muestras = 21
    xi = np.linspace(a,b,muestras)
    fi = fx(xi)
    xc = tabla[:,1]
    yc = tabla[:,4]
    
    plt.plot(xi,fi, label='f(x)')
    plt.plot([a,b],[fx(a),fx(b)],'o',
             color='red',label='[[a,b],[f(a),f(b)]]')
    plt.scatter(xc,yc,color='orange', label='[c,f(c)]')
    plt.axhline(0)
    plt.xlabel('x')
    plt.ylabel('f(x)')
    plt.title('Bisecci n')
    plt.grid()
    plt.legend()
    #plt.show()
    
    # GRAFICA CON ANIMACION ------------
    import matplotlib.animation as animation
    
    xa = tabla[:,0]
    ya = tabla[:,3]
    xb = tabla[:,2]
    yb = tabla[:,5]
    xc = tabla[:,1]
    yc = tabla[:,4]
    
    # Inicializa parametros de trama/foto
    narchivo = 'Biseccion' # nombre archivo
    retardo = 700 # milisegundos entre tramas
    tramas = len(xa)
    
    # GRAFICA animada en fig_ani
    fig_ani, graf_ani = plt.subplots()
    ymax = np.max(fi)
    ymin = np.min(fi)
    deltax = np.abs(b-a)
    deltay = np.abs(ymax-ymin)
    graf_ani.set_xlim([a-0.05*deltax,b+0.05*deltax])
    graf_ani.set_ylim([ymin-0.05*deltay,ymax+0.05*deltay])
    linea0 = graf_ani.axhline(0, color='k')
    # Lineas y puntos base
    lineafx,= graf_ani.plot(xi,fi, label='f(x)')
    
    puntoa, = graf_ani.plot(xa[0], ya[0],'o',
                            color='red', label='a')
    puntob, = graf_ani.plot(xb[0], yb[0],'o',
                            color='green', label='b')
    puntoc, = graf_ani.plot(xc[0], yc[0],'o',
                            color='orange', label='c')
    lineaa, = graf_ani.plot([xa[0],xa[0]],
                            [0,ya[0]],color='red',
                            linestyle='dashed')
    lineab, = graf_ani.plot([xb[0],xb[0]],
                            [0,yb[0]],color='green',
                            linestyle='dashed')
    lineac, = graf_ani.plot([xc[0],xc[0]],
                            [0,yc[0]],
                            color='orange',
                            linestyle='dashed')
    linea_ab, = graf_ani.plot([xa[0],xb[0]],
                            [0,0],
                            color='yellow',
                            linestyle='dotted')
    # Configura gr fica
    graf_ani.set_title('Bisecci n')
    graf_ani.set_xlabel('x')
    graf_ani.set_ylabel('f(x)')
    graf_ani.legend()
    graf_ani.grid()
    
    # Cada nueva trama
    def unatrama(i,xa,ya,xb,yb,xc,yc):
        # actualiza cada punto
        puntoa.set_xdata([xa[i]]) 
        puntoa.set_ydata([ya[i]])
        puntob.set_xdata([xb[i]])
        puntob.set_ydata([yb[i]])
        puntoc.set_xdata([xc[i]])
        puntoc.set_ydata([yc[i]])
        # actualiza cada linea
        lineaa.set_ydata([ya[i], 0])
        lineaa.set_xdata([xa[i], xa[i]])
        lineab.set_ydata([yb[i], 0])
        lineab.set_xdata([xb[i], xb[i]])
        lineac.set_ydata([yc[i], 0])
        lineac.set_xdata([xc[i], xc[i]])
        linea_ab.set_ydata([0, 0])
        linea_ab.set_xdata([xa[i], xb[i]])
        return (puntoa, puntob, puntoc, lineaa, lineab, lineac,linea_ab)
    
    # Limpia trama anterior
    def limpiatrama():
        puntoa.set_ydata(np.ma.array(xa, mask=True))
        puntob.set_ydata(np.ma.array(xb, mask=True))
        puntoc.set_ydata(np.ma.array(xc, mask=True))
        lineaa.set_ydata(np.ma.array([0,0], mask=True))
        lineab.set_ydata(np.ma.array([0,0], mask=True))
        lineac.set_ydata(np.ma.array([0,0], mask=True))
        linea_ab.set_ydata(np.ma.array([0,0], mask=True))
        return (puntoa, puntob, puntoc, lineaa, lineab, lineac,linea_ab)
    
    # contador de tramas
    i = np.arange(0,tramas,1)
    ani = animation.FuncAnimation(fig_ani,unatrama,
                                  i ,
                                  fargs=(xa, ya,
                                         xb, yb,
                                         xc, yc),
                                  init_func=limpiatrama,
                                  interval=retardo,
                                  blit=True)
    # Graba Archivo GIFAnimado y video
    ani.save(narchivo+'_animado.gif', writer='pillow')
    #ani.save(narchivo+'_animado.mp4')
    plt.show()
    
    

    GIF animado: [ Bisección ] [ Posicion Falsa ] [ Newton-Raphson ] [ Punto Fijo ] [ Secante ]

    ..


    Método de la Posición Falsa con gráfico animado en Python

    Posicion Falsa animado GIF

    # 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_tabla(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)
        tabla = []
        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)
                unafila = [a,c,b,fa, fc, fb]
                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
                unafila.append(tramo)
                tabla.append(unafila)
                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
        tabla = np.array(tabla,dtype=float)
        return(respuesta,tabla)
    
    # PROGRAMA ----------------------
    
    # INGRESO
    fx = lambda x: x**3 + 4*x**2 - 10
    
    a = 1
    b = 2
    tolera = 0.001
    
    # PROCEDIMIENTO
    [raiz,tabla] = posicionfalsa_tabla(fx,a,b,tolera,
                                       vertabla=True)
    # SALIDA
    print('raíz en: ', raiz)
    
    # GRAFICA
    import matplotlib.pyplot as plt
    muestras = 21
    xi = np.linspace(a,b,muestras)
    fi = fx(xi)
    xc = tabla[:,1]
    yc = tabla[:,4]
    
    plt.plot(xi,fi, label='f(x)')
    plt.plot([a,b],[fx(a),fx(b)],'o',
             color='red',label='[[a,b],[f(a),f(b)]]')
    plt.scatter(xc,yc,color='orange', label='[c,f(c)]')
    plt.axhline(0)
    plt.xlabel('x')
    plt.ylabel('f(x)')
    plt.title('Posición Falsa')
    plt.grid()
    plt.legend()
    #plt.show()
    
    # GRAFICA CON ANIMACION ------------
    #import matplotlib.pyplot as plt
    import matplotlib.animation as animation
    
    xa = tabla[:,0]
    ya = tabla[:,3]
    xc = tabla[:,1]
    yc = tabla[:,4]
    xb = tabla[:,2]
    yb = tabla[:,5]
    
    # Inicializa parametros de trama/foto
    narchivo = 'PosicionFalsa' # nombre archivo
    retardo = 700 # milisegundos entre tramas
    tramas = len(xa)
    
    # GRAFICA animada en fig_ani
    fig_ani, graf_ani = plt.subplots()
    graf_ani.set_xlim([a,b])
    graf_ani.set_ylim([np.min(fi),np.max(fi)])
    # Lineas y puntos base
    lineafx, = graf_ani.plot(xi,fi,label ='f(x)')
    
    puntoa, = graf_ani.plot(xa[0], ya[0],'o',
                            color='red', label='a')
    puntob, = graf_ani.plot(xb[0], yb[0],'o',
                            color='green', label='b')
    puntoc, = graf_ani.plot(xc[0], yc[0],'o',
                            color='orange', label='c')
    
    lineaab, = graf_ani.plot([xa[0],xb[0]],[ya[0],yb[0]],
                            color ='orange',
                            label='y=mx+b')
    lineac0, = graf_ani.plot([xc[0],xc[0]],[0,yc[0]],
                            color='magenta',
                            linestyle='dashed')
    
    # Configura gráfica
    linea0 = graf_ani.axhline(0, color='k')
    graf_ani.set_title('Posición Falsa')
    graf_ani.set_xlabel('x')
    graf_ani.set_ylabel('f(x)')
    graf_ani.legend()
    graf_ani.grid()
    
    # Cada nueva trama
    def unatrama(i,xa,ya,xc,yc,xb,yb):
        # actualiza cada punto, [] porque es solo un valor
        puntoa.set_xdata([xa[i]])
        puntoa.set_ydata([ya[i]])
        puntob.set_xdata([xb[i]])
        puntob.set_ydata([yb[i]])
        puntoc.set_xdata([xc[i]])
        puntoc.set_ydata([yc[i]])
        # actualiza cada linea
        lineaab.set_ydata([ya[i], yb[i]])
        lineaab.set_xdata([xa[i], xb[i]])
        lineac0.set_ydata([0, yc[i]])
        lineac0.set_xdata([xc[i], xc[i]])  
        return (puntoa, puntob, puntoc, lineaab,lineac0,)
    
    # Cada nueva trama
    def limpiatrama():
        puntoa.set_ydata(np.ma.array(xa, mask=True))
        puntob.set_ydata(np.ma.array(xb, mask=True))
        puntoc.set_ydata(np.ma.array(xc, mask=True))
        lineaab.set_ydata(np.ma.array([0,0], mask=True))
        lineac0.set_ydata(np.ma.array([0,0], mask=True))
        return (puntoa, puntob, puntoc, lineaab,lineac0,)
    
    # contador de tramas
    i = np.arange(0, tramas,1)
    ani = animation.FuncAnimation(fig_ani,unatrama,
                                  i ,
                                  fargs=(xa, ya,
                                         xc, yc,
                                         xb,yb),
                                  init_func=limpiatrama,
                                  interval=retardo,
                                  blit=True)
    # Graba Archivo GIFAnimado y video
    ani.save(narchivo+'_animado.gif', writer='pillow')
    #ani.save(narchivo+'.mp4')
    plt.show()
    

    GIF animado: [ Bisección ] [ Posicion Falsa ] [ Newton-Raphson ] [ Punto Fijo ] [ Secante ]
    ..


    Método de Newton-Raphson con gráfico animado en Python

    Newton Raphson GIF animado

    # Método de Newton-Raphson
    # Ejemplo 1 (Burden ejemplo 1 p.51/pdf.61)
    import numpy as np
    
    def newton_raphson_tabla(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)
        tabla = []
        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]))
    
            tabla.append([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')
        tabla = np.array(tabla,dtype=float)
        return(xi,tabla)
    
    # PROGRAMA ----------------------
    
    # INGRESO
    fx = lambda x: x**3 + 4*x**2 - 10
    dfx = lambda x: 3*(x**2) + 8*x
    
    a = 1
    b = 4
    x0 = 3
    tolera = 0.001
    
    # PROCEDIMIENTO
    [raiz,tabla] = newton_raphson_tabla(fx, dfx, x0,
                                 tolera, vertabla=True)
    # SALIDA
    print('raiz en :', raiz)
    
    # GRAFICA
    import matplotlib.pyplot as plt
    muestras = 21
    xi = np.linspace(a,b,muestras)
    fi = fx(xi)
    xc = tabla[:,3]
    yc = fx(xc)
    
    plt.plot(xi,fi, label='f(x)')
    plt.plot(x0,fx(x0),'o',
             color='red',label='[x0,f(x0)]')
    plt.scatter(xc,yc,color='orange', label='[c,f(c)]')
    plt.axhline(0)
    plt.xlabel('x')
    plt.ylabel('f(x)')
    plt.grid()
    plt.legend()
    #plt.show()
    
    # GRAFICA CON ANIMACION ------------
    #import matplotlib.pyplot as plt
    import matplotlib.animation as animation
    
    xa = tabla[:,0]
    ya = tabla[:,1]
    xb = tabla[:,3]
    dfi = tabla[:,2]
    # Aproximacion con tangente
    b0 = ya[0] - dfi[0]*x0
    tangentei = dfi[0]*xi + b0
    ci = -b0/dfi[0]
    
    # Inicializa parametros de trama/foto
    narchivo = 'NewtonRaphson' # nombre archivo
    retardo = 700 # milisegundos entre tramas
    tramas = len(xa)
    
    # GRAFICA animada en fig_ani
    fig_ani, graf_ani = plt.subplots()
    graf_ani.set_xlim([a,b])
    graf_ani.set_ylim([np.min(fi),np.max(fi)])
    # Lineas y puntos base
    lineafx, = graf_ani.plot(xi,fi,label ='f(x)')
    
    puntoa, = graf_ani.plot(xa[0], ya[0],'o',
                            color='red', label='x[i]')
    puntob, = graf_ani.plot(xb[0], 0,'o',
                            color='green', label='x[i+1]')
    
    lineatanx, = graf_ani.plot(xi,tangentei,
                               color='orange',label='tangente')
    lineaa, = graf_ani.plot([xa[0],xa[0]],
                            [ya[0],0], color='magenta',
                            linestyle='dashed')
    lineab, = graf_ani.plot([xb[0],xb[0]],
                             [0,fx(xb[0])], color='magenta',
                             linestyle='dashed')
    # Configura gráfica
    linea0 = graf_ani.axhline(0, color='k')
    graf_ani.set_title('Newton-Raphson')
    graf_ani.set_xlabel('x')
    graf_ani.set_ylabel('f(x)')
    graf_ani.legend()
    graf_ani.grid()
    
    # Cada nueva trama
    def unatrama(i,xa,ya,xb,dfi):
        # actualiza cada punto, [] porque es solo un valor
        puntoa.set_xdata([xa[i]]) 
        puntoa.set_ydata([ya[i]])
        puntob.set_xdata([xb[i]])
        puntob.set_ydata([0])
        # actualiza cada linea
        lineaa.set_ydata([ya[i], 0])
        lineaa.set_xdata([xa[i], xa[i]])
        lineab.set_ydata([0, fx(xb[i])])
        lineab.set_xdata([xb[i], xb[i]])
        # Aproximacion con tangente
        b0 = ya[i] - dfi[i]*xa[i]
        tangentei = dfi[i]*xi+b0
        lineatanx.set_ydata(tangentei)
        
        return (puntoa, puntob, lineaa, lineab,lineatanx,)
    
    # Limpia trama anterior
    def limpiatrama():
        puntoa.set_ydata(np.ma.array(xa, mask=True))
        puntob.set_ydata(np.ma.array(xb, mask=True))
        lineaa.set_ydata(np.ma.array([0,0], mask=True))
        lineab.set_ydata(np.ma.array([0,0], mask=True))
        lineatanx.set_ydata(np.ma.array(xi, mask=True))
        return (puntoa, puntob, lineaa, lineab,lineatanx,)
    
    # contador de tramas
    i = np.arange(0,tramas,1)
    ani = animation.FuncAnimation(fig_ani,unatrama,
                                  i ,
                                  fargs=(xa, ya,
                                         xb, dfi),
                                  init_func=limpiatrama,
                                  interval=retardo,
                                  blit=True)
    # Graba Archivo GIFAnimado y video
    ani.save(narchivo+'_animado.gif', writer='pillow')
    #ani.save(narchivo+'.mp4')
    plt.show()
    

    GIF animado: [ Bisección ] [ Posicion Falsa ] [ Newton-Raphson ] [ Punto Fijo ] [ Secante ]

    ..


    Método del Punto Fijo con gráfico animado en Python

    Punto Fijo GIF animadoInstrucciones en Python

    # Algoritmo de punto fijo, Tabla
    # Los valores de [a,b] son seleccionados
    # desde la gráfica de la función
    # error = tolera
    
    import numpy as np
    
    def puntofijo_tabla(gx,a,tolera, iteramax=20,
                        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
        b = gx(a)
        tramo = abs(b-a)
        tabla = [[a,b,tramo]]
        if vertabla==True:
            print('método del Punto Fijo')
            print('i', ['xi','gi','tramo'])
            np.set_printoptions(precision)
            print(itera,np.array([a,b,tramo]))
        while(tramo>=tolera and itera<=iteramax):
            a = b
            b = gx(a)
            tramo = abs(b-a)
            itera = itera + 1
            if vertabla==True:
                print(itera,np.array([a,b,tramo]))
            tabla.append([a,b,tramo])
        respuesta = b
        # Valida respuesta
        if itera>=iteramax:
            respuesta = np.nan
            print('itera: ',itera,
                  'No converge,se alcanzó el máximo de iteraciones')
        tabla = np.array(tabla,dtype=float)
        return(respuesta,tabla)
    
    # PROGRAMA ----------------------
    
    # INGRESO
    fx = lambda x: np.exp(-x) - x
    gx = lambda x: np.exp(-x)
    
    a = 0
    b = 1
    tolera = 0.001
    
    # PROCEDIMIENTO
    [raiz,tabla] = puntofijo_tabla(gx,a,tolera, vertabla=True)
    
    # SALIDA
    print('raiz en :', raiz)
    
    # GRAFICA
    import matplotlib.pyplot as plt
    muestras = 21
    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(raiz, color='magenta',
                linestyle='dotted')
    plt.axhline(0)
    plt.xlabel('x')
    plt.ylabel('f(x)')
    plt.title('Punto Fijo')
    plt.grid()
    plt.legend()
    #plt.show()
    
    # GRAFICA CON ANIMACION ------------
    #import matplotlib.pyplot as plt
    import matplotlib.animation as animation
    
    xc = tabla[:,0]
    yc = tabla[:,1]
    
    # Inicializa parametros de trama/foto
    narchivo = 'PuntoFijo' # nombre archivo
    retardo = 700 # milisegundos entre tramas
    tramas = len(xc)
    
    # GRAFICA animada en fig_ani
    fig_ani, graf_ani = plt.subplots()
    dx = np.abs(b-a)
    dy = np.abs(gx(b)-gx(a))
    graf_ani.set_xlim([a-0.05*dx,b+0.05*dx])
    graf_ani.set_ylim([np.min(fi)-0.05*dy,np.max(fi)+0.05*dy])
    
    # Lineas y puntos base
    puntoa,  = graf_ani.plot(xc[0], yc[0], 'ro')
    puntob,  = graf_ani.plot(xc[0], xc[0], 'go') # y=x
    lineafx, = graf_ani.plot(xi,fi, color='blue',
                         linestyle='dashed', label='f(x)')
    lineagx, = graf_ani.plot(xi,gi, color='orange', label='g(x)')
    lineayx, = graf_ani.plot(xi,yi, color='green', label='y=x')
    linearaiz = graf_ani.axvline(raiz, color='magenta',linestyle='dotted')
    
    # Configura gráfica
    linea0 = graf_ani.axhline(0, color='k')
    graf_ani.set_title('Punto Fijo')
    graf_ani.set_xlabel('x')
    graf_ani.set_ylabel('f(x)')
    graf_ani.legend()
    graf_ani.grid()
    
    # Cada nueva trama
    def unatrama(i,xc,yc):
        # actualiza cada punto, [] porque es solo un valor
        puntoa.set_xdata([xc[i]])
        puntoa.set_ydata([yc[i]])
        puntob.set_xdata([xc[i]]) # y=x
        puntob.set_ydata([xc[i]])
        
        # actualiza cada flecha
        dx = xc[i+1]-xc[i]
        dy = yc[i]-xc[i]
        flecha01 = graf_ani.arrow(xc[i],xc[i], 0,dy,
                  length_includes_head = True,
                  head_width = 0.05*abs(dy),
                  head_length = 0.1*abs(dy))
        flecha02 = graf_ani.arrow(xc[i],yc[i], dx,0,
                  length_includes_head = True,
                  head_width = 0.05*abs(dx),
                  head_length = 0.1*abs(dx))
        return (puntoa, puntob, flecha01, flecha02)
    
    # Limpia trama anterior
    def limpiatrama():
        puntoa.set_ydata(np.ma.array([xc[0]], mask=True))
        puntob.set_ydata(np.ma.array([xc[0]], mask=True))
        return (puntoa, puntob)
    
    # contador de tramas
    i = np.arange(0, tramas-1,1)
    ani = animation.FuncAnimation(fig_ani,unatrama,
                                  i ,
                                  fargs=(xc, yc),
                                  init_func=limpiatrama,
                                  interval=retardo,
                                  blit=True)
    # Graba Archivo GIFAnimado y video
    ani.save(narchivo+'_animado.gif', writer='pillow')
    #ani.save(narchivo+'.mp4')
    plt.show()
    

    GIF animado: [ Bisección ] [ Posicion Falsa ] [ Newton-Raphson ] [ Punto Fijo ] [ Secante ]

    ..


    Método de la Secante con gráfico animado en Python

    Secante Metodo animado

    Instrucciones en Python

    # Método de la secante
    # Ejemplo 1 (Burden ejemplo 1 p.51/pdf.61)
    import numpy as np
    
    def secante_raiz_tabla(fx,a,b,tolera, iteramax=20,
                     vertabla=True, precision=6):
        '''fx en forma numérica lambda
        Los valores de [a,b] son seleccionados
        desde la gráfica de la función
        '''
        xi_1 = a
        xi = b
        itera = 0
        tramo = np.abs(xi-xi_1)
        tabla = []
        if vertabla==True:
            print('método de la Secante')
            print('i','[ x[i-1], xi, x[i+1], f[i-1], fi ]','tramo')
            np.set_printoptions(precision)
        while not(tramo<tolera or itera>iteramax):
            fi_1 = fx(xi_1)
            fi = fx(xi)
            xi1 = xi-fi*(xi_1 - xi)/(fi_1-fi)
            tramo = np.abs(xi1-xi)
            if vertabla==True:
                print(itera,np.array([xi_1,xi,xi1,fi_1,fi]),tramo)
            tabla.append([xi_1,xi,xi1,fi_1,fi])
            xi_1 = xi
            xi = xi1
            itera = itera + 1
        if itera>=iteramax:
            xi = np.nan
            print('itera: ',itera,
                  'No converge,se alcanzó el máximo de iteraciones')
        respuesta = xi
        tabla = np.array(tabla,dtype=float)
        return(respuesta,tabla)
    
    # PROGRAMA ----------------------
    
    # INGRESO
    fx = lambda x: x**3 + 4*x**2 - 10
    
    a = 1
    b = 4
    c = (a+b)/2
    tolera = 0.001
    tramos = 100
    
    # PROCEDIMIENTO
    [raiz,tabla] = secante_raiz_tabla(fx,c,b,tolera,
                                      vertabla=True)
    # SALIDA
    print('raiz en :', raiz)
    
    # GRAFICA
    import matplotlib.pyplot as plt
    muestras = 21
    xi = np.linspace(a,b,muestras)
    fi = fx(xi)
    xc = tabla[:,2]
    yc = fx(xc)
    
    plt.plot(xi,fi, label='f(x)')
    plt.plot([a,b],[fx(a),fx(b)],'o',
             color='red',label='[x0,f(x0)]')
    plt.scatter(xc,yc,color='orange', label='[c,f(c)]')
    plt.axhline(0)
    plt.xlabel('x')
    plt.ylabel('f(x)')
    plt.grid()
    plt.legend()
    #plt.show()
    
    
    # GRÁFICO CON ANIMACION ######
    import matplotlib.animation as animation
    
    xa = tabla[:,0]
    ya = tabla[:,3]
    xb = tabla[:,1]
    yb = tabla[:,4]
    xc = tabla[:,2]
    
    # Inicializa parametros de trama/foto
    narchivo = 'SecanteMetodo' # nombre archivo
    retardo = 700 # milisegundos entre tramas
    tramas = len(xa)
    
    # GRAFICA animada en fig_ani
    fig_ani, graf_ani = plt.subplots()
    graf_ani.set_xlim([a,b])
    graf_ani.set_ylim([np.min(fi),np.max(fi)])
    
    # Lineas y puntos base
    lineafx, = graf_ani.plot(xi,fi,label ='f(x)')
    
    puntoa, = graf_ani.plot(xa[0], ya[0],'o',
                            color='red', label='x[i-1]')
    puntob, = graf_ani.plot(xb[0], yb[0],'o',
                            color='green', label='x[i]')
    puntoc, = graf_ani.plot(xc[0], yc[0],'o',
                            color='orange', label='x[i+1]')
    
    lineatan1, = graf_ani.plot([xa[0],xb[0]],
                               [ya[0],yb[0]],
                               color='orange',label='secante ac')
    lineatan2, = graf_ani.plot([xc[0],xb[0]],
                               [0,yb[0]],
                               color='orange',label='secante cb')
    linea_a, = graf_ani.plot([xa[0],xa[0]],
                             [ya[0],0], color='magenta',
                             linestyle='dashed')
    linea_b, = graf_ani.plot([xb[0],xb[0]],
                             [0,yb[0]], color='magenta',
                             linestyle='dashed')
    # Configura gráfica
    linea0 = graf_ani.axhline(0, color='k')
    graf_ani.set_title('Método de la Secante')
    graf_ani.set_xlabel('x')
    graf_ani.set_ylabel('f(x)')
    graf_ani.legend()
    graf_ani.grid()
    
    # Cada nueva trama
    def unatrama(i,xa,ya,xb,yb,xc):
        # actualiza cada punto, [] porque es solo un valor
        puntoa.set_xdata([xa[i]]) 
        puntoa.set_ydata([ya[i]])
        puntob.set_xdata([xb[i]])
        puntob.set_ydata([yb[i]])
        puntoc.set_xdata([xc[i]])
        puntoc.set_ydata([0])
        # actualiza cada linea
        linea_a.set_ydata([ya[i], 0])
        linea_a.set_xdata([xa[i], xa[i]])
        linea_b.set_ydata([0, yb[i]])
        linea_b.set_xdata([xb[i], xb[i]])
        lineatan1.set_ydata([ya[i], 0])
        lineatan1.set_xdata([xa[i], xc[i]])
        lineatan2.set_ydata([0, yb[i]])
        lineatan2.set_xdata([xc[i], xb[i]])
    
        return (puntoa, puntob, puntoc,
                linea_a, linea_b,
                lineatan1, lineatan2,)
    
    # Limpia trama anterior
    def limpiatrama():
        puntoa.set_ydata(np.ma.array(xa, mask=True))
        puntob.set_ydata(np.ma.array(xb, mask=True))
        puntoc.set_ydata(np.ma.array(xc, mask=True))
        linea_a.set_ydata(np.ma.array([0,0], mask=True))
        linea_b.set_ydata(np.ma.array([0,0], mask=True))
        lineatan1.set_ydata(np.ma.array([0,0], mask=True))
        lineatan2.set_ydata(np.ma.array([0,0], mask=True))
        return (puntoa, puntob, puntoc,
                linea_a, linea_b,
                lineatan1, lineatan2,)
    
    # contador de tramas
    i = np.arange(0,tramas,1)
    ani = animation.FuncAnimation(fig_ani,unatrama,
                                  i ,
                                  fargs=(xa, ya,
                                         xb, yb,
                                         xc),
                                  init_func=limpiatrama,
                                  interval=retardo,
                                  blit=True)
    # Graba Archivo GIFAnimado y video
    ani.save(narchivo+'_animado.gif', writer='pillow')
    #ani.save(narchivo+'.mp4')
    plt.show()
    

    GIF animado: [ Bisección ] [ Posicion Falsa ] [ Newton-Raphson ] [ Punto Fijo ] [ Secante ]

  • 2.6 Sistemas de Ecuaciones no lineales - Newton-Raphson

    Referencia: Chapra 6.5 p162, Chapra Ejercicio 6.11 p166

    Con el método de Newton-Raphson para múltiples ecuaciones, determine las raíces para:

    x^2+xy =10 y + 3xy^2 = 57

    Observe que un par correcto de raíces es x=2 y y=3.
    Use como valores iniciales x=1.5, y=3.5

    Planteamiento

    Las ecuaciones se expresan de la forma f(x,y) = 0

    x^2+xy -10 = 0 y + 3xy^2 -57 = 0

    Se puede usar extensiones de los métodos abiertos para resolver ecuaciones simples, por ejemplo Newton-Raphson.

    u_{i+1} = u_i + (x_{i+1}-x_i)\frac{\partial u_i}{\partial x} + (y_{i+1}-y_i) \frac{\partial u_i}{\partial y} v_{i+1} = v_i + (x_{i+1}-x_i)\frac{\partial v_i}{\partial x} + (y_{i+1}-y_i) \frac{\partial v_i}{\partial y}

    ecuaciones que se pueden reordenar y encontrar la solución a partir de la matriz Jacobiano.

    Instrucciones en Python

    Usando un algoritmo para resolver el Jacobiano y estimar los puntos luego de cada iteración se obtienen:

    iteración:  1
    Jacobiano con puntos iniciales: 
    [ 6.5   1.5 ]
    [           ]
    [36.75  32.5]
    determinante:  156.12499999999994
    puntos xi,yi: 2.03602882305845 2.84387510008006
    error: 0.656124899919936
    iteración:  2
    Jacobiano con puntos iniciales: 
    [6.91593274619696  2.03602882305845]
    [                                  ]
    [24.2628767545662  35.7412700376474]
    determinante:  197.78430344142245
    puntos xi,yi: 1.99870060905582 3.00228856292451
    error: 0.158413462844444
    iteración:  3
    Jacobiano con puntos iniciales: 
    [6.99968978103616  1.99870060905582]
    [                                  ]
    [27.0412098452019  37.0040558756713]
    determinante:  204.96962918261596
    puntos xi,yi: 1.99999998387626 2.99999941338891
    error: 0.00228914953559523
    iteración:  4
    Jacobiano con puntos iniciales: 
    [6.99999938114143  1.99999998387626]
    [                                  ]
    [26.9999894410015  36.9999926704397]
    determinante:  204.9999473486533
    puntos xi,yi: 1.99999999999998 3.00000000000008
    error: 5.86611161867978e-7
    Resultado: 
    1.99999999999998 3.00000000000008
    

    Algoritmo presentado para dos ecuaciones y dos incógnitas, en la unidad 3 se puede ampliar la propuesta. Revisar el método de Gauss-Seidel y Jacobi.

    # Ejercicio Chapra Ej:6.11
    # Sistemas de ecuaciones no lineales
    # con método de Newton Raphson para xy
    
    import numpy as np
    import sympy as sym
    
    def matrizJacobiano(variables, funciones):
        n = len(funciones)
        m = len(variables)
        # matriz Jacobiano inicia con ceros
        Jcb = sym.zeros(n,m)
        for i in range(0,n,1):
            unafi = sym.sympify(funciones[i])
            for j in range(0,m,1):
                unavariable = variables[j]
                Jcb[i,j] = sym.diff(unafi, unavariable)
        return Jcb
    
    # PROGRAMA ----------
    # INGRESO
    x = sym.Symbol('x')
    y = sym.Symbol('y')
    
    f1 = x**2 + x*y - 10
    f2 = y + 3*x*(y**2)-57
    
    x0 = 1.5
    y0 = 3.5
    
    tolera = 0.0001
    
    # PROCEDIMIENTO
    funciones = [f1,f2]
    variables = [x,y]
    n = len(funciones)
    m = len(variables)
    
    Jxy = matrizJacobiano(variables, funciones)
    
    # valores iniciales
    xi = x0
    yi = y0
    
    # tramo inicial, mayor que tolerancia
    itera = 0
    tramo = tolera*2
    
    while (tramo>tolera):
        J = Jxy.subs([(x,xi),(y,yi)])
    
        # determinante de J
        Jn = np.array(J,dtype=float)
        determinante =  np.linalg.det(Jn)
    
        # iteraciones
        f1i = f1.subs([(x,xi),(y,yi)])
        f2i = f2.subs([(x,xi),(y,yi)])
    
        numerador1 = f1i*Jn[n-1,m-1]-f2i*Jn[0,m-1]
        xi1 = xi - numerador1/determinante
        numerador2 = f2i*Jn[0,0]-f1i*Jn[n-1,0]
        yi1 = yi -numerador2/determinante
        
        tramo = np.max(np.abs([xi1-xi,yi1-yi]))
        xi = xi1
        yi = yi1
    
        itera = itera +1
        print('iteración: ',itera)
        print('Jacobiano con puntos iniciales: ')
        sym.pprint(J)
        print('determinante: ', determinante)
        print('puntos xi,yi:',xi,yi)
        print('error:',tramo)
        
    # SALIDA
    print('Resultado: ')
    print(xi,yi)
    
  • 2.5.1 Método de la Secante - Ejemplo con Python

    [ Secante ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ]
    ..


    1. Ejercicio

    Referencia: Burden 2.1 ejemplo 1 p38

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

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

    Secante Metodo animadoLa gráfica se realiza iniciando con los valores 2.5 y 4 para observar mejor el efecto gráfico en las iteraciones.

    [ Secante ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [función]
    ..


    2. Desarrollo Analítico

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

    Se probará con valores iniciales dentro del tramo [1,2], tolera = 0.001

    itera = 0

    X-1 = 1,  X0 =2,

    f(1) = (1)^3 + 4(1)^2 -10 = -5 f(2) = (2)^3 + 4(2)^2 -10 =14 x_{1}= 2 - 14\frac{(1 - 2)}{(-5 - 14)} = 1.2632 errado = |1.2632 - 2| = 0.7368

    itera = 1

    X0 = 2,  X1 = 1.2632

    f(2) =14 f(1.2632) = ( 1.2632)^3 + 4( 1.2632)^2 -10 =-1.6023 x_{2}= 1.2632 - (-1.6023)\frac{(2 - 1.2632)}{(14 -(-1.6023))} =1.3388 errado = |1.3388 - 1.2632 | = 0.0756

    itera = 2

    X1 =1.2632, X2 =1.3388

    f(1.2632) = -1.6023 f(1.3388) = (1.3388)^3 + 4(1.3388)^2 -10 =-0.4304 x_{3}= 1.3388 - (-0.4304)\frac{(1.2632 - 1.3388)}{(-1.6023 -(-0.4304))} = 1.3666 errado = | 1.3388- 1.3666| =0.0277

    desarrollando una tabla con el algoritmo se tiene:

    método de la Secante
    i [ x[i-1], xi, x[i+1], f[i-1], fi ] tramo
    0 [ 1.      2.      1.2632 -5.     14.    ] 0.736842105263158
    1 [ 2.      1.2632  1.3388 14.     -1.6023] 0.0756699440909967
    2 [ 1.2632  1.3388  1.3666 -1.6023 -0.4304] 0.027788555891506528
    3 [ 1.3388  1.3666  1.3652 -0.4304  0.0229] 0.001404492087488718
    4 [ 1.3666e+00  1.3652e+00  1.3652e+00  2.2909e-02 -2.9907e-04] 1.809847900258177e-05
    raíz:  1.3652300011108591
    tramo: 1.809847900258177e-05
    itera: 5

    [ Secante ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ]
    ..


    3. Algoritmo en Python

    Para identificar los puntos consecutivos con variables simples, se usa la siguiente conversión:

    xi-1 xi xi+1
    xa xb xc

    Instrucciones en Python

    # Método de secante
    import numpy as np
    
    # INGRESO
    fx = lambda x: x**3 + 4*x**2 - 10
    a = 1 # intervalo de búsqueda
    b = 2
    tolera = 0.001
    iteramax = 20
    
    # PROCEDIMIENTO
    xa = a
    xb = b
    itera = 0
    tramo = 2*tolera
    while tramo>tolera and itera<iteramax:
        fa = fx(xa)
        fb = fx(xb)
        xc = xb - fb*(xa - xb)/(fa - fb)
        tramo = abs(xc - xb)
        xa = xb
        xb = xc
        itera = itera + 1
    
    if itera>=iteramax: # revisa convergencia
        xc = np.nan
    
    # SALIDA
    print('raiz:',xc)
    print('tramo:',tramo)
    print('itera:',itera)
    

    [ Secante ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ]
    ..


    4. Función en Python

    Una función resumida, se recomienda realizar la validación de cambio de signo entre [a,b]

    # Método de secante
    import numpy as np
    
    def secante_raiz(fx,a,b,tolera, iteramax=20,
                     vertabla=True, precision=6):
        '''fx en forma numérica lambda
        Los valores de [a,b] son seleccionados
        desde la gráfica de la función
        '''
        xa = a
        xb = b
        itera = 0
        tramo = np.abs(xb-xa)
        if vertabla==True:
            print('método de la Secante')
            print('i','[ x[i-1], xi, x[i+1], f[i-1], fi ]','tramo')
            np.set_printoptions(precision)
        while not(tramo<tolera or itera>iteramax):
            fa = fx(xa)
            fb = fx(xb)
            xc = xb - fb*(xa - xb)/(fa - fb)
            tramo = abs(xc - xb)
            if vertabla==True:
                print(itera,np.array([xa,xb,xc,fa,fb]),tramo)
            xa = xb
            xb = xc
            itera = itera + 1
        if itera>=iteramax:
            xc = np.nan
            print('itera: ',itera,
                  'No converge,se alcanzó el máximo de iteraciones')
        respuesta = xc
        
        return(respuesta)
    
    # INGRESO
    fx = lambda x: x**3 + 4*x**2 - 10
    a = 1
    b = 2
    tolera = 0.001
    
    # PROCEDIMIENTO
    respuesta = secante_raiz(fx,a,b,tolera,
                             vertabla=True,precision=4)
    # SALIDA
    print('raíz: ', respuesta)
    

    5. Función en librería Scipy.optimize.newton - Secante

    El método de la secante se encuentra implementado en Scipy en la forma de algoritmo de newton, que al no incluir como parámetro de la función para la derivada de f(x), usa el método de la secante:

    >>> import scipy as sp
    >>> sp.optimize.newton(fx,xa, tol=tolera)
    1.3652320383201266

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


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

  • 2.5 Método de la Secante - Concepto

    [ Secante ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ]
    ..


    Método de la Secante

    Referencia: Chapra 6.3 p154, Burden 2.3 p54,

    El método de la secante busca evitar un posible inconveniente en el desarrollo el algoritmo para método de Newton-Raphson, que es el implementar la evaluación de la expresión de la primera derivada.

    Secante Metodo animado

    La derivada de la función f(x) también se puede aproximar mediante una diferencia finita dividida hacia atrás:

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

    método de la secante animado

    la que se sustituye en la ecuación del método de Newton-Raphson para obtener:

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

    Los métodos de la secante y de la falsa posición tienen ecuaciones idénticas, usan dos valores iniciales x[i-1] , x[i], para proyectar el nuevo valor de x[i+1].

    Sin embargo, existe una diferencia importante entre ambos métodos en la forma en que uno de los valores iniciales se reemplaza por la nueva aproximación.

    En el método de la falsa posición, la última aproximación de la raíz reemplaza cualquiera de los valores iniciales que dé un valor de la función con el mismo signo. En consecuencia, las dos aproximaciones siempre encierran a la raíz, es un método cerrado y siempre converge.

    En el método de la secante se reemplaza los valores en secuencia estricta: xi reemplaza a x[i – 1], con el nuevo valor x[i+1] se reemplaza a xi . Por lo que, algunas veces los dos valores están en el mismo lado de la raíz y en ciertos casos esto puede llevar a divergencias.

    Observación: ¿Cuál es la diferencia con el método de Newton-Raphson?

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


    Recta Secante para la gráfica

    Si se quiere dibujar la recta secante, se inicia con la ecuación de la recta usando el valor de la pendiente del triángulo presentado en la gráfica del método, cambie la constante b por b0.

    y = mx + b

    Para identificar los puntos consecutivos con variables simples, se usa la siguiente conversión

    xi-1 xi xi+1
    xa xb xc
    m = \frac{f(x_a)-f(x_b)}{x_a-x_b}

    Se usa un punto conocido, por ejemplo (xb,fb) y se encuentra b0

    f(x_b) = m x_b + b_0 f(x_b) - m x_b = b_0

    reordenando, se usan tres expresiones para disponer de la función de la recta secante.

    m = \frac{f(x_a)-f(x_b)}{x_a-x_b} b_0 = f(x_b) - m x_b rsc(x)= m x + b_0

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

  • 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

    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

    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

    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 ]