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:

Ejemplo cálculo de distancias, euclideana, de manhattan y de chebyshev

En este ejemplo se presenta la resolución de algunos ejercicios matemáticos para calcular la distancia entre dos puntos.

Distancia euclideana

La distancia euclideana entre dos puntos n dimensionales se define de la siguiente manera:

Para la resolución de estor problemas vamos a requerir la librería math de pyrhon la cual vamos a importar:

import math

Ahora para resolver el problema de la distancia euclideana definimos la siguiente función:

def euclidean(x,y):
   total = 0
   for i in range(len(x)):
     diff = x[i] - y[i]
     total = total + math.pow(diff, 2)
   return math.sqrt(total)

Distancia de manhattan

La distancia de manhattan entre dos puntos n dimensionales se define de la siguiente manera:

La función a continuación corresponde al algoritmo que calcula la distancia de manhattan

def manhattan(x,y):
   total = 0
   for i in range(len(x)):
     diff = x[i] - y[i]
     total = total + abs(diff)
   return total

Distancia de chebyshev

La distancia de chebyshev entre dos puntos n dimensionales se define de la siguiente manera:

La función a continuación corresponde al algoritmo que calcula la distancia de chebyshev:

def chebyshev(x, y):
   diffs = []
   for i in range(len(x)):
     diffs.append(abs(x[i] - y[i]))
   return max(diffs)

PONIENDOLO TODO JUNTO...

Ahora vamos a crear una función que nos permita validar que los dos vectores tienen la misma dimensión, es decir el mismo número de elementos. Esta recibirá los dos arreglos y validará la cantidad de elementos así:

def validateParams(x,y):
   if len(x) != len(y):
     return False
   return True

Para dinamizar la elección de el método a utilizar crearemos una función que reciba tres parámetros, el vector x, el vector y, y un string indicando la función a utilizar para calcular la distancia. Además, en caso de que no se especifique el nombre del método a utilizar se tomará por defecto el método euclideano.

Hay que tomar en cuenta que el nombre puede ser escrito incorrectamente, por eso haremos uso de la funcionalidad raise de Python para lanzar una excepción indicando que elk nombre ha sido mal escrito.

La función quedaría de la siguiente manera:

def distance(x, y, method='euclidean'):
   if not validateParams(x,y):
     return 0
   if method == 'euclidean':
     return euclidean(x, y)
   elif method == 'manhattan':
     return manhattan(x, y)
   elif method == 'chebyshev':
     return chebyshev(x, y)
   else:
     raise Exception("Debe seleccionar uno de los métodos [euclidean, manhattan, chebyshev]")

Ahora es tiempo de probar nuestros algoritmos, para lo cual definiremos dos vectores de la siguiente manera:

x = [9,3,5,3,-8]
y = [4,21,-5,7,11]

Ahora llamaremos a nuestra función de entrada distance con cada uno de los métodos e imprimiremos los resultados por consola:

print( f"Distancia euclideana: {distance(x,y,'euclidean')}")
print( f"Distancia de manhattan: {distance(x,y,'manhattan')}")
print( f"Distancia de chebyshev: {distance(x,y,'chebyshev')}")

El resultado de la ejecución se muestra en la siguiente imagen: