Back to tutorial index

# Root-finding

## Introduction

Finding the roots of functions 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.

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 the Bisection, Secant, and Inverse Quadratic Interpolation (IQI) methods. 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

## Lab exercises

1. Find the roots of the following functions. Use function options and the problem structure to call the routines fzero.

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