3

Numerical Solution of Equations of a Single Variable

This chapter focuses on numerical solution of equations of a single variable, which appear
in the general form

f ( x) = 0 (3.1)

Graphically, a solution (or root) of f(x) = 0 refers to the point of intersection of f(x) and
the x-axis. Therefore, depending on the nature of the graph of f(x) in relation to the x-axis,
Equation 3.1 may have a unique solution, multiple solutions, or no solution. A root of an
equation can sometimes be determined analytically resulting in an exact solution in closed
form. For instance, the equation e3x − 2 = 0 can be solved analytically to obtain a unique
solution x = 31 ln 2 . In most situations, however, this is not possible and the root(s) must be
found numerically. As an example, consider the equation 2 − x + cos x = 0. Figure 3.1 shows
that this equation has one solution only, which may be approximated to within a desired
accuracy with the aid of a numerical method.

3.1 Numerical Solution of Equations

As described in Figure 3.2, numerical methods for solving an equation are divided into
three main categories: bracketing methods, open methods, and using the built-in MATLAB
function fzero.
Bracketing methods require that an interval containing the root irst be identiied.
Referring to Figure 3.3, this means an interval [a, b] with the property that f(a)f(b) < 0 so that
a root lies in [a, b]. The length of the interval is then reduced in succession until a desired
accuracy is satisied. Exactly how this interval gets narrowed in each step depends on the
speciic method used. It is readily seen that bracketing methods always converge to the root.
Open methods require an initial estimate of the solution, somewhat close to the intended
root. Subsequently, more accurate estimates are successively generated by a speciic method;
Figure 3.4. Open methods are more eficient than bracketing methods, but do not always
generate a sequence that converges to the root. The built-in function fzero inds the root of
a function f near a speciied point, or in a speciied interval [a, b] such that f(a)f(b) < 0.

3.2 Bisection Method

The bisection method is the simplest bracketing method to ind a root of f(x) = 0. It is
assumed that f(x) is continuous on an interval [a, b] and has a root there so that f(a) and f(b)
55
56 Numerical Methods for Engineers and Scientists Using MATLAB®, Second Edition

1
2 – x + cos(x)

–1

–2

–3 Approximate root
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
x

FIGURE 3.1
Approximate solution of 2 − x + cos x = 0.

Bisection method

Bracketing methods
Regula Falsi method

Numerical solution of f(x) = 0 Fixed-point method

Open methods Newton’s method

Secant method
MATLAB function fzero

FIGURE 3.2
Classiication of methods to solve an equation of one variable.

have opposite signs, hence f(a)f(b) < 0. The procedure goes as follows: Locate the midpoint
of [a, b], that is, c1 = 12 ( a + b), Figure 3.5. If f(a) and f(c1) have opposite signs, the interval [a, c1]
contains the root and will be retained for further analysis; that is, the left end is retained
while the right end is adjusted. If f(b) and f(c1) have opposite signs, we continue with [c1, b];
that is, the right end is kept while the left end is adjusted. In Figure 3.5 it so happens that
the interval [c1, b] brackets the root and is retained. Since the right endpoint is unchanged,
we update the interval [a, b] by resetting the left endpoint a = c1. The process is repeated
until the length of the most recent interval [a, b] satisies the desired accuracy.
The initial interval [a, b] has length b − a. Beyond that, the irst generated interval has
length 12 (b − a), the next interval 14 (b − a), and so on. Thus, the n-th interval constructed
Numerical Solution of Equations of a Single Variable 57

f (x)

Root b
x
0
a
Initial interval

Root b
x
0
a
Second interval

Root b
x
0
a
Third interval

FIGURE 3.3
Philosophy of bracketing methods.

in this manner has length (b − a)/2n−1, and because it brackets the root, the absolute error
associated with the nth iteration satisies

b−a
|en |≤ ( b > a)
2 n −1

This upper bound is usually larger than the actual error at the nth step. If the bisection
method is used to approximate the root of f(x) = 0 within a prescribed tolerance ε > 0, then
it can be shown that the number N of iterations needed to meet the tolerance condition
satisies

ln(b − a) − ln ε
N> (3.2)
ln 2

The user-deined function Bisection shown below generates a sequence of values

(midpoints) that ultimately converges to the true solution. The iterations terminate when
2 (b − a) < ε, where ε is a prescribed tolerance. The output of the function is the last gener-
1

ated estimate of the root at the time the tolerance was met. It also returns a table that com-
prises iteration counter, interval endpoints, and interval midpoint per iteration, as well as
the value of 12 (b − a) to see when the terminating condition is satisied.
58 Numerical Methods for Engineers and Scientists Using MATLAB®, Second Edition

Initial estimate

Root
x
0
f(x)

Second estimate

Root
x
0

Third estimate

Root
x
0

FIGURE 3.4
Philosophy of open methods.

f(x)

Root b
x
0 a

c1
b
First estimate
a c1 = 12 (a + b)
Adjust interval
a = c1
c2
b
Second estimate
a c2 = 12 (a + b) = 1
2 (c1 + b)
Adjust interval
b = c2

b
Third estimate
a c3 c3 = 12 (a + b) = 12 (c1 + c2)

FIGURE 3.5
Bisection method: three iterations shown.
Numerical Solution of Equations of a Single Variable 59

