Uno de estos algoritmos es el método de Newton para intervalos. Esta es una versión del algoritmo estándar de Newton (o Newton-Raphson), el cual es un método iterativo para encontrar raíces (ceros) de funciones. La versión para invervalos, sin embargo, es fundamentalmente diferente que la orginal, en el hecho de que puede (bajo las mejores circustancias) proveer garantías riguross acerca de la presencia o ausencia de raíces únicas de una función dada en un intervalo dado, y nos dice explícitamente cuando no puede proveer dicha garantía.
La idea del método de Newton es calcular una raíz $x^*$ de una función $f$ (i.e., un valor tal que $f(x^*) = 0$) dada una estimación inicial $x$ usando
$$ x^* = x - \frac{f(x)}{f'(\xi)}, $$para algún $\xi$ entre $x$ y $x^*$. Debido a que $\xi$ es desconocido, podemos limitarlo como
$$ f'(\xi) \in F'(X), $$donde $X$ es un intervalo que lo contiene, y $F'(X)$ denota el itervalo de extensión de la función $f$, que consiste en la aplicación de las mismas operaciones como la función $f$ al intervalo $X$.
Definimos un operador de intervalo de Newton $\mathcal{N}$ como:
$$ \mathcal{N}(X) := m(X) - \frac{F(m(X))}{F'(X)}, $$donde $m(X)$ es el punto medio de $X$ convertido en un intervalo.
Resulta que $\mathcal{N}$ nos dice precisamente si existe una raiz de $f$ en el intervalo $X$: no hay raíz si $\mathcal{N}(X) \bigcap X = \emptyset$, y hay una raíz única si $\mathcal{N}(X) \subseteq X$. También hay una extensión a intervalos en los que da derivada $F'(X)$ contiene el 0, en este caso el operador de Newton retorna la unión de dos intervalos.
Iterando el operador de Newton en el conjunto resultante nos da un algoritmo riguroso que garantiza encontrar todas las raíces de una función real en un intervalo dado; en realidad, lo hace si la raíz es única, pero no puede hacerlo si tenemos una raíz múltiple.
In [ ]:
using ValidatedNumerics
In [ ]:
f(x) = x^2 - 2
In [ ]:
newton(f, @interval(-5,5))
A la función newton se le pasa la función y el intervalo en que buscará las raíces; retorna un arreglo de objetos Root, que contienen el intervalo donde la raíz fue encontrada, junto con un símbolo :unique si está garantizado que la raíz es única en ese intervalo, o :unknown si el método no puede garantizarlo.
La función newton en algunos casos puede devolver un resultado en que algunas raíces tengan el símbolo :unique y otras :unknown. Un ejemplo interesante de esto se puede ver con los polimonios de Wilkison, que recordamos que tienen la siguiente forma
Usaremos el paquete Polynomials para encontrar los coeficientes:
In [ ]:
using Polynomials
function wilkinson_coefficients(n)
p = poly(collect(1:n))
p.a
end
function generate_wilkinson_horner(n)
coeffs = wilkinson_coefficients(n)
expr = :($(coeffs[end])) # comenzar con la potencia más alta
for i in length(coeffs)-1:-1:1
expr = :($expr*x + $(coeffs[i]))
end
expr
end
function subscriptify(n::Int)
subscript_digits = [c for c in "₀₁₂₃₄₅₆₇₈₉"]
dig = reverse(digits(n))
join([subscript_digits[i+1] for i in dig])
end
for n in 1:15
fn_name = symbol(string("W", subscriptify(n)))
expr = generate_wilkinson_horner(n)
@eval $(fn_name)(x) = $(expr)
end
Por ejemplo, $W_2 = 2 - 3x + x^2$,
In [ ]:
wilkinson_coefficients(2)
In [ ]:
generate_wilkinson_horner(2)
Veamos entonces un ejemplo intersante de lo que hemos comentado:
In [ ]:
wilkinson_coefficients(12)
In [ ]:
generate_wilkinson_horner(12)
In [ ]:
roots = newton(W₁₂, @interval(-10, 10))
Vemos que en algunos casos la función newton como habíamos dicho no nos dice si la raíz es única en algunos casos, y a veces repite el mismo valor de intervalo para el cual no encontró la unicidad. En algunos casos se le puede dar la vuelta a este problema aumentado la precisión de los intervalos, y utilizando el vector de raíces conocidas, para refinarlas con más precisión, y mejorando el performance del paquete (velocidad):
In [ ]:
set_interval_precision(256)
roots2 = newton(W₁₂, roots) #ver la defición de roots2