• 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. %
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.

Top