Los métodos numéricos validados tienen diversas aplicaciones, en teoría se podría escribir código en términos de artimética de intervalos para la mayoría, sino todos, de los actuales códigos numéricos. Se han creado librerías y paquetes en la mayoría de los lenguajes más importantes, los cuales están disponibles para su uso junto con una documentación detallada (en algunos casos).
Se han utilizado métodos numéricos validados para resolver y mejorar la solución a muchos problemas matemáticos, físicos y computacionales. Entre ellos:
En física se han utilizado para encontrar órbitas periódicas en sistemas dinámicos, integración y evolución de ecuaciones diferenciales ordinarias en sistemas hamiltonianos, entre muchos otros.
Hay un problema común que se encuentra en muchas aplicaciones físicas, y es encontrar raíces de funciones en un intervalo dado. Como hemos demostrado, estos métodos validados permiten encontrar estas raíces rigurosamente, garantizando que no "perderemos" ningún cero en un intervalo dado, y que nos diga también rigurosamente si hay un único cero en una región dada, lamentablemente aún no se logra un algoritmo que pueda asegurar la unicidad de los ceros para cualquier función.
A continuación les mostraré un código que he desarrollado para encontrar raíces a funciones de una variable (actualmente la librería solo trabaja para funciones de una variable) utilizando métodos numéricos validados.
La función roots_interval utiliza el método de Newton de intervalos que vimos anteriormente. Hay que probar la existencia y unicidad de las raíces en el método de Newton de intervalos, así como su estabilidad para que pueda considerarse un método numérico valido y serio. Esto fue hecho por los que plantearon este algoritmo originalmente y algunos otros matemáticos, pueden encontrarse en el libro de Ramon Moore y algunos artículos de Warwick Tucker.
En resumen, se utilizan teoremas de punto fijo para probar la existencia y unicidad. En caso del algoritmo de intervalos de Newton, la existencia está dada por el teorema de punto fijo de Brouwer que nos dice si un conjunto $B$ es homeomorfo a la bola cerrada unitaria en $\mathbb{R}^n$, entonces dado un mapeo continuo $f: B \rightarrow B$ entonces $\exists x \in B$ tal que $f(x) = x$. Esto aplica a la última parte del método de Newton para intervalos, en el operado de Newton $\mathcal{N}$ para verificar la existencia de la raíz. La unicidad viene de que $f'\ne 0$ utilizando los teoremas de Schauder y Banach.
Antes de mostrarle la función y los resultados que obtuvieron, primero veamos algunos ejemplos de funciones polinomiales importantes, que pueden obtenerse utilizando las librerías Polynomials.jl y Jacobi.jl.
In [ ]:
using ValidatedNumerics, Polynomials, Jacobi, PyPlot, PyCall
@pyimport matplotlib.patches as patches
In [ ]:
# Construir polinomio 1 + 2x
Poly([1,2])
In [ ]:
# Construir polinomio 1 + 2x + 10x^2 + 22x^3 + 4x^4
Poly([1,2,10,22,4])
In [ ]:
# Puedo escoger la variable
Poly([1,2,10,22,4], :y)
In [ ]:
# Construir polinomio con sus raíces
# Lo siguiente representa (x-1)*(x-2)*(x-3)
poly([1,2,3])
Jacobi.jlEsta librería permite hacer operaciones y trabajar con los polinomios de orden $n$ de Jacobi, Legendre y los polinomios de Chebyshev de primero y segundo tipo; así como evaluarlos y obtener sus raíces (de nuevo con métodos usuales). Además utilizando la librería Polynomials.jl permite obtener expresiones para todos estos polinomios.
In [ ]:
poly_legendre(4)
In [ ]:
poly_chebyshev(3)
In [ ]:
poly_legendre(20)
In [ ]:
poly_chebyshev(24)
In [ ]:
x = linspace(-1.0,1.0,401)
y1 = legendre(x,1)
y2 = legendre(x,2)
y3 = legendre(x,3)
y4 = legendre(x,4)
y5 = legendre(x,5)
y6 = legendre(x,6)
plot(x,y1,x,y2,x,y3,x,y4,x,y5,x,y6)
axhline(y=0.0, color="black")
La función roots_interval que he creado recibe una función de una variable, un límite inferior y uno superior para definir
el intervalo donde buscará los ceros. La función se le pasa como una función anónima (lambda function). Por adentro utiliza dos funciones de la librería ValidatedNumerics, uno es la ya vista newton que permite encontrar ceros para una función dada, y además utiliza la función find_roots_midpoint que ubica el punto medio del intervalo y da un valor concreto para el cero. Luego utilizando PyPlot grafica la función y sus ceros.
El cero que muestra la gráfica es el valor concreto del punto medio del invervalo donde se encuentra el cero, pero debe recordarse que se están utilizando métodos numéricos validados y por lo tanto nos interesa es el intervalo donde se encuentra el cero, ya que está garantizado que adentro del mismo está el cero real. La función también imprime el intervalo donde están los ceros.
In [ ]:
function roots_interval(fun,infBound,supBound;
steps=Int((supBound-infBound)*100),lineColor="b",zerosColor="rs")
try
x = linspace(infBound,supBound,steps)
result = fun(x)
ro = newton(fun,@interval(infBound,supBound))
z = find_roots_midpoint(fun,infBound,supBound)[1]
plot(x,result,lineColor)
plot(z,zeros(z),zerosColor)
axhline(y=0.0, color="black")
ro
catch err
if isa(err,BoundsError)
print("There are no zeros in the interval ","[",infBound,",",supBound,"]")
elseif isa(err,UndefVarError)
print("Function is not correctly written. Check you variable.")
elseif isa(err,MethodError)
print("Function is not correctly written. Check your methods.")
end
end
end
Ejemplo simple, ceros del polinomio de Legendre de cuarto orden
$$ P_4(x) = \frac{1}{8}(35x^4 - 30x^2 + 3) $$
In [ ]:
# Obtención del polimonio con Jacobi.jl
poly_legendre(4)
In [ ]:
roots_interval(x -> 0.375 - 3.75x.^2 + 4.375x.^4,-1,1)
La función valida que la función que se le ingresa esté bien escrita, y también le dice al usuario si no existen ceros en el intervalo que se le está pasando.
In [ ]:
roots_interval(x -> cos(x) + x.^4,-1,1)
In [ ]:
plot(x,cos(x) + x.^4,-1,1)
In [ ]:
roots_interval(x -> sin(2x) + cos(3x),-5,5)
In [ ]:
roots_interval(x -> x.^7 - 28x.^6 + 322x.^5
- 1960x.^4+6769x.^3-13132x.^2+13068x -5040,-5,10)
roots_interval_rectangleAunque hemos graficado el punto medio del intervalo, que nos da la raíz a una alta precisión, sería interesante poder graficar la función y los intervalos donde se encuentra la raíz, en vez de el punto medio de los mismos. Tomando en consideración esto, nace la función roots_interval_rectangle que permite hacer lo mismo que la función roots_interval pero grafica el intervalo donde se encuentra la raíz y no el punto intermedio donde está.
In [ ]:
function draw_rectangle(x, y, xwidth, yheight, color="green")
rectangle = patches.Rectangle
ax = gca()
ax[:add_patch](rectangle((x, y), xwidth, yheight, facecolor=color, alpha=0.5))
end
In [ ]:
function roots_interval_rectangle(fun,infBound,supBound;
steps=Int((supBound-infBound)*100),lineColor="b",zerosColor="rs")
try
z = linspace(infBound,supBound,steps)
result = fun(z)
ro = newton(fun,@interval(infBound,supBound))
roots = map(i -> (i.interval.lo,i.interval.hi),ro)
for r in roots
xwidth = fun(r[1])+0.07
x = r[1]-xwidth/2
y = - abs(fun(x))
yheight = abs(y)*1.8
draw_rectangle(x,y,xwidth,yheight)
end
plot(z,result,lineColor)
axhline(y=0.0, color="black")
ro
catch err
if isa(err,BoundsError)
print("There are no zeros in the interval ","[",infBound,",",supBound,"]")
elseif isa(err,UndefVarError)
print("Function is not correctly written. Check you variable.")
elseif isa(err,MethodError)
print("Function is not correctly written. Check your methods.")
end
end
end
In [ ]:
roots_interval_rectangle(x -> 0.375 - 3.75x.^2 + 4.375x.^4,-1,1)
In [ ]:
roots_interval(x -> ?? ,a,b)