This preview shows page 1. Sign up to view the full content.
Unformatted text preview: Chapter 4 Zeros and Roots This chapter describes several basic methods for computing zeros of functions and then combines three of the basic methods into a fast, reliable algorithm known as `zeroin'. 4.1 Let's compute 2. We will use interval bisection, which is a kind of systematic trial and error. We know that 2 is between 1 and 2. Try x = 1 1 . Because x2 is greater 2 than 2, this x is too big. Try x = 1 1 . Because x2 is less than 2, this x is too small. 4 Continuing in this way, our approximations to 2 are 1 1 3 5 13 27 1 , 1 , 1 , 1 , 1 , 1 , .... 2 4 8 16 32 64 Here is a Matlab program, including a step counter. M = 2 a = 1 b = 2 k = 0; while ba > eps x = (a + b)/2; if x^2 > M b = x else a = x end k = k + 1; end
February 15, 2008 Bisection 1 2 Chapter 4. Zeros and Roots We are sure that 2 is in the initial interval [a,b]. This interval is repeatedly cut in half and always brackets the answer. The entire process requires 52 steps. Here are the first few and the last few values. b = 1.50000000000000 a = 1.25000000000000 a = 1.37500000000000 b = 1.43750000000000 a = 1.40625000000000 b = 1.42187500000000 a = 1.41406250000000 b = 1.41796875000000 b = 1.41601562500000 b = 1.41503906250000 b = 1.41455078125000 ..... b = 1.41421356237311 a = 1.41421356237299 a = 1.41421356237305 a = 1.41421356237308 a = 1.41421356237309 b = 1.41421356237310 b = 1.41421356237310 Using format hex, here are the final values of a and b. a = 3ff6a09e667f3bcc b = 3ff6a09e667f3bcd They agree up to the last bit. We haven't actually computed 2, which is irrational and cannot be represented in floating point. But we have found two successive floatingpoint numbers, one on either side of the theoretical result. We've come as close as we can using floatingpoint arithmetic. The process takes 52 steps because there are 52 bits in the fraction of an IEEE doubleprecision number. Each step decreases the interval length by about one bit. Interval bisection is a slow but sure algorithm for finding a zero of f (x), a realvalued function of a real variable. All we assume about the function f (x) is that we can write a Matlab program that evaluates it for any x. We also assume that we know an interval [a, b] on which f (x) changes sign. If f (x) is actually a continuous mathematical function, then there must be a point x somewhere in the interval where f (x ) = 0. But the notion of continuity does not strictly apply to floatingpoint computation. We might not be able to actually find a point where f (x) is exactly zero. Our goal is as follows: Find a very small interval, perhaps two successive floatingpoint numbers, on which the function changes sign. The Matlab code for bisection is 4.2. Newton's Method k = 0; while abs(ba) > eps*abs(b) x = (a + b)/2; if sign(f(x)) == sign(f(b)) b = x; else a = x; end k = k + 1; end 3 Bisection is slow. With the termination condition in the above code, it always takes 52 steps for any function. But it is completely reliable. If we can find a starting interval with a change of sign, then bisection cannot fail to reduce that interval to two successive floatingpoint numbers that bracket the desired result. 4.2 Newton's Method Newton's method for solving f (x) = 0 draws the tangent to the graph of f (x) at any point and determines where the tangent intersects the xaxis. The method requires one starting value, x0 . The iteration is xn+1 = xn  The Matlab code is k = 0; while abs(x  xprev) > eps*abs(x) xprev = x; x = x  f(x)/fprime(x) k = k + 1; end As a method for computing square roots, Newton's method is particularly elegant and effective. To compute M , find a zero of f (x) = x2  M. In this case, f (x) = 2x and x2  M xn+1 = xn  n 2xn 1 M = xn + . 2 xn The algorithm repeatedly averages x and M/x. The Matlab code is f (xn ) . f (xn ) 4 while abs(x  xprev) > eps*abs(x) xprev = x; x = 0.5*(x + M/x); end Here are the results for 2, starting at x = 1. 1.50000000000000 1.41666666666667 1.41421568627451 1.41421356237469 1.41421356237309 1.41421356237309 Chapter 4. Zeros and Roots Newton's method takes only six iterations. In fact, it was done in five, but the sixth iteration was needed to meet the termination condition. When Newton's method works as it does for square roots, it is very effective. It is the basis for many powerful numerical methods. But, as a generalpurpose algorithm for finding zeros of functions, it has three serious drawbacks. The function f (x) must be smooth. It might not be convenient to compute the derivative f (x). The starting guess must be close to the final result. In principle, the computation of the derivative f (x) can be done using a technique known as automatic differentiation. A Matlab function f(x) or a suitable code in any other programming language, defines a mathematical function of its arguments. By combining modern computer science parsing techniques with the rules of calculus, especially the chain rule, it is theoretically possible to generate the code for another function fprime(x), that computes f (x). However, the actual implementation of such techniques is quite complicated and has not yet been fully realized. The local convergence properties of Newton's method are very attractive. Let x be a zero of f (x) and let en = xn  x be the error in the nth iterate. Assume f (x) and f (x) exist and are continuous, x0 is close to x . Then it is possible to prove [2] that en+1 = 1 f () 2 e , 2 f (xn ) n where is some point between xn and x . In other words, en+1 = O(e2 ). n 4.3. A Perverse Example
1.5 5 1 0.5 0 0.5 1 1.5 0 0.5 1 1.5 2 2.5 3 3.5 4 Figure 4.1. Newton's method in an infinite loop. This is called quadratic convergence. For nice, smooth functions, once you are close enough to the zero, the error is roughly squared with each iteration. The number of correct digits approximately doubles with each iteration. The results we saw for 2 are typical. When the assumptions underlying the local convergence theory are not satisfied, Newton's method might be unreliable. If f (x) does not have continuous, bounded first and second derivatives, or if the starting point is not close enough to the zero, then the local theory does not apply and we might get slow convergence, or even no convergence at all. The next section provides one example of what might happen. 4.3 A Perverse Example
xn+1 = xn  f (xn ) f (xn ) Let's see if we can get Newton's method to iterate forever. The iteration cycles back and forth around a point a if xn+1  a = (xn  a). This happens if f (x) satisfies xa f (x) = (x  a). f (x) This is a separable ordinary differential equation: f (x) 1 = . f (x) 2(x  a) 6 The solution is Chapter 4. Zeros and Roots f (x) = sign(x  a) x  a. The zero of f (x) is, of course, at x = a. A plot of f (x), Figure 4.1, with a = 2, is obtained with ezplot('sign(x2)*sqrt(abs(x2))',0,4) If we draw the tangent to the graph at any point, it intersects the xaxis on the opposite side of x = a. Newton's method cycles forever, neither converging nor diverging. The convergence theory for Newton's method fails in this case because f (x) is unbounded as x a. It is also interesting to apply the algorithms discussed in the next sections to this function. 4.4 Secant Method The secant method replaces the derivative evaluation in Newton's method with a finite difference approximation based on the two most recent iterates. Instead of drawing a tangent to the graph of f (x) at one point, you draw a secant through two points. The next iterate is the intersection of this secant with the xaxis. The iteration requires two starting values, x0 and x1 . The subsequent iterates are given by sn = xn+1 f (xn )  f (xn1 ) , xn  xn1 f (xn ) = xn  . sn This formulation makes it clear how Newton's f (xn ) is being replaced with the slope of the secant, sn . The formulation in the following Matlab code is a little more compact: while abs(ba) > eps*abs(b) c = a; a = b; b = b + (b  c)/(f(c)/f(b)1); k = k + 1; end For 2, starting with a = 1 and b = 2, the secant method requires seven iterations, compared with Newton's six. 1.33333333333333 1.40000000000000 1.41463414634146 1.41421143847487 1.41421356205732 1.41421356237310 1.41421356237310 4.5. Inverse Quadratic Interpolation 7 The secant method's primary advantage over Newton's method is that it does not require code to compute f (x). Its convergence properties are similar. Again, assuming f (x) and f (x) are continuous, it is possible to prove [2] that en+1 = 1 f ()f (n )f (n1 ) en en1 , 2 f ()3 where is some point between xn and x . In other words, en+1 = O(en en1 ). This is not quadratic convergence, but it is superlinear convergence. It turns out that en+1 = O(e ), n where is the golden ratio, (1 + 5)/2. Once you get close, the number of correct digits is roughly multiplied by 1.6 with each iteration. That's almost as fast as Newton's method and a whole lot faster than the one bit per step produced by bisection. We leave it to exercise 4.8 to investigate the behavior of the secant method on the perverse function from the previous section: f (x) = sign(x  a) x  a. 4.5 Inverse Quadratic Interpolation The secant method uses two previous points to get the next one, so why not use three? Suppose we have three values, a, b, and c, and corresponding function values, f (a), f (b), and f (c). We could interpolate these values by a parabola, a quadratic function of x, and take the next iterate to be the point where the parabola intersects the xaxis. The difficulty is that the parabola might not intersect the xaxis; a quadratic function does not necessarily have real roots. This could be regarded as an advantage. An algorithm known as Muller's method uses the complex roots of the quadratic to produce approximations to complex zeros of f (x). But, for now, we want to avoid complex arithmetic. Instead of a quadratic in x, we can interpolate the three points with a quadratic function in y. That's a "sideways" parabola, P (y), determined by the interpolation conditions a = P (f (a)), b = P (f (b)), c = P (f (c)). This parabola always intersects the xaxis, which is y = 0. So x = P (0) is the next iterate. This method is known as inverse quadratic interpolation. We will abbreviate it with IQI. Here is Matlab code that illustrates the idea. k = 0; while abs(cb) > eps*abs(c) 8 x a b c k end = = = = = Chapter 4. Zeros and Roots polyinterp([f(a),f(b),f(c)],[a,b,c],0) b; c; x; k + 1; The trouble with this "pure" IQI algorithm is that polynomial interpolation requires the abscissae, which in this case are f (a), f (b), and f (c), to be distinct. We have no guarantee that they are. For example, if we try to compute 2 using f (x) = x2  2 and start with a = 2, b = 0, c = 2, we are starting with f (a) = f (c) and the first step is undefined. If we start near this singular situation, say with a = 2.001, b = 0, c = 1.999, the next iterate is near x = 500. So IQI is like an immature racehorse. It moves very quickly when it is near the finish line, but its global behavior can be erratic. It needs a good trainer to keep it under control. 4.6 Zeroin The idea behind the zeroin algorithm is to combine the reliability of bisection with the convergence speed of secant and IQI. T. J. Dekker and colleagues at the Mathematical Center in Amsterdam developed the first version of the algorithm in the 1960s [3]. Our implementation is based on a version by Richard Brent [1]. Here is the outline: Start with a and b so that f (a) and f (b) have opposite signs. Use a secant step to give c between a and b. Repeat the following steps until b  a < b or f (b) = 0. Arrange a, b, and c so that f (a) and f (b) have opposite signs, f (b) f (a), c is the previous value of b. If c = a, consider an IQI step. If c = a, consider a secant step. If the IQI or secant step is in the interval [a, b], take it. If the step is not in the interval, use bisection. This algorithm is foolproof. It never loses track of the zero trapped in a shrinking interval. It uses rapidly convergent methods when they are reliable. It uses a slow, but sure, method when it is necessary. 4.7. fzerotx
J (x)
0 9 1 0.5 0 0.5 0 5 10 15 x 20 25 30 Figure 4.2. Zeros of J0 (x). 4.7 fzerotx The Matlab implementation of the zeroin algorithm is called fzero. It has several features beyond the basic algorithm. A preamble takes a single starting guess and searches for an interval with a sign change. The values returned by the function f(x) are checked for infinities, NaNs, and complex numbers. Default tolerances can be changed. Additional output, including a count of function evaluations, can be requested. Our textbook version of zeroin is fzerotx. We have simplified fzero by removing most of its additional features while retaining the essential features of zeroin. We can illustrate the use of fzerotx with the zerothorder Bessel function of the first kind, J0 (x). This function is available in Matlab as besselj(0,x). The following code finds the first 10 zeros of J0 (x) and produces Figure 4.2 (except for the red `x', which we will add later). J0 = @(x) besselj(0,x); for n = 1:10 z(n) = fzerotx(J0,[(n1) n]*pi); end x = 0:pi/50:10*pi; y = J0(x); plot(z,zeros(1,10),'o',x,y,'') line([0 10*pi],[0 0],'color','black') axis([0 10*pi 0.5 1.0]) You can see from the figure that the graph of J0 (x) is like an amplitude and frequency modulated version of cos x. The distance between successive zeros is close to . 10 Chapter 4. Zeros and Roots The function fzerotx takes two arguments. The first specifies the function F (x) whose zero is being sought and the second specifies the interval [a, b] to search. fzerotx is an example of a Matlab function function, which is a function that takes another function as an argument. ezplot is another example. Other chapters of this bookChapter 6, Quadrature; Chapter 7, Ordinary Differential Equations; and even Chapter 9, Random Numbersalso describe "tx" and "gui" Mfiles that are function functions. A function can be passed as an argument to another function in two different ways: function handle, A function handle uses the '@' character preceding the name of a builtin function or a function defined in an Mfile. Examples include @cos @humps @bessj0 where bessj0.m is the twoline Mfile function y = bessj0(x) y = besselj(0,x) These handles can then be used as arguments to function functions. z = fzerotx(@bessj0,[0,pi]) Note that @besselj is also a valid function handle, but for a function of two arguments. Anonymous functions were introduced in Matlab 7. Examples include F = @(t) cos(pi*t) G = @(z) z^32*z5 J0 = @(x) besselj(0,x) These objects are called anonymous functions because the construction @(arguments) expression defines a function, but does not give it a name. Mfiles and anonymous functions can define functions of more than one argument. In this case, the values of the extra arguments can be passed through fzerotx to the objective function. These values remain constant during the zero finding iteration. This allows us to find where a particular function takes on a specified value y, instead of just finding a zero. For example, consider the equation J0 () = 0.5. Define an anonymous function with two or three arguments. anonymous function. 4.7. fzerotx F = @(x,y) besselj(0,x)y or B = @(x,n,y) besselj(n,x)y Then either xi = fzerotx(F,[0,2],.5) or xi = fzerotx(B,[0,2],0,.5) produces xi = 1.5211 The point (, J0 ()) is marked with an 'x' in Figure 4.2. The preamble for fzerotx is as follows. function b = fzerotx(F,ab,varargin); %FZEROTX Textbook version of FZERO. % x = fzerotx(F,[a,b]) tries to find a zero of F(x) between % a and b. F(a) and F(b) must have opposite signs. % fzerotx returns one endpoint of a small subinterval of % [a,b] where F changes sign. % Additional arguments, fzerotx(F,[a,b],p1,p2,...), % are passed on, F(x,p1,p2,...). 11 The first section of code initializes the variables a, b, and c that characterize the search interval. The function F is evaluated at the endpoints of the initial interval. a = ab(1); b = ab(2); fa = F(a,varargin{:}); fb = F(b,varargin{:}); if sign(fa) == sign(fb) error('Function must change sign on the interval') end c = a; fc = fa; d = b  c; e = d; Here is the beginning of the main loop. At the start of each pass through the loop, a, b, and c are rearranged to satisfy the conditions of the zeroin algorithm. 12 while fb ~= 0 Chapter 4. Zeros and Roots % The three current points, a, b, and c, satisfy: % f(x) changes sign between a and b. % abs(f(b)) <= abs(f(a)). % c = previous b, so c might = a. % The next point is chosen from % Bisection point, (a+b)/2. % Secant point determined by b and c. % Inverse quadratic interpolation point determined % by a, b, and c if they are distinct. if sign(fa) == sign(fb) a = c; fa = fc; d = b  c; e = d; end if abs(fa) < abs(fb) c = b; b = a; a = c; fc = fb; fb = fa; fa = fc; end This section is the convergence test and possible exit from the loop. m = 0.5*(a  b); tol = 2.0*eps*max(abs(b),1.0); if (abs(m) <= tol)  (fb == 0.0), break end The next section of code makes the choice between bisection and the two flavors of interpolation. % Choose bisection or interpolation if (abs(e) < tol)  (abs(fc) <= abs(fb)) % Bisection d = m; e = m; else % Interpolation s = fb/fc; if (a == c) % Linear interpolation (secant) p = 2.0*m*s; q = 1.0  s; else % Inverse quadratic interpolation q = fc/fa; r = fb/fa; 4.8. fzerogui p = s*(2.0*m*q*(q  r)  (b  c)*(r  1.0)); q = (q  1.0)*(r  1.0)*(s  1.0); end; if p > 0, q = q; else p = p; end; % Is interpolated point acceptable if (2.0*p < 3.0*m*q  abs(tol*q)) & (p < abs(0.5*e*q)) e = d; d = p/q; else d = m; e = m; end; end The final section evaluates F at the next iterate. % Next point c = b; fc = fb; if abs(d) > tol b = b + d; else b = b  sign(ba)*tol; end fb = F(b,varargin{:}); end 13 4.8 fzerogui The Mfile fzerogui demonstrates the behavior of zeroin and fzerotx. At each step of the iteration, you are offered a chance to choose the next point. The choice always includes the bisection point, shown in red on the computer screen. When there are three distinct points active, a, b, and c, the IQI point is shown in blue. When a = c, so there are only two distinct points, the secant point is shown in green. A plot of f (x) is also provided as a dotted line, but the algorithm does not "know" these other function values. You can choose any point you like as the next iterate. You do not have to follow the zeroin algorithm and choose the bisection or interpolant point. You can even cheat by trying to pick the point where the dotted line crosses the axis. We can demonstrate how fzerogui behaves by seeking the first zero of the Bessel function. It turns out that the first local minimum of J0 (x) is located near x = 3.83. So here are the first few steps of fzerogui('besselj(0,x)',[0 3.83]) Initially, c = b, so the two choices are the bisection point and the secant point (Figure 4.3). 14 Chapter 4. Zeros and Roots 1 0.8 0.6 0.4 0.2 a 0 c 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 3 3.5 4 b Figure 4.3. Initially, choose secant or bisection. If you choose the secant point, then b moves there and J0 (x) is evaluated at x = b. We then have three distinct points, so the choice is between the bisection point and the IQI point (Figure 4.4).
1 0.8 0.6 0.4 0.2 a 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 3 3.5 4 b c Figure 4.4. Choose IQI or bisection. If you choose the IQI point, the interval shrinks, the GUI zooms in on the reduced interval, and the choice is again between the bisection and secant points, which now happen to be close together (Figure 4.5). You can choose either point, or any other point close to them. After two more steps, the interval shrinks again and the situation shown in Figure 4.6 is reached. This is the typical configuration as we approach convergence. The graph of the function looks nearly like a straight line and the secant or IQI point is much closer 4.9. Value Finding and Reverse Interpolation 15 0.15 0.1 0.05 b 0 c 0.05 a 0.1 0.15 2.2 2.3 2.4 2.5 2.6 2.7 Figure 4.5. Secant and bisection points nearly coincide. to the desired zero than the bisection point. It now becomes clear that choosing secant or IQI will lead to much faster convergence than bisection. After several more steps, the length of the interval containing a change of sign is reduced to a tiny fraction of the original length, and the algorithm terminates, returning the final b as its result. 0.15 0.1 0.05 a 0 b c 0.05 0.1 0.15 2.15 2.2 2.25 2.3 2.35 2.4 Figure 4.6. Nearing convergence. 4.9 Value Finding and Reverse Interpolation These two problems look very similar. Given a function F (x) and a value , find so that F () = . 16 Chapter 4. Zeros and Roots Given data (xk , yk ) that sample an unknown function F (x), and a value , find so that F () = . For the first problem, we are able to evaluate F (x) at any x, so we can use a zero finder on the translated function f (x) = F (x)  . This gives us the desired so that f () = 0, and hence F () = . For the second problem, we need to do some kind of interpolation. The most obvious approach is to use a zero finder on f (x) = P (x)  , where P (x) is some interpolant, such as pchiptx(xk,yk,x) or splinetx(xk,yk,x). This often works well, but it can be expensive. The zero finder calls for repeated evaluation of the interpolant. With the implementations we have in this book, that involves repeated calculation of the interpolant's parameters and repeated determination of the appropriate interval index. A sometimes preferable alternative, known as reverse interpolation, uses pchip or spline with the roles of xk and yk reversed. This requires monotonicity in the yk , or at least in a subset of the yk around the target value . A different piecewise polynomial, say Q(y), is created with the property that Q(yk ) = xk . Now it is not necessary to use a zero finder. We simply evaluate = Q(y) at y = . The choice between these two alternatives depends on how the data are best approximated by piecewise polynomials. Is it better to use x or y as the independent variable? 4.10 Optimization and fmintx The task of finding maxima and minima of functions is closely related to zero finding. In this section, we describe an algorithm similar to zeroin that finds a local minimum of a function of one variable. The problem specification involves a function f (x) and an interval [a, b]. The objective is to find a value of x that gives a local minimum of f (x) on the given interval. If the function is unimodular, that is, has only one local minimum on the interval, then that minimum will be found. But if there is more than one local minimum, only one of them will be found, and that one will not necessarily be minimum for the entire interval. It is also possible that one of the endpoints is a minimizer. Interval bisection cannot be used. Even if we know the values of f (a), f (b), and f ((a + b)/2), we cannot decide which half of the interval to discard and still keep the minimum enclosed. Interval trisection is feasible, but inefficient. Let h = (b  a)/3, so u = a + h and v = b  h divide the interval into three equal parts. Assume we find that f (u) < f (v). Then we could replace b with v, thereby cutting the interval length by a factor of twothirds, and still be sure that the minimum is in the reduced interval. However, u would be in the center of the new interval and would not be useful in the next step. We would need to evaluate the function twice each step. The natural minimization analogue of bisection is golden section search. The idea is illustrated for a = 0 and b = 1 in Figure 4.7. Let h = (b  a), where is a quantity a little bit larger than 1/3 that we have yet to determine. Then the points u = a + h and v = b  h divide the interval into three unequal parts. For the 4.10. Optimization and fmintx 1 17 0 u v 1 0 v' u 1 Figure 4.7. Golden section search. first step, evaluate both f (u) and f (v). Assume we find that f (u) < f (v). Then we know the minimum is between a and v. We can replace b with v and repeat the process. If we choose the right value for , the point u is in the proper position to be used in the next step. After the first step, the function has to be evaluated only once each step. The defining equation for is 1 = , 1 1 2  3 + 1 = 0. The solution is = 2  = (3  5)/2 0.382. or Here is the golden ratio that we used to introduce Matlab in the first chapter of this book. With golden section search, the length of the interval is reduced by a factor of  1 0.618 each step. It would take 52 75 log2 (  1) steps to reduce the interval length to roughly eps, the size of IEEE doubleprecision roundoff error, times its original value. After the first few steps, there is often enough history to give three distinct points and corresponding function values in the active interval. If the minimum of the parabola interpolating these three points is in the interval, then it, rather than the golden section point, is usually a better choice for the next point. This combination of golden section search and parabolic interpolation provides a reliable and efficient method for onedimensional optimization. The proper stopping criteria for optimization searches can be tricky. At a minimum of f (x), the first derivative f (x) is zero. Consequently, near a minimum, f (x) acts like a quadratic with no linear term: f (x) a + b(x  c)2 + . 18 Chapter 4. Zeros and Roots The minimum occurs at x = c and has the value f (c) = a. If x is close to c, say x c + for small , then f (x) a + b 2 . Small changes in x are squared when computing function values. If a and b are comparable in size, and nonzero, then the stopping criterion should involve sqrt(eps) because any smaller changes in x will not affect f (x). But if a and b have different orders of magnitude, or if either a or c is nearly zero, then interval lengths of size eps, rather than sqrt(eps), are appropriate. Matlab includes a function function fminbnd that uses golden section search and parabolic interpolation to find a local minimum of a realvalued function of a real variable. The function is based upon a Fortran subroutine by Richard Brent [1]. Matlab also includes a function function, fminsearch, that uses a technique known as the NelderMeade simplex algorithm to search for a local minimum of a realvalued function of several real variables. The Matlab Optimization Toolbox is a collection of programs for other kinds of optimization problems, including constrained optimization, linear programming, and largescale, sparse optimization. Our NCM collection includes a function fmintx that is a simplified textbook version of fminbnd. One of the simplifications involves the stopping criterion. The search is terminated when the length of the interval becomes less than a specified parameter tol. The default value of tol is 106 . More sophisticated stopping criteria involving relative and absolute tolerances in both x and f (x) are used in the full codes.
humps(x) 0 20 40 60 80 100 1 0.5 0 0.5 x 1 1.5 2 Figure 4.8. Finding the minimum of humps(x). The Matlab demos directory includes a function named humps that is intended to illustrate the behavior of graphics, quadrature, and zerofinding routines. The function is 1 1 h(x) = + . (x  0.3)2 + 0.01 (x  0.9)2 + 0.04 Exercises The statements F = @(x) humps(x); fmintx(F,1,2,1.e4) 19 take the steps shown in Figure 4.8 and in the following output. We see that golden section search is used for the second, third, and seventh steps, and that parabolic interpolation is used exclusively once the search nears the minimizer. step init: gold: gold: para: para: para: gold: para: para: para: para: para: x 0.1458980337 0.8541019662 0.2917960675 0.4492755129 0.4333426114 0.3033578448 0.2432135488 0.3170404333 0.2985083078 0.3003583547 0.3003763623 0.3003756221 f(x) 25.2748253202 20.9035150009 2.5391843579 29.0885282699 33.8762343193 96.4127439649 71.7375588319 93.8108500149 96.4666018623 96.5014055840 96.5014085540 96.5014085603 Exercises
4.1. Use fzerogui to try to find a zero of each of the following functions in the given interval. Do you see any interesting or unusual behavior? x3  2x  5 sin x x3  0.001 log (x + 2/3) sign(x  2) x  2 atan(x)  /3 1/(x  ) [0, 3] [1, 4] [1, 1] [0, 1] [1, 4] [0, 5] [0, 5] 4.2. Here is a little footnote to the history of numerical methods. The polynomial x3  2x  5 was used by Wallis when he first presented Newton's method to the French Academy. It has one real root, between x = 2 and x = 3, and a pair of complex conjugate roots. (a) Use the Symbolic Toolbox to find symbolic expressions for the three roots. Warning: The results are not pretty. Convert the expressions to numerical values. 20 Chapter 4. Zeros and Roots (b) Use the roots function in Matlab to find numerical values for all three roots. (c) Use fzerotx to find the real root. (d) Use Newton's method starting with a complex initial value to find a complex root. (e) Can bisection be used to find the complex root? Why or why not? 4.3. Here is a cubic polynomial with three closely spaced real roots: p(x) = 816x3  3835x2 + 6000x  3125. (a) What are the exact roots of p? (b) Plot p(x) for 1.43 x 1.71. Show the location of the three roots. (c) Starting with x0 = 1.5, what does Newton's method do? (d) Starting with x0 = 1 and x1 = 2, what does the secant method do? (e) Starting with the interval [1, 2], what does bisection do? (f) What is fzerotx(p,[1,2])? Why? What causes fzerotx to terminate? (a) How does fzerotx choose between the bisection point and the interpolant point for its next iterate? (b) Why is the quantity tol involved in the choice? Derive the formula that fzerotx uses for IQI. Hoping to find the zero of J0 (x) in the interval 0 x , we might try the statement z = fzerotx(@besselj,[0 pi],0) This is legal usage of a function handle, and of fzerotx, but it produces z = 3.1416. Why? 4.8. Investigate the behavior of the secant method on the function f (x) = sign(x  a) x  a. 4.4. 4.5. 4.6. 4.7. 4.9. Find the first ten positive values of x for which x = tan x. 4.10. (a) Compute the first ten zeros of J0 (x). You can use our graph of J0 (x) to estimate their location. (b) Compute the first ten zeros of Y0 (x), the zerothorder Bessel function of the second kind. (c) Compute all the values of x between 0 and 10 for which J0 (x) = Y0 (x). (d) Make a composite plot showing J0 (x) and Y0 (x) for 0 x 10, the first ten zeros of both functions, and the points of intersection. 4.11. The gamma function is defined by an integral: (x + 1) = tx et dt.
0 Integration by parts shows that, when evaluated at the integers, (x) interpolates the factorial function (n + 1) = n!. Exercises 21 (x) and n! grow so rapidly that they generate floatingpoint overflow for relatively small values of x and n. It is often more convenient to work with the logarithms of these functions. The Matlab functions gamma and gammaln compute (x) and log (x), respectively. The quantity n! is easily computed by the expression prod(1:n) but many people expect there to be a function named factorial, so Matlab has such a function. (a) What is the largest value of n for which (n + 1) and n! can be exactly represented by a doubleprecision floatingpoint number? (b) What is the largest value of n for which (n + 1) and n! can be approximately represented by a doubleprecision floatingpoint number that does not overflow? 4.12. Stirling's approximation is a classical estimate for log (x + 1): log (x + 1) x log(x)  x + 1 log(2x). 2 Bill Gosper [4] has noted that a better approximation is log (x + 1) x log(x)  x + 1 log(2x + /3). 2 The accuracy of both approximations improves as x increases. (a) What is the relative error in Stirling's approximation and in Gosper's approximation when x = 2? (b) How large must x be for Stirling's approximation and for Gosper's approximation to have a relative error less than 106 ? 4.13. The statements y = 2:.01:10; x = gammaln(y); plot(x,y) produce a graph of the inverse of the log function. (a) Write a Matlab function gammalninv that evaluates this function for any x. That is, given x, y = gammalninv(x) computes y so that gammaln(y) is equal to x. (b) What are the appropriate ranges of x and y for this function? 4.14. Here is a table of the distance, d, that a hypothetical vehicle requires to stop if the brakes are applied when it is traveling with velocity v. 22 v(m/s) 0 10 20 30 40 50 60 d(m) 0 5 20 46 70 102 153 Chapter 4. Zeros and Roots What is the speed limit for this vehicle if it must be able to stop in at most 60 m? Compute the speed three different ways. (a) piecewise linear interpolation, (b) piecewise cubic interpolation with pchiptx, (c) reverse piecewise cubic interpolation with pchiptx. Because these are wellbehaved data, the three values are close to each other, but not identical. 4.15. Kepler's model of planetary orbits includes a quantity E, the eccentricity anomaly, that satisfies the equation M = E  e sin E, where M is the mean anomaly and e is the eccentricity of the orbit. For this exercise, take M = 24.851090 and e = 0.1. (a) Use fzerotx to solve for E. You can assign the appropriate values to M and e and then use them in the definition of a function of E. M = 24.851090 e = 0.1 F = @(E) E  e*sin(E)  M Use F for the first argument to fzerotx. (b) An "exact" formula for E is known: E =M +2 1 Jm (me) sin(mM ), m m=1 where Jm (x) is the Bessel function of the first kind of order m. Use this formula, and besselj(m,x) in Matlab , to compute E. How many terms are needed? How does this value of E compare to the value obtained with fzerotx? 4.16. Utilities must avoid freezing water mains. If we assume uniform soil conditions, the temperature T (x, t) at a distance x below the surface and time t after the beginning of a cold snap is given approximately by T (x, t)  Ts x = erf . Ti  Ts 2 t Here Ts is the constant surface temperature during the cold period, Ti is the initial soil temperature before the cold snap, and is the thermal conductivity Exercises 23 4.17. 4.18. 4.19. 4.20. 4.21. of the soil. If x is measured in meters and t in seconds, then = 0.138 106 m2 /s. Let Ti = 20 C, and Ts = 15 C, and recall that water freezes at 0 C. Use fzerotx to determine how deep a water main should be buried so that it will not freeze until at least 60 days' exposure under these conditions. Modify fmintx to provide printed and graphical output similar to that at the end of section ??. Reproduce the results shown in Figure 4.8 for humps(x). Let f (x) = 9x2  6x + 2. What is the actual minimizer of f (x)? How close to the actual minimizer can you get with fmintx? Why? Theoretically, fmintx(@cos,2,4,eps) should return pi. How close does it get? Why? On the other hand, fmintx(@cos,0,2*pi) does return pi. Why? If you use tol = 0 with fmintx(@F,a,b,tol), does the iteration run forever? Why or why not? Derive the formulas for minimization by parabolic interpolation used in the following portion of fmintx: r = (x  w)*(fx  fv); q = (x  v)*(fx  fw); p = (x  v)*q  (x  w)*r; s = 2.0*(q  r); if s > 0.0, p = p; end s = abs(s); % Is the parabola acceptable? para = ( (abs(p)<abs(0.5*s*e)) & (p > s*(a  x)) & (p < s*(b  x)) ); if para e = d; d = p/s; newx = x + d; end 4.22. Let f (x) = sin (tan x)  tan (sin x), 0 x . (a) Plot f (x). (b) Why is it difficult to compute the minimum of f (x)? (c) What does fmintx compute as the minimum of f (x)? (d) What is the limit as x /2 of f (x)? (e) What is the glb or infimum of f (x)? 24 Chapter 4. Zeros and Roots Bibliography
[1] R. P. Brent, Algorithms for Minimization Without Derivatives, PrenticeHall, Englewood Cliffs, NJ, 1973. [2] G. Dahlquist and A. Bjorck, Numerical Methods, PrenticeHall, Englewood Cliffs, NJ, 1974. [3] T. J. Dekker, Finding a zero by means of successive linear interpolation, in Constructive Aspects of the Fundamental Theorem of Algebra, B. Dejon and P. Henrici (editors), WileyInterscience, New York, 1969, pp. 3748. [4] E. Weisstein, World of Mathematics, Stirling's Approximation. http://mathworld.wolfram.com/StirlingsApproximation.html 25 ...
View
Full
Document
This note was uploaded on 03/11/2012 for the course MAT 421 taught by Professor Staff during the Fall '11 term at ASU.
 Fall '11
 Staff

Click to edit the document details