5 - 1/25/2010 Incremental Search (incsearch.m) function [...

Info iconThis preview shows page 1. Sign up to view the full content.

View Full Document Right Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: 1/25/2010 Incremental Search (incsearch.m) function [ brackets ] = incsearch ( f, min, max, points ) % INCSEARCH: locates roots by incremental search % Inputs: f = a function of one variable (need not be able to deal with vectors) % min = lower bound of range to be searched % max = upper bound of range to be searched % points = number of search steps % Outputs: brackets(i, 1) = lower bound for ith b k %O b k (i 1) l b d f i h bracket % brackets(i, 2) = upper bound for ith bracket % **** if no brackets found, brackets = **** nb = 0; brackets = ; % brackets is initially 0 by 0 x = linspace (min, max, points); flo = f(x(1)); for i = 2: points fhi = f(x(i)); = f(x(i)); if sign(flo) ~= sign(fhi) nb = nb + 1; brackets(nb, 1) = x(i 1); brackets(nb, 2) = x(i); end flo = fhi; end end The script below finds and outputs all roots of f between 0 and 20. f = @(x) 50 * sin(0.5 * x) x.^2 17 * x + 60; % from last problem brackets = incsearch ( f, 0, 20, 100 ); [n m] = size(brackets); % note vector on left hand side if n == 0 fprintf ('No roots found.\n'); else for k = 1 : n x = fzero(f, brackets(k, :)); % select whole row fprintf ('There is a root at x = %f\n', x); end end Output: There is a root at x = 1.583780 There is a root at x = 6.636039 There is a root at x = 12.825489 There is a root at x = 16.266248 1 1/25/2010 Bisection Search (basic idea) start with xLOW (less than root) and xHIGH (greater than root) while true pick x pick xROOT half way between xLOW and xHIGH half way between x and x if termination conditions satisfied stop endif if f(xMID) has same sign as f(xLOW) xLOW = xROOT else xHIGH = xROOT endif endwhile Possible Termination Conditions: 1/. The error in xROOT has become acceptably small (i.e. we have got close enough to the actual root). 2/. The function at xROOT is zero or acceptably close to zero. 3/. Some maximum number of iterations has been reached. Condition (1) is the normal termination condition. Condition (2) makes senses in situations in which accuracy in x is unimportant as long as f(x) is close enough to zero. Condition (3) is really more applicable to iterative processes that can go Condition (3) is really more applicable to iterative processes that can go wrong. As the bisection process is essentially guaranteed to succeed there is no real need to guard against disaster. 2 1/25/2010 Max possible error = EMAX = (xHIGH xLOW) / 2 Reduces by a factor of two on every iteration. Let x0 be the original interval width. Max error after n it ti M ft iterations = x 0 2n Number of iterations required to reduce EMAX to EDESIRED = x0 x0 x0 / ln 2 / log 2 ln log log 2 E E E DESIRED DESIRED DESIRED Bisection search unusual in that the number of iterations required to achieve a desired accuracy can be predicted in advance. desired accuracy can be predicted in advance Also unusual in that xROOT can be guaranteed to be within some amount of the correct answer. Text Notes (1) Text uses "approximate error" = EA = |xROOTNEW xROOTOLD| Small changes in xROOT are assumed to mean that error is correspondingly small For a bisection search EA and EMAX are in fact exactly the same thing. y g |xROOTNEW xROOTOLD| is always equal to (xHIGH xLOW) / 2 OLD Assume xROOT xLOW NEW OLD xROOT xROOT xLOW Then xHIGH xLOW x xHIGH x xLOW LOW HIGH 2 2 2 OLD Assume xROOT xHIGH NEW OLD xROOT xROOT xHIGH Then xHIGH xLOW x xLOW HIGH 2 2 3 1/25/2010 Text Notes (2): Text also works with "relative" errors. Relative error = absolute error / magnitude of xROOT The idea is that a given error is less significant when dealing with larger numbers. An error of 0.01 is insignificant if xROOT is 100000 but matters if xROOT is 1. Relative error undefined when xROOT is zero, problematic when xROOT is small. Bisection Search (bisect.m) function [xr] = bisect (f, xl, xh, Edes, display) % BISECT Finds a root by performing a bisection search. % Inputs: f = a function of one variable % xl = lower bound of range to be searched % xh = upper bound of range to be searched % **** f(a) and f(b) must have different signs **** % f(a) and f(b) must have different signs % Edes = tolerance in x (function stops when xr is guaranteed to % be within Edes of a root) % display = display option (0 = no output, 1 = output, defaults to 0) % Outputs: xr = estimate of root if nargin < 5; display = 0; end y ( ); y yl = f(xl); yh = f(xh); ( ); if sign(yl) == sign(yh), error ('no sign change'), end if display fprintf (' step xl xh end xr yr Emax\n'); signOfYl = sign(yl); % keep track of sign of function at xl 4 1/25/2010 for k = 1:1000 % 1000 max iterations xr = (xl + xh) /2; yr = f(xr); Emax = (xh xl) / 2; if display fprintf ('%5d %12.6f %12.6f %12.6f %12.6f %12.6f\n', k, xl, xh, xr, yr, Emax); end if Emax <= Edes || yr == 0 % error acceptably small or direct hit return; end if sign(yr) == signOfYl xl = xr; else xh = xr; xr; end end error ('Bisection search has not converged'); % most unlikely end Bisection with a Casio Calculator Example: root of f(x) = x35 between 1 and 2 Store low limit in memory "A" ( 1 SHIFT STO A ) Store high limit in memory "B" ( 2 SHIFT STO B ) Enter function using M as x ( ALPHA M SHIFT x2 5 ) Evaluate at low limit ( CALC 1 = ), remember sign at low limit Set memory M = (A + B) /2 ( (ALPHA A + ALPHA B) / 2 ) SHIFT STO M ) while true Use up arrow to recall function Evaluate function at M (CALC = ) if sign of result same as sign at low limit Store M in A ( ALPHA M SHIFT STO A ) else Store M in B ( ALPHA M SHIFT STO B ) Store M in B ( ALPHA M SHIFT STO B ) endif Use up arrow to recall (A+B)/2 > M Update M ( = ) endwhile 5 1/25/2010 varargin: Bisection code in text uses vargin. function [root, ea, iter] = bisect (func, xl, xh, es, maxit, varargin) Allows for functions that require parameters in addition to the independent All f f ti th t i t i dditi t th i d d t variable. f = @(x, p1, p1) ....... Parameter values are tacked on to end of argument list when using bisect. root = bisect (f, 0, 10, 0.0001, 1000, 15, 2); % 15 and 2 are values for p1 and p2 bisect passes additional values on to the function when function called test = func (xl, varargin{:}) * func (xr, varargin{:}); For details see text section 3.5.3 False Position (aka Regula Falsi) Similar to bisection but uses different rule for picking xROOT xROOT xHIGH instead of f ( xHIGH )( xLOW xHIGH ) f ( xLOW ) f ( xHIGH ) ( xLOW xHIGH ) 2 xROOT Advantage: Typically (but not always) faster than bisection Disadvantage: Does not allow maximum possible error to be reduced to some Does not allow maximum possible error to be reduced to some chosen amount. Termination must be based on approximate error (stop when changes in xROOT from one iteration to the next become acceptably small). 6 ...
View Full Document

Ask a homework question - tutors are online