• API
• FAQ
• Tools
• Archive
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. %
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