Find first root of a black box function, or any ne

Question

I have a black box function, f(x) and a range of values for x. I need to find the lowest value of x for which f(x) = 0.

I know that for the start of the range of x, f(x) > 0, and if I had a value for which f(x) < 0 I could use regula falsi, or similar root finding methods, to try determine f(x)=0.

I know f(x) is continuous, and should only have 0,1 or 2 roots for the range in question, but it might have a local minimum.

f(x) is somewhat computationally expensive, and I'll have to find this first root a lot.

I was thinking some kind of hill climbing with a degree of randomness to avoid any local minimums, but then how do you know if there was no minimum less than zero or if you just haven't found it yet? I think the function shouldn't have more than two minimum points, but I can't be absolutely certain of that enough to rely on it.

If it helps, x in this case represents a time, and f(x) represents the distance between a ship and a body in orbit (moon/planet) at that time. I need the first point where they are a certain distance from each other.

Answer 1

Your function has only 0, 1 or 2 roots, so it can be done using an algorithm it doesn't ensure the first root.

• Find one root using Newton's method or other method. If it can't find any root, this algorithm also give up.
• Let the found root is r and beginning of the range of the x is x0. let d = (r-x0)/2.
• While d > 0, calculate f(r-d). if f(r-d) > 0, half d (d := d / 2) and loop. if f(r-d) <= 0, escape the loop.
• if loop is finished by d = 0, report r as the first root. if d > 0, find a root between x0 and r-d by using any other method and report it.

I assumed two prerequiesite conditions.

1. f(x) takes x of floating point numbers
2. At each point of the roots of f(x), the graph of f(x) crosses to x-axis. They are not touching root like x=0 in f(x)=x^2.

Using condition 2, you can prove that if there is no point such that f(r-d) < 0, ∀ x: x0 < x < r, f(x) > 0.

Answer 2

My method will sound pretty complicated, but in the end the computation time of the method will be far smaller than the distance calculations (evaluation of your f(x)). Also, there are quite many implementations of it already written up in existing libraries.

So what I would do:

• approximate f(x) with a Chebychev polynomial
• find the real roots of that polynomial
• If any are found, use those roots as initial estimates in a more precise rootfinder (if needed)

Given the nature of your function (smooth, continuous, otherwise well-behaved) and the information that there's 0,1 or 2 roots, a good Chebychev polynomial can already be found with 3 evaluations of f(x).

Then find the eigenvalues of the companion matrix of the Chebychev coefficients; these correspond to the roots of the Chebychev polynomial.

If all are imaginary, there's 0 roots. If there are some real roots, check if two are equal (that "rare" case you spoke of). Otherwise, all real eigenvalues are roots; the lowest one of which is the root you seek.

Then use Newton-Raphson to refine (if necessary, or use a better Chebychev polynomial). Derivatives of f can be approximated using central differences

f'(x) = ( f(x+h)-f(h-x) ) /2/h     (for small h)

I have an implementation of the Chebychev routines in Matlab/Octave (given below). Use like this:

R = FindRealRoots(@f, x_min, x_max, 5, true,true);

with [x_min,x_max] your range in x, 5 the number of points to use for finding the polynomial (the higher, the more accurate. Equals the amount of function evaluations needed), and the last true will make a plot of the actual function and the Chebychev approximation to it (mainly for testing purposes).

Now, the implementation:

% FINDREALROOTS     Find approximations to all real roots of any function 
%                   on an interval [a, b].
%
% USAGE:
%    Roots = FindRealRoots(funfcn, a, b, n, vectorized, make_plot)
%
% FINDREALROOTS() approximates all the real roots of the function 'funfcn' 
% in the interval [a,b]. It does so by finding the roots of an [n]-th degree 
% Chebyshev polynomial approximation, via the eignevalues of the associated 
% companion matrix. 
%
% When the argument [vectorized] is [true], FINDREALROOTS() will evaluate 
% the function 'funfcn' at all [n] required points in the interval 
% simultaneously. Otherwise, it will use ARRAFUN() to calculate the [n] 
% function values one-by-one. [vectorized] defaults to [false]. 
%
% When the argument [make_plot] is true, FINDREALROOTS() plots the 
% original function and the Chebyshev approximation, and shows any roots on
% the given interval. Also [make_plot] defaults to [false]. 
%
% All [Roots] (if any) will be sorted.
% 
% First version 26th May 2007 by Stephen Morris, 
% Nightingale-EOS Ltd., St. Asaph, Wales.
%
% Modified 14/Nov (Rody Oldenhuis)
%
% See also roots, eig.

