Consultar el Capítulo 19: Root Finding.
El objetivo del método de bisección es encontrar una raíz, o solución, para una ecuación de la forma \(f (x) = 0\), para una función \(f\) dada. Una raíz de esta ecuación también recibe el nombre de cero de la función \(f\).
La primera técnica, basada en el teorema de valor intermedio, recibe el nombre de bisección o método de búsqueda binaria.
Supongamos que \(f: [a, b]\subset \mathbb{R}\rightarrow \mathbb{R}\) es una función continua con \(f(a)\) y \(f(b)\) de signos opuestos. El teorema del valor intermedio implica que existe un número \(p\) en \((a,b)\) con \(f (p)=0\). Inicialmente suponemos que la raíz en este intervalo es única. El método realiza repetidamente una reducción a la mitad (o bisección) de los subintervalos de \([a, b]\) y, en cada paso, localizar la mitad que contiene \(p\).
Para comenzar, sean \(a_1=a\), \(b_1=b\) y \(p_1\) el punto medio de \([a,b]\), es decir,
\[\begin{equation} p_1=a_1+\frac{b_1-a_1}{2}=\frac{a_1+b_1}{2} \end{equation}\]
Si \(f(p_1)=0\), entonces \(p=p_1\) y terminamos.
Si \(f(p_1)\neq 0\), entonces \(f(p_1)\) tiene el mismo signo que ya sea \(f(a_1)\) o \(f(b_1)\).
Si \(f(p_1)\) y \(f(a_1)\) tienen el mismo signo, \(p\in (p_1, b_1)\). Sea \(a_2=p_1\) y \(b_2=b_1\).
Si \(f(p_1)\) y \(f(a_1)\) tienen signos opuestos, \(p\in (a_1, p_1)\). Sea \(a_2=a_1\) y \(b_2=p_1\).
Entonces, volvemos a aplicar el proceso al intervalo \([a_2, b_2]\).
Ejemplo 4.1 Demostrar que \(f(x)=x^3+4x^2-10=0\) tiene una raíz en \([1,2]\) y utiliza el método de la bisección para determinar una aproximación para la raíz que sea precisa por lo menos dentro de \(10^{-4}\)
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(1,2, 100)
y = x**3+4*x**2-10
plt.figure(figsize=(6,4.5))
plt.plot(x,y, color = "blue")
plt.axhline(0, color = "black", linewidth = 0.8, linestyle = "--") #Eje x
#plt.axvline(0, color = "black", linewidth = 0.8, linestyle = "--") #Eje y
# Personalizar gráfica
plt.xlabel("x")
plt.ylabel("y")
plt.grid()
plt.show()
Se podría seleccionar una tolerancia \(\varepsilon >0\) y generar \(p_1,\dots , p_N\) hasta que se cumpla una de las siguientes condiciones:
\[\begin{eqnarray} |p_N-p_{N-1}|<\varepsilon, &\\ \frac{|p_N-p_{N-1}|}{|p_N|}<\varepsilon,& \qquad p_N\neq 0\\ |f(p_N)|<\varepsilon& \end{eqnarray}\]
Se pueden presentar dificultades al utilizar algunos de estos criterios de tolerancia. Por ejemplo, existen sucesiones \(\{p_n\}_{n =0}^{\infty}\) con la propiedad de que las diferencias \(p_n − p_{n−1}\) convergen a cero mientras la sucesión misma diverge. También es posible que \(f(p_n)\) esté cerca de cero mientras \(p_n\) difiere significativamente de \(p\). Sin conocimiento adicional sobre \(f\) o \(p\), la segunda desigualdad es el mejor criterio de paro que se puede aplicar porque verifica el error relativo.
Teorema 4.1 Supongamos que \(f\in C[a,b]\) y \(f(a)\cdot f(b)<0\). El método de bisección genera una sucesión \(\{p_n\}_{n=1}^\infty\) que aproxima \(f(p)=0\) con
\[\begin{equation} |p_n-p|\leq \frac{b-a}{2^n}, \quad \text{ cuando } n\geq 1. \end{equation}\]
Al implementar el método en una computadora, necesitamos considerar los efectos del error de redondeo. Por ejemplo, el cálculo del punto medio del intervalo \([a_n, b_n]\) se debería encontrar a partir de la ecuación
\[\begin{equation} p_n=a_n+\frac{b_n-a_n}{2}\quad \text{ en lugar de } \quad p_n=\frac{a_n+b_n}{2} \end{equation}\]
Para determinar si el subintervalo de \([a_n, b_n]\) contiene una raíz de \(f\), es mejor utilizar la función signo, que se define como
\[\begin{equation} \operatorname{sign}(x)=\left\{\begin{array}{ll} -1, & \text { si } x<0 \\ 0, & \text { si } x=0 \\ 1, & \text { si } x>0 \end{array}\right. \end{equation}\]
En Python, scipy.optimize proporciona funciones para minimizar (o maximizar) funciones objetivo, posiblemente sujetas a restricciones. Incluye aproximación de soluciones para problemas no lineales (con soporte para algoritmos de optimización locales y globales), programación lineal, mínimos cuadrados restringidos y no lineales, búsqueda de raíces y ajuste de curvas.
En particular tiene implementada la rutina bisect para calcular la raíz de una función \(f\) por medio del método de bisección. Se sugiere revisar la documentación para consultar los parámetros y resultados de la función. A continuación se muestra el resultado para la función del Ejemplo 4.1.
from scipy import optimize
f = lambda x: x**3+4*x**2-10
root, info = optimize.bisect(f, 1, 2, full_output = True)
print("Información de la convergencia del método:")
print(info)
print("------------------------------------")
print(f"La raíz de la función es: {root}")Información de la convergencia del método:
converged: True
flag: converged
function_calls: 41
iterations: 39
root: 1.3652300134144753
method: bisect
------------------------------------
La raíz de la función es: 1.3652300134144753Supongamos que \(f\in C^2[a,b]\). Sea \(p_0\in [a,b]\) una aproximación para la raíz \(p\) de la función, de tal forma que \(f'(p_0)\neq 0\) y \(|p-p_0|\) es pequeño. Considérese el primer polinomio de Taylor para \(f(x)\) expandido alrededor de \(p_0\) y evaluado en \(x=p\) :
\[\begin{equation} f(p)=f(p_0)+(p-p_0)f'(p_0)+\frac{(p-p_0)^2}{2}f''(\xi(p)), \end{equation}\]
donde \(\xi(p)\) se encuentra entre \(p\) y \(p_0\). Puesto que \(f(p)=0\), de la última ecuación resulta
\[\begin{equation} 0=f(p_0)+(p-p_0)f'(p_0)+\frac{(p-p_0)^2}{2}f''(\xi(p)), \end{equation}\]
El método de Newton-Raphson se deriva al suponer que como \(|p-p_0|\) es pequeño, el término \((p-p_0)^2\) lo es aun más, entonces
\[\begin{equation} p\approx p_0-\frac{f(p_0)}{f'(p_0)}=p_1. \end{equation}\]
Esto constituye la base para el método de Newton, que empieza con una aproximación inicial \(p_0\) y genera la sucesión \(\{p_n\}_{n=0}^\infty\) mediante
\[\begin{equation} p_n=p_{n-1}-\frac{f(p_{n-1})}{f'(p_{n-1})}, \qquad \text{ para } n\geq 1. \end{equation}\]
Teorema 4.2 Sea \(f\in C^2[a,b]\). Si \(p\in (a,b)\) es tal que \(f(p)=0\) y \(f'(p)\neq 0\), entonces existe una \(\delta >0\) tal que el método de Newton genera una sucesión \(\{p_n\}_{n=1}^\infty\) que converge a \(p\) para cualquier aproximación inicial \(p_0\in[p-\delta, p+\delta]\).
Ejemplo 4.2 Consideremos la función \(f(x)=cos(x)-x\). Aproximar la raíz en el intervalo \([0, \pi/2]\) por medio del método de Newton-Raphson.
Inicialmente se muestra la gráfica de la función.
f = lambda x: np.cos(x)-x
x = np.linspace(- np.pi , np.pi, 100)
y = f(x)
plt.figure(figsize=(6,4.5))
plt.plot(x,y, color = "darkred")
plt.axhline(0, color = "black", linewidth = 0.8, linestyle = "--") #Eje x
plt.axvline(0, color = "black", linewidth = 0.8, linestyle = "--") #Eje y
# Personalizar gráfica
plt.xlabel("x")
plt.ylabel("y")
plt.grid()
plt.show()
En este caso la derivada de la función es \(f'(x)=-sen(x)-1\), la cual se anula para \(x=\pi/2+2n\pi\), \(n\in \mathbb{Z}\). Por lo tanto, podemos elegir como aproximación inicial \(p_0=0\). Utilizaremos la función newton de scipy.optimize, se sugiere consultar la respectiva documentación.
# Se brinda la derivada de la función
f_der = lambda x: - np.sin(x)-1
# Se aplica el método de Newton-Raphson
root_newton, info_newton = optimize.newton(f, 0, f_der, full_output = True)
print("Información de la convergencia del método de Newton-Raphson:")
print(info_newton)
print("------------------------------------")
print(f"La raíz de la función es: {root_newton}")Información de la convergencia del método de Newton-Raphson:
converged: True
flag: converged
function_calls: 10
iterations: 5
root: 0.7390851332151607
method: newton
------------------------------------
La raíz de la función es: 0.7390851332151607El método de la secante.
El método de Newton tiene una debilidad importante: la necesidad de conocer el valor de la derivada de \(f\) en cada aproximación. Para evitar la evaluación de la derivada en el método de Newton se presenta una ligera variación. Por definición,
\[\begin{equation} f'(p_{n-1})=\lim_{x\rightarrow p_{n-1}}\frac{f(x)-f(p_{n-1})}{x-p_{n-1}} \end{equation}\]
Si \(p_{n-2}\) está cerca de \(p_{n-1}\) entonces
\[\begin{equation} f^{\prime}\left(p_{n-1}\right) \approx \frac{f\left(p_{n-2}\right)-f\left(p_{n-1}\right)}{p_{n-2}-p_{n-1}}=\frac{f\left(p_{n-1}\right)-f\left(p_{n-2}\right)}{p_{n-1}-p_{n-2}} \end{equation}\]
Usando esta aproximación para \(f'(p_{n-1})\) en la fórmula de Newton obtenemos
\[\begin{equation} p_n=p_{n-1}-\frac{f(p_{n-1})(p_{n-1}-p_{n-2})}{f(p_{n-1})-f(p_{n-2})} \end{equation}\]
Ejemplo 4.3 Aproximar la raíz de la función \(f(x)=cos(x)-x\) en el intervalo \([0, \pi/2]\) por medio del método de la secante. Comparar el resultado con el que se obtuvo utilizando Newton-Raphson.
También se utiliza la función newton de scipy.optimize, pero omitiendo la derivada de la función y brindando los valores de \(x_0=0\) y \(x_1=1\):
# Se aplica el método de la secante
root_secante, info_secante = optimize.newton(f, x0 = 0, x1= 1, full_output = True)
print("Información de la convergencia del método de la secante:")
print(info_secante)
print("------------------------------------")
print(f"La raíz de la función es: {root_secante}")Información de la convergencia del método de la secante:
converged: True
flag: converged
function_calls: 7
iterations: 6
root: 0.7390851332151607
method: secant
------------------------------------
La raíz de la función es: 0.7390851332151607Observamos que en ambos casos hubo convergencia, en el método de la secante fueron necesarias 6 iteraciones, mientras que con Newton-Raphson fueron 5. Además, con la tolerancia dada por defecto (tol=1.48e-08) ambas aproximaciones coinciden:
root_newton == root_secantenp.True_fsolvePara aproximar la raíz de una función de manera general se puede utilizar la optimize.solve de la biblioteca SciPy, la cual utiliza algoritmos basados en el método de Newton-Raphson para aproximar las soluciones, por ejemplo, para la función \(f(x)=cos(x)-x\), se tiene:
f = lambda x: np.cos(x)-x
root_fsolve = optimize.fsolve(f, 0)
print(f"La raíz de la función es: {root_fsolve[0]}")La raíz de la función es: 0.7390851332151607Con esta función también es posible aproximar la solución de sistemas de ecuaciones no lineales con varias variables sin la necesidad de brindar el Jacobiano de la función.
Ejemplo 4.4 Aproximar la solución del siguiente sistema de dos ecuaciones cuadráticas:
\[\begin{eqnarray} x^2 +y^2-4 & = & 0\\ x^2 -y-1 & = & 0\\ \end{eqnarray}\]
En este caso es posible graficar las curvas de nivel, para la primera ecuación es una circunferencia de radio 2 con centro en el origen; para la segunda es una parábola:
# Valores para la ecuación 1
radio = 2
theta = np.linspace(0, 2*np.pi, 500) # Valores del ángulo
x1= radio * np.cos(theta)
y1= radio * np.sin(theta)
# Valores para la ecuación 2
x2= np.linspace(-2, 2, 500)
y2= x2 ** 2 - 1
# Gráfica de las curvas de nivel
plt.figure(figsize=(6, 6))
plt.plot(x1,y1, label = "Ecuación 1", color = "darkred")
plt.plot(x2,y2, label = "Ecuación 2", color = "steelblue")
plt.axhline(0, color = "black", linewidth = 0.8, linestyle = "--") #Eje x
plt.axvline(0, color = "black", linewidth = 0.8, linestyle = "--") #Eje y
plt.gca().set_aspect("equal", adjustable="box") # Escala igual en ambos ejes
# Personalizar gráfica
plt.xlabel("x")
plt.ylabel("y")
plt.grid()
plt.legend()
plt.show()
Observamos que existen dos soluciones (los puntos de intersección de las curvas de nivel) y en este caso, la solución a la cual converja el algoritmo dependerá de la estimación inicial:
def sistema(variables): # variables debe ser una lista de dos elementos
x, y = variables
return [x**2 + y**2 -4, x**2-y-1]
solucion_1 = optimize.fsolve(sistema, x0=[0.5,0])
solucion_2 = optimize.fsolve(sistema, x0=[-0.5,0])
print(f"Las soluciones del sistema son: ({round(solucion_1[0], 4)}, {round(solucion_1[1], 4)}) y ({round(solucion_2[0], 4)}, {round(solucion_2[1], 4)})")Las soluciones del sistema son: (1.5175, 1.3028) y (-1.5175, 1.3028)A continuación se muestran las curvas de nivel con los puntos de intersección:
# Valores para la ecuación 1
radio = 2
theta = np.linspace(0, 2*np.pi, 500) # Valores del ángulo
x1= radio * np.cos(theta)
y1= radio * np.sin(theta)
# Valores para la ecuación 2
x2= np.linspace(-2, 2, 500)
y2= x2 ** 2 - 1
# Gráfica de las curvas de nivel
plt.figure(figsize=(6, 6))
plt.plot(x1,y1, label = "Ecuación 1", color = "darkred")
plt.plot(x2,y2, label = "Ecuación 2", color = "steelblue")
plt.plot(solucion_1[0], solucion_1[1], marker = 'o', ms=8, color= "darkorange")
plt.plot(solucion_2[0], solucion_2[1], marker = 'o', ms=8, color= "darkorange")
plt.axhline(0, color = "black", linewidth = 0.8, linestyle = "--") #Eje x
plt.axvline(0, color = "black", linewidth = 0.8, linestyle = "--") #Eje y
plt.gca().set_aspect("equal", adjustable="box") # Escala igual en ambos ejes
# Personalizar gráfica
plt.xlabel("x")
plt.ylabel("y")
plt.grid()
plt.legend()
plt.show()