Recursive Function Bisector

One-dimensional function bisector

Motivation

I often deal with functions that vary rapidly at certain inputs but are largely smooth otherwise. For example, this is the spectral reflection response from a device known as a Guided Mode Resonance Filter:

If we wanted to accurately resolve this resonance with linearly spaced values of wavelength we would need roughly 1500 points. Approximately 95% of these wavelengths don't require nearly this high of a resolution, though. Suppose also that the function to calculate the response takes 60 seconds for every wavelength. That's a lot of wasted computation time on smooth function regions!

The bisect() function aims to decrease the number of required samples but still accurately resolve rapidly varying features in the supplied function. Instead of using ~1500 points, bisect() allowed me to plot this spectrum with only 181 points!

Who can benifit from this function?

If you are someone who:

• deals with computationally expensive functions and want to get away with fewer points,
• has rapidly changing functions where you need to resolve certain features with fewer points than others,
• doesn't know if and where the function values rapidly change, or
• Just tired of finding yourself saying "I guess I needed 2*N points to resolve this function rather than N points..."

this function is definitely for you!

How do I use the function?

The only input bisect() needs is:

1. A handle to the function (func) of interest.
2. Lower and upper bounds on the input to func.
3. Maximum alloweable tolerance
4. Optional maximum number of function evaluations maxEvals

The output of the function is, in general, a non-linearly spaced array of inputs (x) and their corresponding function values (f) that satisfy f = func(x).

How does the function work?

The bisect() function works recursively, which means the function calls itself until some criteria is met. Given a lower (lb) and upper (ub) bound and a tolerance (tol) on func, it:

• Calculates the midpoint (xm) between lb and ub
• Calculates the function value (fm) at xm
• Checks if the difference between fm and interpolated function value against tol. If it is within tolerance it returns. If not, it recalls bisect() to the left and right of the midpoint where these steps are repeated.

The graphic below shows an iteration where the difference between fm and interpolated function value is greater than tol requiring new midpoints to the left and right

The blue line represents func(x), the red line represents the interpolation between the upper and lower bounds, and the green points represent new function values to be ran through the bisect() algorithm at another iteration.

Example Usage

Suppose we want to resolve a sine wave on the interval where adjacent function values don't vary by more than 0.05. Implementing this would look something like:

% BISECT OPTIONS
opts.lb  = 0;
opts.ub  = 2*pi;
opts.tol = 0.05;

% BISECT FUNCTION HANDLE
func = @(x) sin(x);

% CALL BISECT()
[x,f] = bisect(func,opts);

Plotting the result gives us:

Well... that's obviously not correct! It's because the function value at the midpoint was exactly equal to the interpolation between the upper and lower bounds. This issue is resolved by setting minIter to something higher, which would look something like:

% BISECT OPTIONS
opts.lb      = 0;
opts.ub      = 2*pi;
opts.tol     = 0.05;
opts.minIter = 4;

% BISECT FUNCTION HANDLE
func = @(x) sin(x);

% CALL BISECT()
[x,f] = bisect(func,opts);

Now we get:

which is what we expected the first time around.

Some parting advice

1. The bisect() function can produce erroneous results if the midpoint function value and interpolated value are equivalent or less than the set tolerance, similar to the example show in the previous section. To fix this, adjust minEvals to something that will avoid this problem.

2. If your function values span several orders of magnitude, consider using

        func = @(x) log10(func(x));

and adjust tol as needed to avoid using more points than necessary.