function Roots = FindRealRoots(funfcn, a, b, n, vectorized, make_plot)

    % parse input and initialize.
    inarg = nargin; 
    if n <= 2, n = 3; end                    % Minimum [n] is 3:    
    if (inarg < 5), vectorized = false; end  % default: function isn't vectorized
    if (inarg < 6), make_plot = false; end   % default: don't make plot

    % some convenient variables
    bma = (b-a)/2;  bpa = (b+a)/2;  Roots = [];

    % Obtain the Chebyshev coefficients for the function
    %
    % Based on the routine given in Numerical Recipes (3rd) section 5.8;
    % calculates the Chebyshev coefficients necessary to approximate some
    % function over the interval [a,b]

    % initialize 
    c = zeros(1,n);  k=(1:n)';  y = cos(pi*((1:n)-1/2)/n); 
    % evaluate function on Chebychev nodes
    if vectorized
        f = feval(funfcn,(y*bma)+bpa);
    else
        f = arrayfun(@(x) feval(funfcn,x),(y*bma)+bpa);
    end

    % compute the coefficients
    for j=1:n, c(j)=(f(:).'*(cos((pi*(j-1))*((k-0.5)/n))))*(2-(j==1))/n; end       

    % coefficients may be [NaN] if [inf]
    % ??? TODO - it is of course possible for c(n) to be zero...
    if any(~isfinite(c(:))) || (c(n) == 0), return; end

    % Define [A] as the Frobenius-Chebyshev companion matrix. This is based
    % on the form given by J.P. Boyd, Appl. Num. Math. 56 pp.1077-1091 (2006).
    one = ones(n-3,1);
    A = diag([one/2; 0],-1) + diag([1; one/2],+1);
    A(end, :) = -c(1:n-1)/2/c(n);
    A(end,end-1) = A(end,end-1) + 0.5;

    % Now we have the companion matrix, we can find its eigenvalues using the
    % MATLAB built-in function. We're only interested in the real elements of
    % the matrix:
    eigvals = eig(A);  realvals = eigvals(imag(eigvals)==0);

    % if there aren't any real roots, return
    if isempty(realvals), return; end

    % Of course these are the roots scaled to the canonical interval [-1,1]. We
    % need to map them back onto the interval [a, b]; we widen the interval just
    % a tiny bit to make sure that we don't miss any that are right on the 
    % boundaries.
    rangevals = nonzeros(realvals(abs(realvals) <= 1+1e-5));

    % also sort the roots
    Roots = sort(rangevals*bma + bpa);

    % As a sanity check we'll plot out the original function and its Chebyshev
    % approximation: if they don't match then we know to call the routine again
    % with a larger 'n'.
    if make_plot
        % simple grid
        grid = linspace(a,b, max(25,n));
        % evaluate function
        if vectorized
            fungrid = feval(funfcn, grid);
        else
            fungrid = arrayfun(@(x) feval(funfcn,x), grid);
        end        
        % corresponding Chebychev-grid (more complicated but essentially the same)
        y = (2.*grid-a-b)./(b-a); d = zeros(1,length(grid)); dd = d;
        for j = length(c):-1:2, sv=d; d=(2*y.*d)-dd+c(j); dd=sv; end, chebgrid=(y.*d)-dd+c(1);
        % Now make plot
        figure(1), clf,  hold on
        plot(grid, fungrid ,'color' , 'r');
        line(grid, chebgrid,'color' , 'b'); 
        line(grid, zeros(1,length(grid)), 'linestyle','--')
        legend('function', 'interpolation')
    end % make plot

end % FindRealRoots

Answer 3

You could use the secant method which is a discrete version of Newton's method.

The root is estimated by calculating the line between two points (= the secant) and its crossing of the X axis.

Answer 4

You could make a small change to the uniroot.all function from the R library rootSolve.

uniroot.all <- function (f, interval, lower= min(interval),
                         upper= max(interval), tol= .Machine$double.eps^0.2,
                         maxiter= 1000, n = 100, nroots = -1, ... ) {

  ## error checking as in uniroot...
  if (!missing(interval) && length(interval) != 2)
    stop("'interval' must be a vector of length 2")
  if (!is.numeric(lower) || !is.numeric(upper) || lower >=
      upper)
    stop("lower < upper  is not fulfilled")

  ## subdivide interval in n subintervals and estimate the function values
  xseq <- seq(lower,upper,len=n+1)
  mod  <- f(xseq,...)

  ## some function values may already be 0
  Equi <- xseq[which(mod==0)]

  ss   <- mod[1:n]*mod[2:(n+1)]  # interval where functionvalues change sign
  ii   <- which(ss<0)

  for (i in ii) {
    Equi <- c(Equi, uniroot(f, lower = xseq[i], upper = xseq[i+1] ,...)$root)
    if (length(Equi) == nroots) {
      return(Equi)
    }
  }
  return(Equi)
}

And run it like this:

uniroot.all(f = your_function, interval = c(start, stop), nroots = 1)