commit 7b2a9a6d4aaa30aacbe7139d10bc266799e3dd47 Author: Jérôme Date: Wed Apr 3 08:52:54 2019 +0200 Initial commit. TOMS748 1st algorithm coded following strictly the paper. diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..079dad0 --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +*.pdf +__pycache__ diff --git a/findEnclosingBracket.py b/findEnclosingBracket.py new file mode 100644 index 0000000..283f2ff --- /dev/null +++ b/findEnclosingBracket.py @@ -0,0 +1,50 @@ +# method to find a bracket that encloses a root of f +from math import * + +def findEnclosingBracket(f, x, factor): + """ Finds a bracket [a b] that encloses a root of f : f(a)*f(b) < 0. + A typical value for the factor is 2. + The initial guess must be on the correct side of the real line : the algorithm will never cross the zero during the search. + The factor is increased every 10 iterations in order to speed up convergence for when the root is orders of magnitude away from the initial guess. + If the bracket returned by the function is too wide, try reducing the factor (while always keeping it > 1). + """ + a = x + b = x*factor # try higher + fa = f(a) + fb = f(b) + + i = 0 + while i in range(100) and fa*fb > 0: + if fabs(fa) < fabs(fb): # fa is closer to 0 than fb, extend the interval in the a direction (reducing it) + a = a / factor + fa = f(a) + else: # fb is closer to 0 than fa, extend the interval in the b direction (augmenting it) + b = b * factor + fb = f(b) + i = i+1 + + if i % 10 == 0: # every 10 iterations, bump up the factor to speed up convergence + factor = factor * 5 + + return (a, b, fa, fb) + +f = lambda x : exp(sin(x))-2 +a, b, fa, fb = findEnclosingBracket(f, 0.005, 2) +print((a,b)); print((fa, fb)) + +f = lambda x : exp(-x)-1e-200 +a, b, fa, fb = findEnclosingBracket(f, 0.001, 2) +print((a,b)); print((fa, fb)) + +a, b, fa, fb = findEnclosingBracket(f, 0.001, 5) +print((a,b)); print((fa, fb)) + +f = lambda x : x**(1/3)-1000 +a, b, fa, fb = findEnclosingBracket(f, 1e-10, 2) +print((a,b)); print((fa, fb)) + +# from toms748 import TOMS748_solve1 +# def f(x): +# print(x) +# return x**(1/3)-1000 +# TOMS748_solve1(f,a,b) diff --git a/toms748.py b/toms748.py new file mode 100644 index 0000000..6520219 --- /dev/null +++ b/toms748.py @@ -0,0 +1,141 @@ +# Implementation of TOMS748 +# https://na.math.kit.edu/alefeld/download/1995_Algorithm_748_Enclosing_Zeros_of_Continuous_Functions.pdf + +from math import sin, fabs, cos + +# finite differences 1st derivative given a, b, fa, fb +def fbracket1(a, b, fa, fb): + return (fb-fa)/(b-a) + +# finite differences 2nd derivative given a, b, d, fa, fb, fd +def fbracket2(a, b, d, fa, fb, fd): + return (fbracket1(b, d, fb, fd) - fbracket1(a, b, fa, fb))/(d-a) + +# standard bracketing routine +# returns ahat, bhat, d +# d is a point outside the new interval +def bracket(a, b, c, fa, fb, fc): + # print('bracket') + if fabs(fc) < 1e-12: + print("root found : %f ; f(root) = %f" % (c, fc)) + if fa*fc < 0: + return (a, c, b, fa, fc, fb) + else: + return (c, b, a, fc, fb, fa) + +# finds an approximate solution to the quadratic P(x) = fa + f[a,b]*(x-a) + f[a,b,d]*(x-a)(x-b) +# with f[a,b] = fbracket1(a,b) and f[a,b,d] = fbracket2(a,b,d) +def NewtonQuadratic(a, b, d, fa, fb, fd, k): + A = fbracket2(a,b,d,fa,fb,fd) + B = fbracket1(a,b,fa,fb) + if A == 0: + return a - fa/B + if A*fa > 0: + r = a + else: + r = b + + for i in range(k): + r = r - B*r/(B + A*(2*r - a - b)) + return r + +# Inverse cubic interpolation evaluated at 0 (modified Aitken-Neville interpolation) +def ipzero(a, b, c, d, fa, fb, fc, fd): + Q11 = (c-d)*fc/(fd-fc) + Q21 = (b-c)*fb/(fc-fb) + Q31 = (a-b)*fa/(fb-fa) + D21 = (b-c)*fc/(fc-fb) + D31 = (a-b)*fb/(fb-fa) + Q22 = (D21-Q11)*fb/(fd-fb) + Q32 = (D31-Q21)*fa/(fc-fa) + D32 = (D31-Q21)*fc/(fc-fa) + Q33 = (D32-Q22)*fa/(fd-fa) + return a + Q31 + Q32 + Q33 + +def TOMS748_solve1(f, a, b): + mu = 0.5 + fa, fb = f(a), f(b) + c = a - fa/fbracket1(a, b, fa, fb) # 4.1.1 secant method + fc = f(c) + a, b, d, fa, fb, fd = bracket(a, b, c, fa, fb, fc) # 4.1.2 + e = d + fe = fd + # --- + for n in range(2,10): # 4.1.3 + if n == 2 or (fa == fb or fb == fd or fd == fe): + c = NewtonQuadratic(a, b, d, fa, fb, fd, 2) + else: + c = ipzero(a, b, d, e, fa, fb, fd, fe) + if (c-a)*(c-b) >= 0: + c = NewtonQuadratic(a, b, d, fa, fb, fd, 2) + # --- + fc = f(c) + a, b, dbar, fa, fb, fdbar = bracket(a, b, c, fa, fb, fc) # 4.1.4 + # --- + if fabs(fa) < fabs(fb): # 4.1.5 + u = a + fu = fa + else: + u = b + fu = fb + # --- + c = u - 2*fu/fbracket1(a, b, fa, fb) # 4.1.6 + # --- + if fabs(c - u) > 0.5*(b - a): # 4.1.7 + c = 0.5*(b + a) + # --- + fc = f(c) + a, b, dhat, fa, fb, fdhat = bracket(a, b, c, fa, fb, fc) # 4.1.8 + # --- + if b - a < mu*(b - a): # 4.1.9 + d = dhat + e = dbar + fd = fdhat + fe = fdbar + else: + e = dhat + fe = fdhat + c = 0.5*(a+b) + fc = f(c) + a, b, d, fa, fb, fd = bracket(a, b, c, fa, fb, fc) + +# TEST +f = lambda E : E - 0.5*sin(E) - 0.3 +a, b, c, d = 0, 1, 0.3, 1.2 +fa, fb, fc, fd = f(a), f(b), f(c), f(d) + +print(fbracket1(a, b, fa, fb)) +print(fbracket2(a, b, d, fa, fb, fd)) +print(bracket(a, b, c, fa, fb, fc)) +print(NewtonQuadratic(a, b, d, fa, fb, fd, 2)) +print(ipzero(a, b, c, d, fa, fb, fc, fd)) + +# solution using bracket only : +print('bracket only') +for i in range(50): + c = (a+b)/2 + fc = f(c) + a, b, d, fa, fb, fd = bracket(a, b, c, fa, fb, fc) + if fabs(b-a) < 1e-12: + break + print(i) + +print((a,b), (a+b)/2, b-a) + +print('Newton with derivative') +df = lambda E : 1 - 0.5*cos(E) +x = 0.5 +for i in range(50): + delta = f(x)/df(x) + x = x - delta + if fabs(delta) < 1e-12: + break + print(i) +print(x) + +print('TOMS748_solve1 only') +def fct1(E): + print("x = %20.15f" % E) + return E - 0.5*sin(E) - 0.3 + +TOMS748_solve1(fct1, 0, 1)