%% Root-finding and minimization %% Introduction % % Finding the roots of functions, and minimum or maximum values is a % fundamental problem in computational mathematics. Even problems that % appear quite simple cannot be solved analytically. For example, consider % solving this scalar equation % $$xe^x = 7$$ % Despite the fact that the left hand side involves only elementary % functions, the solution cannot be written in terms of elementary % functions. Yet we can easily show graphically that this equation has a % solution. We do this by defining % $f(x) = xe^x - 7$ and checking that % $f(x)$ crosses the x-axis. % f = @(x) x.*exp(x) - 7; x = linspace(0,2,100); plot(x,f(x),'k','linewidth',2); hold on; plot(x,0*x,'k'); xlabel('x','fontsize',16); ylabel('f(x)','fontsize',16); title('Graphical solution to root-finding','fontsize',18); snapnow; %% % % To get a better approximation to the root, however, we will need to % explore numerical root finding techniques. For scalar equations, Matlab % has a built in root-finding function called fzero. % %% % % From this graph, we estimate that the solution to our original equation % is roughly $x \approx 1.5$. % %% The fzero function % %
% help fzero
%  fzero  Single-variable nonlinear zero finding.
%     X = fzero(FUN,X0) tries to find a zero of the function FUN near X0,
%     if X0 is a scalar.  It first finds an interval containing X0 where the
%     function values of the interval endpoints differ in sign, then searches
%     that interval for a zero.  FUN is a function handle.  FUN accepts real
%     scalar input X and returns a real scalar function value F, evaluated
%     at X. The value X returned by fzero is near a point where FUN changes
%     sign (if FUN is continuous), or NaN if the search fails.
%
%     X = fzero(FUN,X0), where X0 is a vector of length 2, assumes X0 is a
%     finite interval where the sign of FUN(X0(1)) differs from the sign of
%     FUN(X0(2)). An error occurs if this is not true.  Calling fzero with a
%     finite interval guarantees fzero will return a value near a point where
%     FUN changes sign.
%
%     X = fzero(FUN,X0), where X0 is a scalar value, uses X0 as a starting
%     guess. fzero looks for an interval containing a sign change for FUN and
%     containing X0.  If no such interval is found, NaN is returned.
%     In this case, the search terminates when the search interval
%     is expanded until an Inf, NaN, or complex value is found. Note: if
%     the option FunValCheck is 'on', then an error will occur if an NaN or
%     complex value is found.
%
%     X = fzero(FUN,X0,OPTIONS) solves the equation with the default optimization
%     parameters replaced by values in the structure OPTIONS, an argument
%     created with the OPTIMSET function.  See OPTIMSET for details.  Used
%     options are Display, TolX, FunValCheck, OutputFcn, and PlotFcns.
%
%     X = fzero(PROBLEM) finds the zero of a function defined in PROBLEM.
%     PROBLEM is a structure with the function FUN in PROBLEM.objective,
%     the start point in PROBLEM.x0, the options structure in PROBLEM.options,
%     and solver name 'fzero' in PROBLEM.solver. The structure PROBLEM must have
%     all the fields.
%
%     [X,FVAL]= fzero(FUN,...) returns the value of the function described
%     in FUN, at X.
%
%     [X,FVAL,EXITFLAG] = fzero(...) returns an EXITFLAG that describes the
%     exit condition of fzero. Possible values of EXITFLAG and the corresponding
%     exit conditions are
%
%       1  fzero found a zero X.
%      -1  Algorithm terminated by output function.
%      -3  NaN or Inf function value encountered during search for an interval
%           containing a sign change.
%      -4  Complex function value encountered during search for an interval
%           containing a sign change.
%      -5  fzero may have converged to a singular point.
%      -6  fzero can not detect a change in sign of the function.
%
%     [X,FVAL,EXITFLAG,OUTPUT] = fzero(...) returns a structure OUTPUT
%     with the number of function evaluations in OUTPUT.funcCount, the
%     algorithm name in OUTPUT.algorithm, the number of iterations to
%     find an interval (if needed) in OUTPUT.intervaliterations, the
%     number of zero-finding iterations in OUTPUT.iterations, and the
%     exit message in OUTPUT.message.
%
%     Examples
%       FUN can be specified using @:
%          X = fzero(@sin,3)
%       returns pi.
%          X = fzero(@sin,3,optimset('Display','iter'))
%       returns pi, uses the default tolerance and displays iteration information.
%
%       FUN can be an anonymous function:
%          X = fzero(@(x) sin(3*x),2)
%
%       FUN can be a parameterized function.  Use an anonymous function to
%       capture the problem-dependent parameters:
%          f = @(x,c) cos(c.*x);  % The parameterized function.
%          c = 2;                 % The parameter.
%          X = fzero(@(x) myfun(x,c),0.1)
%
%     Limitations
%          X = fzero(@(x) abs(x)+1, 1)
%       returns NaN since this function does not change sign anywhere on the
%       real axis (and does not have a zero as well).
%          X = fzero(@tan,2)
%       returns X near 1.5708 because the discontinuity of this function near the
%       point X gives the appearance (numerically) that the function changes sign at X.
%
%
%     Reference page in Help browser
%        doc fzero
% 
%
% This function, like all other numerical root finding techniques, finds % the root by making a series of intelligent guesses. The guesses are % based on one of several algorithms, including Bisection, the Secant % method and Inverse Quadratic Interpolation. Because fzero is iterative, we must supply it with a good % initial guess to our solution, or an interval that we believe brackets % the root. Here are two examples of how to call fzero to solve the problem above. % f = @(x) x.*exp(x) - 7; x0 = 1.5; x_soln = fzero(f,x0); fprintf('Solution is : %24.16f\n',x_soln); fprintf('Residual error is : %16.8e\n',f(x_soln)); %% % % Alternatively, we can supply an interval that we believe contains the % root : % x0 = [1.4 1.6]; x_soln = fzero(f,x0); fprintf('Solution is : %24.16f\n',x_soln); fprintf('Residual error is : %16.8e\n',f(x_soln)); %% Function minimization % % To find the minimum values of a function in an interval, we can use the % Matlab command fminbnd. This function has a % very similar structure to the fzero function in % that it takes a function handle. But it will require an interval over % which to search for the minimum value. %
%
% help fminbnd
%  fminbnd Single-variable bounded nonlinear function minimization.
%     X = fminbnd(FUN,x1,x2) attempts to find  a local minimizer X of the function
%     FUN in the interval x1 < X < x2.  FUN is a function handle.  FUN accepts
%     scalar input X and returns a scalar function value F evaluated at X.
%
%     X = fminbnd(FUN,x1,x2,OPTIONS) minimizes with the default optimization
%     parameters replaced by values in the structure OPTIONS, created with
%     the OPTIMSET function. See OPTIMSET for details. fminbnd uses these
%     options: Display, TolX, MaxFunEval, MaxIter, FunValCheck, PlotFcns,
%     and OutputFcn.
%
%     X = fminbnd(PROBLEM) finds the minimum for PROBLEM. PROBLEM is a
%     structure with the function FUN in PROBLEM.objective, the interval
%     in PROBLEM.x1 and PROBLEM.x2, the options structure in PROBLEM.options,
%     and solver name 'fminbnd' in PROBLEM.solver. The structure PROBLEM must
%     have all the fields.
%
%     [X,FVAL] = fminbnd(...) also returns the value of the objective function,
%     FVAL, computed in FUN, at X.
%
%     [X,FVAL,EXITFLAG] = fminbnd(...) also returns an EXITFLAG that describes
%     the exit condition of fminbnd. Possible values of EXITFLAG and the
%     corresponding exit conditions are
%
%      1  fminbnd converged with a solution X based on OPTIONS.TolX.
%      0  Maximum number of function evaluations or iterations reached.
%     -1  Algorithm terminated by the output function.
%     -2  Bounds are inconsistent (that is, ax > bx).
%
%     [X,FVAL,EXITFLAG,OUTPUT] = fminbnd(...) also returns a structure
%     OUTPUT with the number of iterations taken in OUTPUT.iterations, the
%     number of function evaluations in OUTPUT.funcCount, the algorithm name
%     in OUTPUT.algorithm, and the exit message in OUTPUT.message.
%
%     Examples
%       FUN can be specified using @:
%          X = fminbnd(@cos,3,4)
%        computes pi to a few decimal places and gives a message upon termination.
%          [X,FVAL,EXITFLAG] = fminbnd(@cos,3,4,optimset('TolX',1e-12,'Display','off'))
%        computes pi to about 12 decimal places, suppresses output, returns the
%        function value at x, and returns an EXITFLAG of 1.
%
%       FUN can be an anonymous function:
%          X = fminbnd(@(x) sin(x)+3,2,5)
%
%       FUN can be a parameterized function.  Use an anonymous function to
%       capture the problem-dependent parameters:
%          f = @(x,c) (x-c).^2;  % The parameterized function.
%          c = 1.5;              % The parameter.
%          X = fminbnd(@(x) f(x,c),0,1)
%
%
%     Reference page in Help browser
%        doc fminbnd
% 
% %% Lab exercises % %
%
%
1. Find the roots and the minimum values of the following functions. % Use function options and the problem structure to call the routines % fzero and fminbnd. %

