Back to tutorial index

Root-finding and minimization


Topics in this lab

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$.

Back to the top

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.
 
    See also roots, fminbnd, function_handle.
 
    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));
Solution is       :       1.5243452049841444
Residual error is :   0.00000000e+00

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));
Solution is       :       1.5243452049841444
Residual error is :   0.00000000e+00

Back to the top

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)
 
    See also optimset, fminsearch, fzero, function_handle.
 
    Reference page in Help browser
       doc fminbnd

Back to the top

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. $\cos(5x) = x$ over $-2\pi \le x \le 2\pi$.
    3. $\tan(x) = 0$ over $-2\pi \le x \le 2\pi$.
    4. $\sin(1/x) = 0$ over $-\pi/50 \le x \le \pi/3$.
  2. 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. $\epsilon/D = 0.0001$, $\mbox{Re}_D = 3 \times 10^5$
    3. $\epsilon/D = 0.03$, $\mbox{Re}_D = 1 \times 10^4$
    4. $\epsilon/D = 0.01$, $\mbox{Re}_D = 3 \times 10^5$
    (This problem was taken from G. Recktenwald, Numerical Methods with Matlab, Prentice Hall, 2000. Page 289.)

Back to the top

Get the code

Do you want to try the above code fragments on your own? Download the Matlab script that produces this page here. (lab_16.m)

Powered by MathJax

Published with MATLAB® 8.3