SHARE
TWEET

Untitled

a guest Jul 17th, 2019 80 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. function [x, fx, exitFlag] = bisection(f,lb,ub,target,options)
  2. % BISECTION Fast and robust root-finding method that handles n-dim arrays.
  3. %
  4. %   [x,fVal,ExitFlag] = BISECTION(f,LB,UB,target,options) finds x +/- TolX
  5. %   (LB < x < UB) such that f(x) = target +/- TolFun.
  6. %
  7. %   x = BISECTION(f,LB,UB) finds the root(s) of function f on the interval
  8. %   [LB, UB], i.e. finds x such that f(x) = 0 where LB <= x <= UB. f will
  9. %   never be evaluated outside of the interval specified by LB and UB. f
  10. %   should have only one root and f(UB) and f(LB) must bound it. Elements
  11. %   of x are NaN for instances where a solution could not be found.
  12. %
  13. %   x = BISECTION(f,LB,UB,target) finds x such that f(x) = target.
  14. %
  15. %   x = BISECTION(f,LB,UB,target,TolX) will terminate the search when the
  16. %   search interval is smaller than TolX (TolX must be positive).
  17. %
  18. %   x = BISECTION(f,LB,UB,target,options) solves with the default
  19. %   parameters replaced by values in the structure OPTIONS, an argument
  20. %   created with the OPTIMSET function. Used options are TolX and TolFun.
  21. %   Note that OPTIMSET will not allow arrays for tolerances, so set the
  22. %   fields of the options structure manually for non-scalar TolX or TolFun.
  23. %
  24. %   [x,fVal] = BISECTION(f,...) returns the value of f evaluated at x.
  25. %
  26. %   [x,fVal,ExitFlag] = BISECTION(...) returns an ExitFlag that describes
  27. %   the exit condition of BISECTION. Possible values of elements of
  28. %   ExitFlag and the corresponding exit conditions are
  29. %
  30. %       1   Search interval smaller than TolX.
  31. %       2   Function value within TolFun of target.
  32. %       3   Search interval smaller than TolX AND function value within
  33. %           TolFun of target.
  34. %      -1   No solution found.
  35. %
  36. %   Any or all of f(scalar), f(array), LB, UB, target, TolX, or TolFun may
  37. %   be scalar or n-dim arrays. All non-scalar arrays must be the same size.
  38. %   All outputs will be this size.
  39. %
  40. %   Default values are target = 0, TolX = 1e-6, and TolFun = 0.
  41. %
  42. %   There is no iteration limit. This is because BISECTION (with a TolX
  43. %   that won't introduce numerical issues) is guaranteed to converge if f
  44. %   is a continuous function on the interval [UB, LB] and f(x)-target
  45. %   changes sign on the interval.
  46. %
  47. %   The <a href="http://en.wikipedia.org/wiki/Bisection_method">bisection method</a> is a very robust root-finding method. The absolute
  48. %   error is halved at each step so the method converges linearly. However,
  49. %   <a href="http://en.wikipedia.org/wiki/Brent%27s_method">Brent's method</a> (such as implemented in FZERO) can converge
  50. %   superlinearly and is as robust. FZERO also has more features and input
  51. %   checking, so use BISECTION in cases where FZERO would have to be
  52. %   implemented in a loop to solve multiple cases, in which case BISECTION
  53. %   will be much faster because of vectorization.
  54. %
  55. %   Define LB, UB, target, TolX, and TolFun for each specific application
  56. %   using great care for the following reasons:
  57. %     - There is no iteration limit, so given an unsolvable task, BISECTION
  58. %       may remain in an unending loop.
  59. %     - There is no initial check to make sure that f(x) - target changes
  60. %       sign between LB and UB.
  61. %     - Very large or very small numbers can introduce numerical issues.
  62. %
  63. %   Example 1: Find cube root of array 'target' without using NTHROOT and
  64. %   compare speed to using FZERO.
  65. %       options = optimset('TolX', 1e-9);
  66. %       target = [(-100:.1:100)' (-1000:1:1000)'];
  67. %
  68. %       tic;
  69. %       xfz = zeros(size(target));
  70. %       for ii = 1:numel(target)
  71. %           xfz(ii) = fzero(@(x) x.^3-target(ii), [-20 20], options);
  72. %       end
  73. %       fzero_time = toc
  74. %
  75. %       tic;
  76. %       xbis = bisection(@(x) x.^3, -20, 20, target, options);
  77. %       bisection_time = toc
  78. %
  79. %       fprintf('FZERO took %0.0f times longer than BISECTION.\n',...
  80. %                   fzero_time/bisection_time)
  81. %
  82. %   Example 2: Find roots by varying the function coefficients.
  83. %       [A, B] = meshgrid(linspace(1,2,6), linspace(4,12,10));
  84. %       f = @(x) A.*x.^0.2 + B.*x.^0.87 - 15;
  85. %       xstar = bisection(f,0,5)
  86. %
  87. %   See also FZERO, FMINBND, OPTIMSET, FUNCTION_HANDLE.
  88. %
  89. %   [x,fVal,ExitFlag] = BISECTION(f,LB,UB,target,options)
  90. %   FEX URL: http://www.mathworks.com/matlabcentral/fileexchange/28150
  91. %   Copyright 2010-2015 Sky Sartorius
  92. %   Author  - Sky Sartorius
  93. %   Contact - www.mathworks.com/matlabcentral/fileexchange/authors/101715
  94. % --- Process inputs. ---
  95. % Set default values
  96. tolX    = 1e-6;
  97. tolFun  = 0;
  98. if nargin == 5
  99.     if isstruct(options)
  100.         if isfield(options,'TolX') && ~isempty(options.TolX)
  101.             tolX = options.TolX;
  102.         end
  103.         if isfield(options,'TolFun') && ~isempty(options.TolFun)
  104.             tolFun = options.TolFun;
  105.         end
  106.     else
  107.         tolX = options;
  108.     end      
  109. end
  110. if nargin<4 || isempty(target); target=0; end
  111. ub_in = ub; lb_in = lb;
  112. f = @(x) f(x) - target;
  113. % --- Flip UB and LB if necessary. ---
  114. isFlipped = lb>ub;
  115. if any(isFlipped(:))
  116.     ub(isFlipped) = lb_in(isFlipped);
  117.     lb(isFlipped) = ub_in(isFlipped);
  118.     ub_in = ub; lb_in = lb;
  119. end
  120. % --- Make sure everything is the same size for a non-scalar problem. ---
  121. if isscalar(lb) && isscalar(ub)
  122.     % Test if f returns multiple outputs for scalar input.
  123.     if ~isscalar(target)
  124.         ub = ub + zeros(size(target));
  125.     else
  126.         jnk = f(ub);
  127.         if ~isscalar(jnk)
  128.             ub = ub + zeros(size(jnk));
  129.         end
  130.     end
  131. end
  132. % Check if lb and/or ub need to be made into arrays.
  133. if isscalar(lb) && ~isscalar(ub)    
  134.     lb = lb + zeros(size(ub));
  135. elseif ~isscalar(lb) && isscalar(ub)    
  136.     ub = ub + zeros(size(lb));
  137. end
  138. % In newer versions of Matlab, variables should be initialized in the parent
  139. % function.
  140. stillNotDone = [];
  141. outsideTolX = [];
  142. outsideTolFun = [];
  143. testconvergence();
  144. % --- Iterate ---
  145. while any(stillNotDone(:))
  146.     bigger  = fx.*f(ub) > 0;
  147.     ub(bigger)= x(bigger);
  148.     lb(~bigger)= x(~bigger);
  149.    
  150.     testconvergence();
  151. end
  152.     function testconvergence()
  153.         x=(ub+lb)/2;
  154.         fx=f(x);
  155.         outsideTolFun =  abs(fx)  > tolFun;
  156.         outsideTolX =   (ub - lb) > tolX;
  157.         stillNotDone = outsideTolX & outsideTolFun;
  158.     end
  159. % --- Check that f(x+tolX) and f(x-tolX) have opposite sign. ---
  160. fu = f(min(x+tolX,ub_in));
  161. fl = f(max(x-tolX,lb_in));
  162. unboundedRoot = (fu.*fl) > 0;
  163. % Throw out unbounded results if not meeting TolFun convergence criteria.
  164. x(unboundedRoot & outsideTolFun) = NaN;
  165. % --- Catch NaN elements of UB, LB, target, or other funky stuff. ---
  166. x(isnan(fx)) = NaN;
  167. % --- Characterize results. ---
  168. fx = fx + target;
  169. if nargout > 2
  170.     exitFlag                                    = +~outsideTolX;
  171.     exitFlag(~outsideTolFun)                    =  2;
  172.     exitFlag(~outsideTolFun & ~outsideTolX)     =  3;
  173.     exitFlag(isnan(x))                          = -1;
  174. end
  175. end
  176. % V2: July     2010
  177. % V3: December 2012
  178. % don't remember when
  179. %   typo line 39; added fn handle to see also; made array in example 2
  180. %   smaller; changed wording in example 1
  181. % 2013-08-23
  182. %  -changed scalar*ones(...) calls to scalar+zeros(...) calls based on
  183. %   http://undocumentedmatlab.com/blog/allocation-performance-take-2/
  184. %  -rearranged help block and formatted a tiny bit
  185. %  -uploaded to FEX
  186. % 2015-08-20 - updated help; uploaded to FEX
  187. % 2019-02-19 - initialization of variables in parent function
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top