Code covered by the BSD License

### Highlights from Vectorized Bisection Search

4.875
4.9 | 8 ratings Rate this file 28 Downloads (last 30 days) File Size: 7.09 KB File ID: #28150 Version: 1.13

# Vectorized Bisection Search

### Sky Sartorius (view profile)

08 Jul 2010 (Updated )

Robust root-finding method that handles n-dim inputs and therefore has key advantages over FZERO.

File Information
Description

BISECTION is a fast and robust root-finding method that handles n-dim arrays.
[x,fVal,ExitFlag] = BISECTION(f,LB,UB,target,options) finds x +/- TolX (LB < x < UB) such that f(x) = target +/- TolFun.
x = BISECTION(f,LB,UB) finds the root(s) of function f on the interval [LB, UB], i.e. finds x such that f(x) = 0 where LB <= x <= UB. f will never be evaluated outside of the interval specified by LB and UB. f should have only one root and f(UB) and f(LB) must bound it. Elements of x are NaN for instances where a solution could not be found.
x = BISECTION(f,LB,UB,target) finds x such that f(x) = target.

x = BISECTION(f,LB,UB,target,TolX) will terminate the search when the search interval is smaller than TolX (TolX must be positive).
x = BISECTION(f,LB,UB,target,options) solves with the default parameters replaced by values in the structure OPTIONS, an argument created with the OPTIMSET function. Used options are TolX and TolFun. Note that OPTIMSET will not allow arrays for tolerances, so set the fields of the options structure manually for non-scalar TolX or TolFun.
[x,fVal] = BISECTION(f,...) returns the value of f evaluated at x.

[x,fVal,ExitFlag] = BISECTION(...) returns an ExitFlag that describes the exit condition of BISECTION. Possible values of elements of ExitFlag and the corresponding exit conditions are

1 Search interval smaller than TolX.
2 Function value within TolFun of target.
3 Search interval smaller than TolX AND function value within TolFun of target.
-1 No solution found.

Any or all of f(scalar), f(array), LB, UB, target, TolX, or TolFun may be scalar or n-dim arrays. All non-scalar arrays must be the same size. All outputs will be this size.

Default values are target = 0, TolX = 1e-6, and TolFun = 0.

There is no iteration limit. This is because BISECTION (with a TolX that won't introduce numerical issues) is guaranteed to converge if f is a continuous function on the interval [UB, LB] and f(x)-target changes sign on the interval.

The bisection method is a very robust root-finding method. The absolute error is halved at each step so the method converges linearly. However, Brent's method (such as implemented in FZERO) can converge superlinearly and is as robust. FZERO also has more features and input checking, so use BISECTION in cases where FZERO would have to be implemented in a loop to solve multiple cases, in which case BISECTION will be much faster because of vectorization.

Define LB, UB, target, TolX, and TolFun for each specific application using great care for the following reasons:
- There is no iteration limit, so given an unsolvable task, such as TolX = TolFun = 0, BISECTION remains in an unending loop.
- There is no initial check to make sure that f(x) - target changes sign between LB and UB.
- Very large or very small numbers can introduce numerical issues.

Example 1: Find cube root of array 'target' without using NTHROOT and compare speed to using FZERO (typical result is FZERO taking ~70 times as long as BISECTION).
options = optimset('TolX', 1e-9);
target = [(-100:.1:100)' (-1000:1:1000)'];

tic;
xfz = zeros(size(target));
for ii = 1:numel(target)
xfz(ii) = fzero(@(x) x.^3-target(ii), [-20 20], options);
end
fzero_time = toc

tic;
xbis = bisection(@(x) x.^3, -20, 20, target, options);
bisection_time = toc

fprintf('FZERO took %0.0f times longer than BISECTION.\n', fzero_time/bisection_time)

Example 2: Find roots by varying the function coefficients.
[A, B] = meshgrid(linspace(1,2,6), linspace(4,12,10));
f = @(x) A.*x.^0.2 + B.*x.^0.87 - 15;
xstar = bisection(f,0,5)

Acknowledgements

This file inspired Fuel Fraction Sizing.

MATLAB release MATLAB 8.5 (R2015a)
MATLAB Search Path
`/`
20 Aug 2015 Herbert

06 May 2015 Ole

### Ole (view profile)

If there are several root in the interval does it find the first closes to LB only ?

Comment only
27 Mar 2015 Chi-Fu

### Chi-Fu (view profile)

great program, saved me a lot of time, thank you!

16 Dec 2014 SiddharthK

### SiddharthK (view profile)

The vectorization feature is really really helpful. I was vexed with having to put fzero into for loops.

25 Aug 2014 Philip Ohnewein

### Philip Ohnewein (view profile)

Works brilliantly in my case. Replaces a loop with ~1 million iterations, brings down execution time by several orders of magnitude.
Plus it is well-written and well-documented and a numerically robust method.

03 May 2014 Fabian

### Fabian (view profile)

Excellent file. Much faster than using fzero in a long loop!

10 Mar 2014 Siegmar W

### Siegmar W (view profile)

Thanks! Had a similar file of my own, but yours is better!

25 Jul 2013 Umberto Picchini

### Umberto Picchini (view profile)

I am so glad I found this submission and I'm very grateful to the author for providing an excellent, well-documented code. I had my custom Newton-Raphson algorithm (with provided analytical gradient) invoked thousands of times inside a for loop. I substituted the loop with a single invocation to bisection.m and achieved a 15x acceleration! Awesome.

01 Dec 2012 Sky Sartorius

### Sky Sartorius (view profile)

I just uploaded an entirely new function with almost all new code and documentation and a lot of added features. With so much new code, please let me know if you find a bug.

This is about as far as I'll take this function. I would love to see MathWorks or someone in the community develop a vectorized implementation of Brent's method, i.e. make FZERO vectorized to be able handle array problems. A vectorized FZERO (with a TolFun feature) would be superior to this in every way.

Comment only
29 Nov 2012 Yoshiaki

### Yoshiaki (view profile)

There seems to be a typo on line 80:
jnk = f(UB+LB)/2; % test if f returns multiple outputs

It should be
f((UB+LB)/2)

Comment only
04 Nov 2012 Matthew

### Matthew (view profile)

Excellent documentation with example. Simple function that works as advertised.

24 Jul 2010 1.1

Vectorized

24 Jul 2010 1.2

Vectorized; fixed bug for decreasing functions; some check

02 Nov 2012 1.5

better example, fixed some help typos, tested in 2012b

12 Nov 2012 1.6

made it possible to handle a function that returns array results for scalar input; changed help a bit and added an example for new awesome feature

19 Nov 2012 1.7

fixed a bug that was very rarely throwing out some valid results

27 Nov 2012 1.8

fixed bug that sometimes caused premature convergence; redid funcntion halding to get rid of one of the input checks and simplify the code and make it more understandable; changed some of the help documentation

03 Dec 2012 1.9

Practically all new code and documentation with added features. A good deal of testing done, but with so much new code, please let me know if you find errors.

26 Dec 2012 1.11

New tagline

23 Aug 2013 1.12

minor speed and formatting improvements

20 Aug 2015 1.13

Small changes to help.