[ create a new paste ] login | about

Link: http://codepad.org/Xe17Hcuf    [ raw code | output | fork ]

ploffie - Python, pasted on Mar 18:
from math import cos
from time import clock

EPS = 2.22e-16

def regfalsi(fun, x1, x2, args=(), xtol=1e-12, maxiter=1000):
    """ not so robust - also not fast (linear)"""
    a, b   = float(x1), float(x2)
    fa, fb = fun(a, *args), fun(b, *args)
    ncalls = 2
    if fa * fb > 0.0:
        msg = "No bracket: f(%f)=%f, f(%f)=%f" % (a, fa, b, fb)
        raise ValueError(msg)
    for i in xrange(maxiter):
        dx = -(a - b) / (1 - fb / fa)
        x = a + dx
        y = fun(x, *args)
        ncalls += 1
        if abs(dx) < xtol or y == 0.0:
            return x
        a, fa, b, fb = (x, y, b, fb) if y * fa > 0 else (a, fa, x, y)
    raise ValueError("too many iterations")

def secant(fun, x1, x2, args=(), xtol=1e-12, maxiter=1000):
    """the fastest (superlinear), but will often fail"""
    a, b   = float(x1), float(x2)
    fa, fb = fun(a, *args), fun(b, *args)
    for i in xrange(maxiter):
        if fa == fb:
            raise ZeroDivisionError("fa == fb in secant")
        dx = -(a - b) / (1 - fb / fa)
        x = a + dx
        y = fun(x, *args)
        if abs(dx) < xtol or y == 0.0:
            return x
        a, fa, b, fb = b, fb, x, y
    raise ValueError("too many iterations")

def dekker(fun, x1, x2, args=(), xtol=1e-8, disp=False):
    """see wiki Brent
       c is last iterate of b
       a, b bracket solution abs(fb) < abs(fa)
       if f(s) has same sign as fa, then old b beomes new contrapoint
       the solution is within 6*EPS + 2 * xtol
    """
    a, b   = float(x1), float(x2)
    fa, fb = fun(a, *args), fun(b, *args)
    if fa * fb > 0.0: raise ValueError("f(%f)=%f - f(%f)=%f" % (a, b, fa, fb))
    c, fc = a, fa
    feval = 2
    while 1:
        if abs(fa) < abs(fb):
            a, fa, b, fb, c, fc = b, fb, a, fa, b, fb
        xm = (a - b) / 2  # bisection step
        tol1 = 2 * EPS * abs(b) + xtol
        if abs(xm) <= tol1 or fb == 0:
            return b
        # now calculate secant step - use points b and c
        s = fb / fc
        p = (c - b) * s
        q = 1 - s
        if p > 0:
            q = -q
        p = abs(p)
        d = p / q if p < (q * xm) else xm
        c, fc = b, fb
        # do not take steps smaller than tol1
        b += (d if abs(d) > tol1 else sign(tol1, xm))
        fb = fun(b, *args)
        feval += 1
        if disp:
            msg = "b" if abs(d) == abs(xm) else "s"
            fmt = "%3d a=%7.4g b=%12.9g c=%7.4g fb=%10.3g xm=%8.2g %s"
            print fmt % (feval, a, b, c, fb, xm, msg)
        if fa * fb > 0:     # if new iterate has same sign as a - old b take s place of a
            a, fa = c, fc
    
def sign(x, y):
    if y > 0.0  :
        return x
    elif y < 0.0:
        return -x
    else:
        return 0.0
        
def dirty(x):
    """zero close to 1.0"""
    D = 1e-2
    epsi = 1e-10
    if x > 1 - D:
        res = ((x - (1 - D)) / D) ** 2
    else:
        res = 0
    return res - epsi
    
def f(x):
    return cos(x) - x
    
disp = False
precision = 12
_xtol = 10 ** -(precision + 1)
fmt = "%20." + str(precision) + "f"

fun, xmin, xmax = dirty, 0, 1
fun, xmin, xmax = f, 0, 1
    
t0 = clock()
print "Secant      ",
try:
    print fmt % secant(fun, xmin, xmax, xtol=_xtol), clock() - t0
except ValueError, e:
    print e
except ZeroDivisionError, e2:
    print e2
t0 = clock()
print "Regula Falsi",
try:
    print fmt % regfalsi(fun, xmin, xmax, xtol=_xtol, maxiter=1000000), clock() - t0
except ValueError, e:
    print e
t0 = clock()
print "Dekker      ",
print fmt % dekker(fun, xmin, xmax, xtol=_xtol, disp=disp), clock() - t0


Output:
1
2
3
Secant             0.739085133215 0.0
Regula Falsi       0.739085133215 0.0
Dekker             0.739085133215 0.0


Create a new paste based on this one


Comments: