Reto de métodos numéricos: Día 2
Durante octubre (2017) estaré escribiendo un programa por día para algunos métodos numéricos famosos en Python y Julia. Esto está pensado como un ejercicio, no esperen que el código sea lo suficientemente bueno para usarse en la "vida real". Además, también debo mencionar que casi que no tengo experiencia con Julia, así que probablemente no escriba un Julia idiomático y se parezca más a Python.
Regula falsi
El segundo método a considerar es el método de posición falsa, o regula falsi. Este método se utiliza para resolver la ecuación \(f(x) = 0\) para \(x\) real, y \(f\) continua. Se comienza con un intervalo \([a, b]\), donde \(f (a)\) y \(f (b)\) deben tener signos opuestos. El método es similar al método de bisección pero en lugar de dividir a la mitad el original intervalo toma como nuevo punto del intervalo la intersección de la línea que conecta la función evaluada en los dos puntos extremos. Entonces, el nuevo el punto se calcula a partir de
Como en el método de bisección, mantenemos el intervalo donde la solución aparece, es decir, donde cambia el signo de la función.
Usaremos la función \(f(x) = \cos(x) - x\) para probar los códigos, y como intervalo inicial \([0.5, \pi/4]\).
A continuación se presentan los códigos.
Python
from __future__ import division, print_function from numpy import abs, cos, pi def regula_falsi(fun, a, b, niter=50, ftol=1e-12, verbose=False): if fun(a) * fun(b) > 0: c = None msg = "The function should have a sign change in the interval." else: for cont in range(niter): qa = fun(a) qb = fun(b) c = (a*qb - b*qa)/(qb - qa) qc = fun(c) if verbose: print("n: {}, c: {}".format(cont, c)) msg = "Maximum number of iterations reached." if abs(qc) < ftol: msg = "Root found with desired accuracy." break elif qa * qc < 0: b = c elif qb * qc < 0: a = c return c, msg def fun(x): return cos(x) - x print(regula_falsi(fun, 0.5, 0.25*pi))
Julia
function regula_falsi(fun, a, b, niter=50, ftol=1e-12, verbose=false) if fun(a) * fun(b) > 0 c = nothing msg = "The function should have a sign change in the interval." else for cont = 1:niter qa = fun(a) qb = fun(b) c = (a*qb - b*qa)/(qb - qa) qc = fun(c) if verbose println("n: $(cont), c: $(c)") end if abs(fun(c)) < ftol msg = "Root found with desired accuracy." break elseif qa * qc < 0 b = c elseif qb * qc < 0 a = c end msg = "Maximum number of iterations reached." end end return c, msg end function fun(x) return cos(x) - x end println(regula_falsi(fun, 0.5, 0.25*pi))
Comparación
Respecto al número de líneas de código tenemos: 29 en Python y 32 en Julia.
La comparación en tiempo de ejecución se realiza con el comando mágico
de IPython %timeit y con @benchmark en Julia.
Para Python:
con resultado
Para Julia:
con resultado
