HP50G Bisection methode

[Albert Chan]

(10-03-2020 02:03 AM)Albert Chan Wrote:  >>> f = lambda x: x**5 - 5
>>> solve_bisect(f,1,2)
1 2
1 1.5
1.25 1.5
1.375 1.5
1.375 1.4375
1.375 1.40625
1.375 1.390625
1.375 1.3828125
1.37890625 1.3828125
...

For ε=1e-15, solve_bisect() required log2(1/ε) ≈ 50 bisects (assuming it did not hit the root head on).

I just learned Ridder's method.
It has the advantage of bisection (guaranteed convergence), worst-case convergence of bisection.
But, if function is well behaved, each loop approximately doubled convergence.
(since it take 2 function calls to do 1 loop, convergence is √2 order)

All code before getting (c, fc) is the same as solve_bisect(), then we add point (d, fd):

Code:

def solve_ridder(f,a,b, eps=1e-8, verbal=True):
    fa, fb = f(a), f(b)
    if fa==0: return a
    if fb==0: return b
    if fa > fb: a, fa, b, fb = b, fb, a, fa
    assert fa < 0 and fb > 0, "bad interval"
    while abs(a-b) > eps:
        if verbal: print a, b
        c = (a+b)/2
        fc = f(c)
        if fc == 0: return c
        d = c - (c-a) * fc/(fc*fc-fa*fb)**0.5
        fd = f(d)
        if fd == 0: return d
        if fd < 0:
            a, fa = d, fd
            if fc > 0: b, fb = c, fc
        else:
            b, fb = d, fd
            if fc < 0: a, fa = c, fc
    return (a+b)/2

>>> f = lambda x: x**5 - 5
>>> solve_ridder(f, 1,2, eps=1e-15)
1 2
1.3789222674 1.5
1.37972953756 1.4394611337
1.37972966146 1.40959533563
1.37972966146 1.37972966146
1.37972966146 1.37972966146
1.3797296614612149