Root finding using Python closures

roots.py

"""
Root finder as a closure.  Has benefits over the class 
approach because you are not carrying around self.
This demo is not really for the root finder.  It's mainly 
to demonstrate the use of closures in solving this kind 
of problem.  So, I don't implement all the methods.  
Also, the methods that are implemented don't have all 
the safeties.  Things like checking for zero slope or 
max iteration limit reached.
Bruce Wernick
13 August 2025
"""


def rootfind(**kwargs):
  """kwargs is optional.  This allows you to 
  control the method, iteration limit and tolerance.
  """  
  # unpack kwargs, using defaults if not supplied
  method = kwargs.get("method", "newton")
  maxits = kwargs.get("maxits", 48)
  tol = kwargs.get("tol", 1e-6)
  dx = kwargs.get("dx", 1e-3)  # for starting step

  
  def newton(f, x):
    """f must return f(x) and slope(x)"""
    for its in range(maxits):
      fx, dfx = f(x)  # note fx and dfx
      x -= fx / dfx  # newton step
      if abs(fx) <= tol:
        return x
    return x

  def secant(f, x):
    fx = f(x)
    x1 = x + dx  # small step
    f1 = f(x1)
    for its in range(maxits):
      x2 = x1 - f1 * (x1 - x) / (f1 - fx)  # secant step
      f2 = f(x2)
      if abs(f2) <= tol:
        return x2
      x, x1 = x1, x2
      fx, f1 = f1, f2
    return x2

  def falsi(f, x):
    """x must be a tuple that brackets the root"""
    ...

  def bisect(f, x):
    """x must be a tuple that brackets the root"""
    a, b = x
    fa = f(a)
    fb = f(b)
    for its in range(maxits):
      x = 0.5*(a + b)  # bisection step
      fx = f(x)
      if fb*fx <= 0:
        a, fa = x, fx
      else:
        b, fb = x, fx
      if abs(a-b) <= tol*(1+x):
        return x
    return x

  def broydn(f, x):
    """Broyden's method reduced to 1D.
    f returns only f(x) and x is a single guess
    """
    ...

  def brent(f, x):
    """Brent's method, generally recommended as a default.
    This very efficient version is by Kiusalaas "Numerical 
    methods in Engineering with Python".
    f returns f(x) and x is a tuple that brackets the root.
    """ 
    ...

  def iqi(f, x):
    """Inverse Quadratic Interpolation.
    f = f(x)
    x = tuple(a, b) that brackets the root
    """
    ...

  # The outer function will return the method you want.
  match method:
    case "newton": return newton
    case "secant": return secant
    case "falsi": return falsi
    case "bisect": return bisect
    case "broydn": return broydn
    case "brent": return brent
    case "iqi": return iqi
    case _: return broydn


# Examples
# Notice that the froot is just getting a
# solver from the rootfind closure.

# Example 1. Solve with Newton's method.
froot = rootfind(method="newton")
f = lambda x: ((x-3)*(x+5), 2*(x+1))
x = froot(f, 1)
print(f"Newton: root = {x:0.6g}")

# Example 2. Solve with secant method.
froot = rootfind(method="secant")
f = lambda x: (x-3)*(x+5)
x = froot(f, 1)
print(f"Secant: root = {x:0.6g}")

# Example 3. Solve with bisection method.
froot = rootfind(method="bisect")
f = lambda x: (x-3)*(x+5)
x = froot(f, (1,3.5))
print(f"Bisect: root = {x:0.6g}")

In another gist, I experimented with the Newton method as a decorator to a function. It works but doesn't really offer a benefit over a simple function approach. If there are a lot of methods, the class approach might seem efficient. But the self reference to all the common variables is a bit messy. This idea is like the class without the self. It looks very readable and efficient. I think I will use this idea.

Another good use of a Python closure would be for an interpolation function.
Given data points xdata=[] and ydata=[].
You could call the interpolator:
f = interp(xdata, ydata)
then simply use f(x) to get any value from the interpolator.

The more I see this, the more I like it.