Pybonacci

Computación Científica con Python en castellano

Ecuaciones no lineales: método de bisección y método de Newton en Python

with 7 comments

En este artículo vamos a ver cómo implementar en Python el método de bisección y el método de Newton, dos métodos iterativos clásicos para hallar raíces de ecuaciones no lineales de la forma f(x) = 0, con f: [a, b] \longrightarrow \mathbb{R} y f \in C^1([a, b]). Estos métodos y muchos otros más refinados están ya implementados en multitud de bibliotecas muy utilizadas, sin ir más lejos en el módulo optimize del paquete Scipy (referencia).

Crearemos un módulo ceros.py en el que incluiremos los dos métodos que vamos a desarrollar aquí, y así veremos un ejemplo de código limpio y fácilmente reutilizable.

Módulo ceros.py

Vamos a ver la anatomía de un módulo en Python. Este es el código del archivo:

# -*- coding: utf-8 -*-

"""Búsqueda de raíces

Este módulo contiene métodos para la búsqueda de raíces de ecuaciones de la
forma f(x) = 0, con f función real de variable real, continua y de derivada
continua.

"""

def biseccion():
    """Método de bisección"""
    pass

def newton():
    """Método de Newton"""
    pass

Analicemos el código:

  • La primera línea es fundamental en Python 2 si vamos a usar caracteres que se salen de ASCII. En ella especificamos que la codificación sea UTF-8, como viene recogido en la PEP 263.
  • La línea entre triples comillas dobles es lo que se llama docstring o cadena de documentación del módulo, y siempre es una cadena que se pone al principio del archivo, como viene especificado en la PEP 257. Lo de usar triples comillas nos aseguramos de que podemos incluir saltos de línea.
  • Definimos dos funciones, también con su docstring. La palabra clave pass se utiliza para que cuadre el sangrado del código: no hace nada.
  • Las convenciones en cuanto a espacios, longitud de las líneas, etc. están definidas en la PEP 8, que viene a ser como el manual de estilo oficial para programas en Python. Puedes comprobar si tu código se adhiere a esta convención con la herramienta pep8 disponible en el Índice de Paquetes de Python (PyPI).

Para utilizar este módulo, simplemente lenzaríamos un intérprete Python en la carpeta donde esté el archivo ceros.py y escribiríamos:

>>> import ceros
>>> dir(ceros)
['__builtins__', '__doc__', '__file__', '__name__', '__package__', 'biseccion', 'newton']
>>> print ceros.__doc__
Búsqueda de raíces

Este módulo contiene métodos para la búsqueda de raíces de ecuaciones de la
forma f(x) = 0, con f función real de variable real, continua y de derivada
continua.
>>> print ceros.biseccion
<function biseccion at 0x7fc17efb5668>

Ahora no hay más que implementar estos métodos.

Método de la bisección

Descripción y algoritmo

