Métodos numéricos, cálculo de raíces mediante método de la bisección

En este ejemplo se presenta el método de la bisección para calcular las raíces de una ecuación de una variable.

Enunciado del problema

Supongamos que tenemos una función continua f definida en el intervalo [a, b], con f(a) y f(b) de signos distintos. Entonces por el corolario V.1 del Teorema del Valor Intermedio, existe un punto c, a < c < b, tal que f(c) = 0. Aunque el procedimiento sirve para el caso en el que f(a) y f(b) tienen signos opuestos y hay más de una raíz en el intervalo [a, b], por simplicidad se supondrá que la raíz en este intervalo es única.

Corolario V.1

Si f ∈ C[a, b] asume valores de signo opuesto en los extremos de un intervalo [α, β], es decir, f(α) · f(β) < 0, entonces el intervalo contendrá al menos una raíz de la ecuación f(x) = 0; en otras palabras, habrá al menos un número c ∈ (α, β) tal que f(c) = 0.

Resolución:

Para la resolución de estos problemas vamos a requerir la librería cexprtk de Python que nos va a permitir interpretar las funciones mátematicas a partir de Strings, la importamos de la siguiente manera:

import cexprtk as cx

Ahora vamos a definir una función que haga uso de la librería que acabamos de importar para evaluar las expresiones matemáticas en un unto dado.

def comp(exp, point, var="x"):
   return cx.evaluate_expression(exp, {var:point})

Definiremos tres funciones para calcular las raíces mediante el método de la bisección, que se van a diferenciar por los parámetros que reciben y validaciones que realizan así:

1. bisectionFull: Esta función recibe cuatro parámetros, la función que se va a evaluar, el intervalo en el que se va a evaluar la función, el número de veces que se va a ejecutar el método y la tolerancia para determinar el cálculo de la raíz. Las valicaciones se hacen o bien sobre la cantidad máxima de veces qe se va a ejecutar el método o la tolerancia, el que se cumpla primero.
2. bisectionSteps: Esta función recibe también la función a evaluar, el intervalo de evaluación y la cantidad de veces que se va a repetir el proceso. No se recibe ni valida la tolerancia. Aquí solamente se valida la cantidad de veces a repetir el método.
3. bisectionTol: Esta función recibe la función a evaluar, el intervalo de evaluación y la tolerancia. Aquí solamente se valida la tolerancia para encontrar la raíz. Se usa esta función por defecto con la tolerancia 0.0001 por defecto en caso de que no se defina ni la cantidad de pasos o la tolerancia.

Función bisectionFull

def bisectionFull(exp, interval, steps, tol):
   a = interval[0]
   b = interval[1]
   for i in range(steps):
     mean = (b-a)/2
     c = a + mean
     resC = comp(exp, c)
     if resC == 0:
       print(f"Se encontró la raíz [{c}] en la iteración [{i+1}]")
       break
     if mean < tol:
       print(f"Se encontró la raíz [{c}] en con el valor [{mean}] menor que la tolerancia [{tol}] en la iteración [{i+1}]")
       break
     if comp(exp, a) > comp(exp, c) > 0:
       a = c
     else:
       b = c
   else:
     print(f"No se pudo encontrar el valor de la raíz luego de las [{steps}] iteraciones")

Función bisectionSteps

def bisectionSteps(exp, interval, steps):
   a = interval[0]
   b = interval[1]
   for i in range(steps):
     mean = (b-a)/2
     c = a + mean
     resC = comp(exp, c)
     if resC == 0:
       print(f"Se encontró la raíz [{c}] en la iteración [{i+1}]")
       break
     if comp(exp, a) > comp(exp, c) > 0:
       a = c
     else:
       b = c
   else:
     print(f"No se pudo encontrar el valor de la raíz luego de las [{steps}] iteraciones")

Función bisectionTol

def bisectionTol(exp, interval, tol=0.0001):
   a = interval[0]
   b = interval[1]
   mean = tol+1
   c = 0
   i = 0
   while mean > tol:
     i += 1
     mean = (b-a)/2
     c = a + mean
     resC = comp(exp, c)
     if resC == 0:
       print(f"Se encontró la raíz [{c}] en la iteración [{i+1}]")
       break
     if mean < tol:
       print(f"Se encontró la raíz [{c}] en con el valor [{mean}] menor que la tolerancia [{tol}] en la iteración [{i}]")
       break
     if comp(exp, a) > comp(exp, c) > 0:
       a = c
     else:
       b = c
   else:
     print(f"No se pudo encontrar el valor de la raíz para la tolerancia [{tol}]")

Ahora vamos a definir la función bisection que serña la puerta de entrada ara el cálculo de las raices en las funciones previamente definidas. En esta función se van a realizar validaciones de datos, para que no sean vacíos o nulos. También se verificará que el producto de los resultados de la función evaluada en los puntos del intervalo sea menor a cero.

Se verificará también que la expresión ingresada sea válida, en caso de no cumplir estos requisitos se lanzará una excepción en la que se indica el requerimiento que no se h acumplido. Si los requisitos se cumplen, se validarán los parámetros recibidos para enviar a resolver el problema en alguna de las funciones previamente definidas. La función se define de la siguiente manera:

def bisection(exp, interval, steps=0, tol=0):
   if not exp:
     raise Exception("El valor de exp debe ser diferente de nulo o vacío")
   elif not interval or len(interval) < 2:
     raise Exception("El valor de interval debe ser diferente de nulo o vacío")
   elif steps < 0 or steps == None:
     raise Exception("El valor de interval debe ser mayor que 0 y diferente de nulo")
   elif tol < 0 or tol == None:
     raise Exception("El valor de tol debe ser mayor que 0 y diferente de nulo")
   else:
     print(f"Se va a evaluar la expresión [{exp}]")
   try:
     if comp(exp, interval[0]) * comp(exp, interval[1]) > 0:
       raise Exception(f"No se puede ejecutar el método con el intervalo {interval}")
   except:
     raise Exception(f"No se pudo evaluar la explresión [{exp}]")
   if steps > 0 and tol > 0:
     bisectionFull(exp, interval, steps, tol)
   elif steps > 0 and tol == 0:
     bisectionSteps(exp, interval, steps)
   elif tol > 0 and steps == 0:
     bisectionTol(exp, interval, tol)
   else:
     bisectionTol(exp, interval, tol)

Con esto tendremos definida la funcionalidad de nuestro algoritmo completo, ahora vamos a definir datos ejemplo para probar el funcionamiento del algoritmo:

exp = "6x^2 - 8x"
interval = [0, 5]
steps = 2000
tol = 0.0001
bisection(exp, interval, steps)

La salida de la ejeución de nuestro algoritmo se muestra en la siguiente imagen: