```1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 ``` ```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 ```
 ```1 2 3 ``` ```Secant 0.739085133215 0.0 Regula Falsi 0.739085133215 0.0 Dekker 0.739085133215 0.0 ```