El método de bisección es un método geométricamente muy intuitivo: partiendo de un intervalo [a, b] tal que f(a) f(b) < 0 (es decir: la función cambia de signo en el intervalo), se va dividiendo en dos generando una sucesión de intervalos encajados hasta que se converge con la precisión deseada a la raíz de la ecuación que, como asegura el teorema de Bolzano, tiene que existir; es decir, el \alpha \in ]a, b[ tal que f(\alpha) = 0.

El algoritmo del método de bisección sería el siguiente:

  1. Sean fa y b.
    1. Sea c \leftarrow \frac{a + b}{2}.
    2. Si f(c) = 0, terminar.
    3. Si f(c) f(a) < 0, entonces b \leftarrow c y volver al paso 1.
    4. Si no, entonces a \leftarrow c y volver al paso 1.

Este método, aunque es lento, tiene convergencia garantizada.

Método de bisección en Python

A la vista del algoritmo anterior, ya podemos implementar el método de bisección en Python. Nos vamos a saltar la parte de escribir el pseudocódigo.

Nota: Cambiado por sugerencia de David para evitar errores de precisión.

import numpy as np

def biseccion(f, a, b, tol=1.0e-6):
    """Método de bisección

    Halla una raíz de la función f en el intervalo [a, b] mediante el método
    de bisección.

    Argumentos
    ----------
    f - Función, debe ser tal que f(a) f(b) < 0
    a - Extremo inferior del intervalo
    b - Extremo superior del intervalo
    tol (opcional) - Cota para el error absoluto de la x

    Devuelve
    --------
    x - Raíz de f en [a, b]

    Excepciones
    -----------
    ValueError - Intervalo mal definido, la función no cambia de signo en el
                 intervalo o cota no positiva

    Ejemplos
    --------
    >>> def f(x): return x ** 2 - 1
    ...
    >>> biseccion(f, 0, 2)
    1.0
    >>> biseccion(f, 0, 5)
    1.000000238418579
    >>> biseccion(f, -2, 2)
    Traceback (most recent call last):
        File "<stdin>", line 1, in <module>
        File "ceros.py", line 35, in biseccion
        raise ValueError("La función debe cambiar de signo en el intervalo")
    ValueError: La función debe cambiar de signo en el intervalo
    >>> biseccion(f, -3, 0, tol=1.0e-12)
    -1.0000000000004547

    """
    if a > b:
        raise ValueError("Intervalo mal definido")
    if f(a) * f(b) >= 0.0:
        raise ValueError("La función debe cambiar de signo en el intervalo")
    if tol <= 0:
        raise ValueError("La cota de error debe ser un número positivo")
    x = (a + b) / 2.0
    while True:
        if b - a < tol:
            return x
        # Utilizamos la función signo para evitar errores de precisión
        elif np.sign(f(a)) * np.sign(f(x)) > 0:
            a = x
        else:
            b = x
        x = (a + b) / 2.0

Vamos a señalar algunas cosas:

  • Como se puede observar, el código se parece bastante al algoritmo.
  • Podemos pasar una función como argumento sin ningún esfuerzo.
  • Hemos incluido en la documentación información sobre los argumentos que recibe la función, de manera que si escribimos en el intérpretehelp(ceros.biseccion) podremos consultarla.
  • No he incluido un número máximo de iteraciones: fijada una precisión, el método convergerá siempre.

Por último, nótese que el programa produce un error si el límite inferior de intervalo es mayor que el límite superior, si la función no cambia de signo o si la cota de error es un número negativo. Python hace muy sencillas este tipo de construcciones: se denomina manejo de excepciones.

Podíamos haber tomado la decisión de no incluir este manejo de errores, y dejar que el programa falle inesperadamente si el usuario hace «cosas raras». Esto que lo decida cada uno.

Para gestionar los errores que se pueden producir, utilizamos el bloque try...except:

try:
    biseccion(f, a, b)
except ValueError:
    pass  # Este bloque se ejecuta si se produce un error

Gracias a esta característica de Python podemos evitar cosas como que una división por cero mate al programa, por ejemplo.

Método de Newton

Descripción y algoritmo

El método de Newton o método de Newton-Raphson linealiza la función a cada paso utilizando su derivada, que se debe proporcionar como argumento, para hallar la raíz de la ecuación en las proximidades de un punto inicial x_0. Este método puede no converger, pero si el punto inicial está lo suficientemente próximo a la raíz, la convergencia será muy rápida.

El algoritmo del método de Newton es este:

  1. Sean ff' y x0.
  2. Sea x \leftarrow x_0.
    1. Sea x \leftarrow x - \frac{f(x)}{f'(x)}.
    2. Si f(x) = 0, terminar.
    3. Si no, volver al paso 2.

Método de Newton en Python

Aquí, debido a que el método no tiene garantizada la convergencia, habrá que considerar un número máximo de iteraciones y cotas de error tanto para la raíz como para el valor de la función.

El código sería este:

def newton(f, df, x_0, maxiter=50, xtol=1.0e-6, ftol=1.0e-6):
    """Método de Newton

    Halla la raíz de la función f en el entorno de x_0 mediante el método de
    Newton.

    Argumentos
    ----------
    f - Función
    df - Función, debe ser la función derivada de f
    x_0 - Punto de partida del método
    maxiter (opcional) - Número máximo de iteraciones
    xtol (opcional) - Cota para el error relativo para la raíz
    ftol (opcional) - Cota para el valor de la función

    Devuelve
    --------
    x - Raíz de la ecuación en el entorno de x_0

    Excepciones
    -----------
    RuntimeError - No hubo convergencia superado el número máximo de
                   iteraciones
    ZeroDivisionError - La derivada se anuló en algún punto
    Exception - El valor de x se sale del dominio de definición de f

    Ejemplos
    --------
    >>> def f(x): return x ** 2 - 1
    ...
    >>> def df(x): return 2 * x
    ...
    >>> newton(f, df, 2)
    1.000000000000001
    >>> newton(f, df, 5)
    1.0
    >>> newton(f, df, 0)
    Traceback (most recent call last):
      File "<stdin>", line 1, in <module>
      File "ceros.py", line 102, in newton
        dx = -f(x) / df(x)  # ¡Aquí se puede producir una división por cero!
    ZeroDivisionError: float division by zero

    """
    x = float(x_0)  # Se convierte a número de coma flotante
    for i in xrange(maxiter):
        dx = -f(x) / df(x)  # ¡Aquí se puede producir una división por cero!
                            # También x puede haber quedado fuera del dominio
        x = x + dx
        if abs(dx / x) < xtol and abs(f(x)) < ftol:
            return x
    raise RuntimeError("No hubo convergencia después de {} \
                        iteraciones".format(maxiter))

Se deja como ejercicio (qué placentero es decir esto) programar el método de la secante que se utilice como alternativa si no se dispone de la derivada de la función.

¡Y con eso terminamos el artículo de hoy! Espero que te haya resultado provechoso, no olvides comentar, retwittear y demás zarandajas :)

Referencias

  • RIVAS, Damián y VÁZQUEZ, Carlos. Elementos de Cálculo Numérico. ADI, 2010.
About these ads

Written by Juanlu001

18 de abril de 2012 a 19:53

Publicado en Tutoriales

Etiquetado con , , ,

7 comentarios

Suscríbete a los comentarios mediante RSS.

  1. Hay dos problemas en el método de la bisección:

    elif f(a) * f(x) < 0:

    Siendo f(a) y f(x) ambos no nulos, esa expresión se puede evaluar como cierta. Y es que si multiplicamos dos números pequeños, el resultado se puede ir por debajo de la precisión y dar un cero. Es mejor comparar el signo. Algo así como:

    elif np.sign(f(a))!= np.sign(f(x)):

    o creando tu propia función signo.

    La otra es también aplicable al método de Newton. Lo que nos interesa es el valor de la x, no f(x). Se debe iterar hasta que a-b sea menor que una cierta tolerancia. Esto se puede implementar en el método de Newton de la misma manera que en la bisección, pero en lugar de elegir la nueva raíz prueba en el medio, se hace usando la derivada.

    Davidmh

    30 de abril de 2012 at 21:51

    • David, tienes toda la razón con lo que has dicho del método de bisección y lo he cambiado. En lo del método de Newton no te acabo de pillar la idea, pero creo que al comprobar el error relativo de la solución y el valor de la función puede ser suficiente (así lo hicimos en clase de Cálculo Numérico).

      Juanlu001

      30 de abril de 2012 at 22:53

  2. Por cierto, el mejor libro de métodos numéricos que conozco es el Numerical Recipes, de H. Press et al. Ahí tiene todo discutido desde un enfoque eminentemente práctico, y sólo los que funcionan, sin zaranjadas académicas.

    Davidmh

    30 de abril de 2012 at 21:52

  3. [...] con Python. Ya vimos hace tiempo cómo encontrar el mínimo de una función con SciPy, y también cómo implementar los métodos de la bisección y de Newton en Python. Ahora, además, exploraremos el caso de sistemas de [...]

  4. Hola hay un error en el método de la bisección, por ejemplo, al correr tu programa para

    def f(x):

        return x ** 2 - 1

    y biseccion(f, 0, 2) te devuelve 1.9999995231628418, cuando debería ser 0.9999995231628418.

    Esto se puede corregir cambiando


    elif np.sign(f(a)) == np.sign(f(x)):

        a=x

    else:

        b=x

    o con


    elif np.sign(f(a)) * np.sign(f(x)) > 0:

        a = x

    else:

        b = x

    Muchas gracias por los post, saludos.

    Miguel Ángel G

    22 de febrero de 2014 at 7:29

    • Hola, he editado tu comentario para que se vea correctamente la indentación.

      Juanlu (autor del post) anda por el otro lado del canal aprendiendo magia negra de diferentes expertos por lo que la semana que viene, cuando vuelva, lo corregirá.

      Muchas gracias por el aviso.

      Kiko

      22 de febrero de 2014 at 12:34

    • ¡Hola Miguel Ángel, y perdón por el retraso! Tienes toda la razón, está ya cambiado. ¡Saludos y gracias!

      Juanlu001

      19 de marzo de 2014 at 16:34


¡Deja un comentario!

Introduce tus datos o haz clic en un icono para iniciar sesión:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s

Seguir

Recibe cada nueva publicación en tu buzón de correo electrónico.

Únete a otros 1.103 seguidores

%d personas les gusta esto: