quadinterp.py

#!/usr/bin/env python3
#
#   quadinterp.py  -  by Don Cross  -  https://github.com/cosinekitty
#
#   Sample code for finding approximate roots of
#   an expensive function f(t) using quadratic interpolation.
#   This source code is public domain. Use at your own risk.
#
import math

def QuadraticInterpolate(f, t0, t1):
    # Find the midpoint between t0 and t1:
    t2 = (t0 + t1) / 2

    # Evaluate the function at all three points.
    f0 = f(t0)
    f2 = f(t2)
    f1 = f(t1)

    # Compute the coefficients of the parabola p(x) that passes
    # through the three points (t0,f0), (t2,f2), (t1,f1).
    Q = (f1 + f0)/2 - f2
    R = (f1 - f0)/2
    S = f2

    if Q == 0:
        # Special case: the three points are collinear.
        if R == 0:
            # There is no root because the line is horizontal.
            return None
        # Solve for the place where the line passes through the x-axis.
        x = -S/R
        # Is the root within the search interval?
        if -1 <= x <= +1:
            # Convert x back to the original independent variable t:
            t = (t2 - t0)*x + t2
            return t
        # The x value was outside the search interval.
        # Therefore it is not valid.
        return None

    # The approximation curve is a parabola.
    # Calculate the radicand u. Then we can determine
    # how many real roots there are.
    u = R*R - 4*Q*S

    # We require a non-tangent real root.
    if u <= 0:
        # There are no real roots for the parabolic curve,
        # or there is a single tangent root.
        return None

    ru = math.sqrt(u)
    x1 = (-R + ru) / (2*Q)
    x2 = (-R - ru) / (2*Q)

    # Form a list of the x values that are within
    # the normalized interval [-1, +1].
    xlist = [x for x in [x1, x2] if -1 <= x <= +1]

    # There must be exactly one real root inside
    # the normalized interval in order for the
    # solution to be considered valid.
    if len(xlist) == 1:
        # Translate the normalized parameter x back
        # into the independent variable t:
        t = (t2 - t0)*xlist[0] + t2
        return t

    # Either there were no valid roots, or
    # there were 2 roots (interval was too large).
    return None

# Here is a sample function we will find an approximate root for.
def Func(t):
    return math.cos(t) - t

if __name__ == '__main__':
    # Use quadratic interpolation to find an approximate
    # solution to the equation cos(t) = t.

    # Suppose we know the root is inside a small interval (t0, t1).
    t0 = 0.7
    t1 = 0.8

    # Use the quadratic interpolator to find the approximate root t.
    t = QuadraticInterpolate(Func, t0, t1)

    if t is None:
        # This happens if the solver could not find exactly one
        # root inside the given interval (t0, t1).
        print('Could not find a unique root in the given interval.')
    else:
        # Here is the known correct solution:
        correct_root = 0.7390851332151607
        print('Found the approximate root:    {}'.format(t))
        print('The correct value of the root: {}'.format(correct_root))
        error = (t - correct_root) / correct_root
        print('The relative error is: {}'.format(error))