# Root-finding

## Topics in this lab

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

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

$\cos(x) = x$ over $-2\pi \le x \le 2\pi$. $\cos(5x) = x$ over $-2\pi \le x \le 2\pi$. $\tan(x) = 0$ over $-2\pi \le x \le 2\pi$. $\sin(1/x) = 0$ over $-\pi/50 \le x \le \pi/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$.

- $\epsilon/D = 0$, $\mbox{Re}_D = 5000$
- $\epsilon/D = 0.0001$, $\mbox{Re}_D = 3 \times 10^5$
- $\epsilon/D = 0.03$, $\mbox{Re}_D = 1 \times 10^4$
- $\epsilon/D = 0.01$, $\mbox{Re}_D = 3 \times 10^5$

*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_21.m)

Published with MATLAB® 9.0