function c = Bisection(f, a, b, kmax, tol)

%
% Bisection uses the bisection method to find a root of f(x) = 0
% in the interval [a,b].
%
% c = Bisection(f, a, b, kmax, tol), where
%
% f is an anonymous function representing f(x),
% a and b are the endpoints of interval [a,b],
% kmax is the maximum number of iterations (default 20),
% tol is the scalar tolerance for convergence (default 1e-4),
%
% c is the approximate root of f(x) = 0.
%
if nargin < 5 || isempty(tol), tol = 1e-4; end
if nargin < 4 || isempty(kmax), kmax = 20; end
if f(a)*f(b) > 0
c = 'failure';
return
end
disp(' k a b c (b-a)/2')
for k = 1:kmax,
c = (a+b)/2; % Find the first midpoint
if f(c) == 0, % Stop if a root has been found
return
end
fprintf('%3i %11.6f%11.6f%11.6f%11.6f\n',k,a,b,c,(b-a)/2)
if (b-a)/2 < tol, % Stop if tolerance is met
return
end
if f(b)*f(c) > 0 % Check sign changes
b = c; % Adjust the endpoint of interval
else a = c;
end
end

EXAMPLE 3.1: BISECTION METHOD

The equation x cos x + 1 = 0 has a root in the interval [−2, 4], as shown in Figure 3.6:

>> f = @(x)(x*cos(x)+1);
>> ezplot(f,[-2,4])

Deine f(x) = x cos x + 1 so that f(−2) > 0 and f(4) < 0. We will perform two steps of the
bisection method as follows. The irst midpoint is found as c1 = 12 (−2 + 4) = 1. Since the
function value at this point is f(c1) = f(1) > 0, the root must be in [1, 4]. This means the left
end is adjusted as a = c1 = 1 while b = 4 remains unchanged. The next midpoint is then
calculated as c2 = 12 (1 + 4) = 2.5 . Since f(c2) = f(2.5) < 0, the root must lie in [1, 2.5]. This
process continues until a desired accuracy is achieved. In particular, if we execute the
user-deined function Bisection with ε = 10−2 and maximum 20 iterations, the follow-
ing results are obtained.
60 Numerical Methods for Engineers and Scientists Using MATLAB®, Second Edition

1.5

0.5
Root
x cos(x) +1

–0.5

–1

–1.5

–2

–2.5
–2 –1 0 1 2 3 4
x

FIGURE 3.6
Location of the root of x cos x + 1 − 0 in [−2, 4].

>> c = Bisection(f, -2, 4, [], 1e-2)

k a b c (b-a)/2
1 -2.000000 4.000000 1.000000 3.000000
2 1.000000 4.000000 2.500000 1.500000
3 1.000000 2.500000 1.750000 0.750000
4 1.750000 2.500000 2.125000 0.375000
5 1.750000 2.125000 1.937500 0.187500
6 1.937500 2.125000 2.031250 0.093750
7 2.031250 2.125000 2.078125 0.046875
8 2.031250 2.078125 2.054688 0.023438
9 2.054688 2.078125 2.066406 0.011719
10 2.066406 2.078125 2.072266 0.005859
c =
2.0723

The data in the irst three rows conirm the hand calculations. Iterations stopped
when 12 (b − a) = 0.005859 < ε = 0.01 . Note that by Equation 3.2, we have

ln(b − a) − ln ε a =−2 , b = 4 ln 6 − ln 0.01

N> = = 9.23
ln 2 ε = 0.01 ln 2

which means at least 10 iterations are required for convergence. This is in agree-
ment with the indings here, as we saw that tolerance was met after 10 iterations.
The accuracy of the solution estimate will improve if a smaller tolerance is imposed.

3.2.1 MATLAB Built-In Function fzero

The fzero function in MATLAB inds the roots of f(x) = 0 for a real function f(x).
FZERO Scalar nonlinear zero finding.

X = FZERO(FUN,X0) tries to find a zero of the function FUN near X0,

if X0 is a scalar.
Numerical Solution of Equations of a Single Variable 61

The fzero function uses a combination of the bisection, secant, and inverse quadratic
interpolation methods. If we know two points where the function value differs in sign, we
can specify this starting interval using a two-element vector for x0. This algorithm is guar-
anteed to return a solution. If we specify a scalar starting point x0, then fzero initially
searches for an interval around this point where the function changes sign. If an interval
is found, then fzero returns a value near where the function changes sign. If no interval
is found, fzero returns a NaN value.
The built-in function fzero can be used to conirm the approximate root in Example 3.1.
This can be done in one of two ways. As a irst option, we may specify a point near which
the root must be found. For instance, selecting x0 = 1 leads to

>> fzero(f,1)

ans =

2.0739

As a second option, we may identify two points where the function value differs in sign.
For instance, choosing the interval [1, 3] leads to

>> fzero(f,[1,3])

ans =

2.0739

The accuracy of the approximate root (2.0723) returned by the user-deined func-
tion Bisection can be improved by choosing a smaller tolerance. For example, the
reader can verify that executing Bisection with tol = 1e-8 returns a root esti-
mate (2.0739) that agrees with fzero to at least 4 decimal places, but requires 30
iterations.