%
%
1. $\cos(x) = x$ over $-2\pi \le x \le 2\pi$.
2. %
3. $\cos(5x) = x$ over $-2\pi \le x \le 2\pi$.
4. %
5. $\tan(x) = 0$ over $-2\pi \le x \le 2\pi$.
6. %
7. $\sin(1/x) = 0$ over $-\pi/50 \le x \le \pi/3$.
8. %
2. %
3. Consider the Colebrook equation for the friction factor in a fully % developed pipe flow % $$% \frac{1}{\sqrt{f}} = -2\log_{10}\left(\frac{\epsilon/D}{3.7} + % \frac{2.51}{\mbox{Re}_D \sqrt{f}}\right) %$$ % where $f$ is the Darcy friction factor, % $\epsilon/D$ is the relative roughness of % the pipe material, and $\mbox{Re}$ is the Reynolds % number based on a pipe diameter D. Use fzero to find % the roots of the Colebrook equation for the following combinations of % $\mbox{Re}$ and % $\epsilon/D$. %

%
%
1. $\epsilon/D = 0$, $\mbox{Re}_D = 5000$
2. %
3. $\epsilon/D = 0.0001$, $\mbox{Re}_D = 3 \times 10^5$
4. %
5. $\epsilon/D = 0.03$, $\mbox{Re}_D = 1 \times 10^4$
6. %
7. $\epsilon/D = 0.01$, $\mbox{Re}_D = 3 \times 10^5$
8. %
% (This problem was taken from G. Recktenwald, Numerical Methods with Matlab, Prentice % Hall, 2000. Page 289.) %
%
%