Return to computing page for the first course APMA0330
Return to computing page for the second course APMA0340
Return to computing page for the fourth course APMA0360
Return to Mathematica tutorial for the first course APMA0330
Return to Mathematica tutorial for the second course APMA0340
Return to Mathematica tutorial for the fourth course APMA0360
Return to the main page for the course APMA0330
Return to the main page for the course APMA0340
Return to the main page for the course APMA0360
Return to Part III of the course APMA0330

Glossary

Bracketing Methods

---

In this section we develop in detail two bracketing methods for finding a zero/null of a real and continuous function f(x). All bracketing methods always converge, whereas open methods (discussed in the next section) may sometimes diverge. We must start with an initial interval [𝑎,b], where f(𝑎) and f(b) have opposite signs. Since the graph y = f(x) of a continuous function is unbroken, it will cross the abscissa at a zero x = α that lies somewhere within the interval [𝑎,b]. One of the ways to test a numerical method for solving the equation \( f(x) =0 \) is to check its performance on a polynomial whose roots are known. It is a custom to use some famous polynomials such as Chebyshev, Legendre or the Wilkinson polynomials: that have the roots \( \{ 1, 2, \ldots , n \} , \) but provide surprisingly tough numerical root-finding problems. For example,

 

Let f be a real single-valued function of a real variable. If f(α) = 0, then α is said to be a zero of f or null or, equivalently, a root of the equation f(x) = 0. It is customary to say that α is a root or zero of an algebraic polynomial f, but just a zero if f is not a polynomial. In this subsection, we discuss an algorithm for finding a root of a function, called the bisection method.

The bisection method is a very simple and robust algorithm, but it is also relatively slow. The method was invented by the Bohemian mathematician, logician, philosopher, theologian and Catholic priest of Italian extraction Bernard Bolzano (1781--1848), who spent all his life in Prague (Kingdom of Bohemia, now Czech republic). Most of Bolzano's works remained in manuscript and did not become noticed and therefore did not influence the development of the subject. Many of his works were not published until 1862 or later.

Assume that a continuous function f(x) on a given interval [𝑎, b] satisfies According to the intermediate value theorem, the function f(x) must have at least one root in [𝑎, b]. Usually [𝑎, b] is chosen to contain only one root α; but the following algorithm for the bisection method will always converge to some root α in [𝑎, b].

The bisection method requires two initial guesses 𝑎 = x₀ and b = x₁ satisfying the bracket condition f(x₀)·f(x₁) < 0. As the values of f(x₀) and f(x₁) are on opposite sides of the x-axis y = 0, the solution α at which f(α) = 0 must reside somewhere in between of these two guesses, i.e., x₀ < α < x₁. While any point between the two endpoints can be chosen for the next iteration, we want to avoid the worst possible case in which the solution always happens to be in the larger of the two sections. The bisection method suggests choosing the middle point in the interval [𝑎, b] as the next approximation to the root α. The decision step for this process of interval halving is first to choose the midpoint \( c= (a+b)/2 = \left( x_0 + x_1 \right) /2 \) and then to analyze the three possibilities that might arise:

Algorithm	Bisection (f, 𝑎, b, root, ϵ)
	Define c:= (𝑎 + b)/2. If b - c ≤ ϵ, then accept root:= c, and exit. If sign f(b) sign f(c) ≤ 0 (so f(b) and f(c) have opposite signs), then 𝑎:= c; otherwise b:= c. Return to step 1.

		plot = Plot[1*x^2, {x, 0, 1}, AspectRatio -> 1, PlotStyle -> {Thickness[0.01], Blue}, Axes -> False, Epilog -> {PointSize[0.03], Pink, Point[{{0, 0}, {2, 2.8}}]}, PlotRange -> {{-0.5, 1.2}, {-0.2, 1.4}}]; arrow = Graphics[{Arrowheads[0.07], Arrow[{{-0.2, 0.4}, {1.2, 0.4}}]}]; l1 = Graphics[{Thick, Dashed, Line[{{0, 0}, {0, 0.4}}]}]; l2 = Graphics[{Thick, Dashed, Line[{{0.5, 0.25}, {0.5, 0.4}}]}]; l3 = Graphics[{Thick, Dashed, Line[{{0.75, 0.54}, {0.75, 0.4}}]}]; l4 = Graphics[{Thick, Dashed, Line[{{1, 0.4}, {1, 1}}]}]; t1 = Graphics[{Black, Text[Style[Subscript[x, 0], 16], {0, 0.47}]}]; t2 = Graphics[{Black, Text[Style[Subscript[x, 1], 16], {1, 0.35}]}]; t3 = Graphics[{Black, Text[Style[Subscript[x, 2], 16], {0.5, 0.47}]}]; t4 = Graphics[{Black, Text[Style[Subscript[x, 3], 16], {0.75, 0.35}]}]; t5 = Graphics[{Black, Text[Style[\[Alpha], 16], {0.6, 0.45}]}]; t8 = Graphics[{Black, Text[Style["f(b)", {1.08, 1}]]}]; t9 = Graphics[{Black, Text[Style["f(x)", Bold, 16], {0.5, 0.7}]}]; Show[arrow, plot, l1, l2, l3, l4, t1, t2, t3, t4, t5, t8, t9, Epilog -> {PointSize[0.02], Purple, Point[{{0, 0}, {0, 0.4}, {0.5, 0.4}, {0.5, 0.25}, {0.75, 0.555}, {0.75, 0.4}, {1, 0.4}, {1, 1}}]}]
Bisection iteration.		Mathematica code

If either of the products \( f(a)\,f(c) \quad\mbox{or} \quad f(c)\,f(b) \) is negative, we have found an interval half as wide as the original interval that contains the root, and we are “squeezing down on it.” To continue the process, relabel the new smaller interval [𝑎, b] and repeat the process until the interval is as small as desired. Since the bisection process involves sequences of nested intervals and their midpoints, we will denote by \( [a_n , b_n ] \) the subinterval obtained on the n-th stage of the bisection algorithm. Therefore, \( [a_0 , b_0 ] \) is just the starting original interval [𝑎, b]. Its midpoint we denote by \( c_0 = (a+b)/2 . \) The next interval \( [a_1 , b_1 ] \) has one of the ends c₀ and the interval is half as wide as the starting one, [𝑎, b].

After arriving at the n-th interval \( [a_n , b_n ] , \) which brackets r and has midpoint c_n, the new interval \( [a_{n+1} , b_{n+1} ] \) is constructed, which also brackets r and is half as wide as \( [a_n , b_n ] . \) So

Although bisection is generally slower than other methods, the neatness of the procedure (its error analysis, in particular) is a positive feature that makes bisection attractive. This method is often used in the built-in root-finding routine of programmable calculators. One disadvantage of bisection is that it cannot be used to find roots when the function is tangent to the abscissa and does not pass through it.

Pseudocode for Bisection Method:

input a,v,eps
external f
fa = f(a)
fb = f(b)
if f(a)*f(b) > 0 then stop
n = fix((log(b-a) - log(eps))/log(2)) + 1
for i=1 to n do
c = a+0.5*(b-a)
fc = f(c)
if fa*fc < 0 then
b = c
fb = fc
else
if fa*fc > 0 then
a = c
fa = fc
else
alpha = c
return
endif
endif
endfor Example 1: bisection method find a root for cubic equation Example 2: \( f(x) = e^x\, \cos x - x\, \sin x \)

We present some scripts for the application of the bisection method. First, we type in a subroutine for halving the interval.

Then we apply this subroutine to perform bisection algorithm with tolerance ξ. Here is a demonstration of how to use the code above to find the null of the function \( f(x) = e^x\, \cos x - x\, \sin x \) based on the bisection method: Another code: Another code:

 

The poor convergence of the bisection method as well as its poor adaptability to higher dimensions (i.e., systems of two or more non-linear equations) motivate the use of better techniques. Another popular algorithm is the method of false position or the regula falsi method. A better approach is obtained if we find the point c on the abscissa where the secant line crosses it. To determine the value of c, we write down two versions of the slope m f the line that connects the points \( (a, f(a)) \quad\mbox{and} \quad (b, f(b)) : \)

where points \( (a, f(a)) ,\quad (b, f(b)) , \quad\mbox{and} \quad (c,0) \) were used. Solving for c, we get The three possibilities are the same as before:

• If f(a) and f(c) have opposite signs, a zero lies in [a,c] .
• If f(b) and f(c) have opposite signs, a zero lies in [c,b] .
• If \( f(c) =0 , \) then the zero is c.

		plot = Plot[1*x^2, {x, 0, 1}, AspectRatio -> 1, PlotStyle -> {Thickness[0.01], Blue}, Axes -> False, Epilog -> {PointSize[0.03], Pink, Point[{{0, 0}, {1, 1}}]}, PlotRange -> {{-0.2, 1.2}, {-0.2, 1.4}}]; line = Graphics[{Thick, Red, Line[{{0, 0}, {1, 1}}]}]; arrow = Graphics[{Arrowheads[0.07], Arrow[{{-0.2, 0.3}, {1.2, 0.3}}]}]; l1 = Graphics[{Thick, Dashed, Line[{{0, -0.1}, {0, 0.3}}]}]; l3 = Graphics[{Thick, Dashed, Line[{{1, 0.3}, {1, 1}}]}]; t1 = Graphics[{Black, Text[Style["a", 16], {0, 0.36}]}]; t3 = Graphics[{Black, Text[Style["f(a)", 16], {-0.08, 0.0}]}]; t4 = Graphics[{Black, Text[Style[\[Alpha], 16], {0.55, 0.23}]}]; t5 = Graphics[{Black, Text[Style["c", 16], {0.3, 0.36}]}]; t6 = Graphics[{Black, Text[Style["b", 16], {1, 0.23}]}]; t8 = Graphics[{Black, Text[Style["f(b)", {1.08, 1}]]}]; t9 = Graphics[{Black, Text[Style["f(x)", Bold, 16], {0.87, 0.56}]}]; Show[arrow, plot, line, l1, l3, t1, t3, t4, t5, t6, t8, t9, Epilog -> {PointSize[0.02], Purple, Point[{{0, 0}, {0, 0.3}, {0.3, 0.3}, {0.55, 0.3}, {1, 1}, {1, 0.3}}]}]
Regula Falsi.		Mathematica code

The decision process is the same as in bisection method and we generate a sequence of intervals \( \left\{ [a_n , b_n ]\right\} \) each of which brackets the null. At each step the approximation of the zero r is This method is called regular falsi (false position method) and it was known to ancient Babylonian mathematicians. It has also linear rate of convergence of the sequence \( \left\{ c_n \right\}_{n\ge 0} \) to the null r of the function f. But beware: although the integral width \( b_n - a_n \) is getting smaller, it is possible that it may not go to zero. If the graph of \( y= f(x) \) is concave near (r,0), one of the endpoints becomes fixed and the other does not march into the solution. Example 3: regula falsi for \( f(x) = e^x\, \cos x - x\, \sin x \)

 

A powerful variant on false position was proposed in 1979 by C.J.F. Ridders. The idea of the method is to apply the regula falsi to the equation h(x) = 0 instead of f(x) = 0, where \( h(x) = f(x)\, e^{mx} \) because both functions,f and h, have the same roots. In many cases, multiplication by an exponential function may facilitate determination of the root of f(x) because its bracketing could be highlighted. One might expect a faster convergence of this method because of a better approximation of the function f(x). However, there is a price to pay: the algorithm is a bit more complex and an additional calculation of the function beyond the bracketing points is required.

		plot = Plot[1*x^2, {x, 0, 1}, AspectRatio -> 1, PlotStyle -> {Thickness[0.01], Blue}, Axes -> False, PlotRange -> {{-0.2, 1.2}, {-0.2, 1.4}} ] line = Graphics[{Thick, Red, Line[{{0, -0.1}, {1, 0.8}}]}, Epilog -> {PointSize[0.03], Pink, Point[{{0.5, 0.5}}]}]; arrow = Graphics[{Arrowheads[0.07], Arrow[{{-0.2, 0.3}, {1.2, 0.3}}]}]; l1 = Graphics[{Thick, Dashed, Line[{{0, -0.1}, {0, 0.3}}]}]; l2 = Graphics[{Thick, Dashed, Line[{{0.5, 0.25}, {0.5, 0.35}}]}]; l3 = Graphics[{Thick, Dashed, Line[{{1, 0.3}, {1, 1}}]}]; t1 = Graphics[{Black, Text[Style[Subscript[x, 1], 16], {0, 0.36}]}]; t2 = Graphics[{Black, Text[Style[Subscript[h, 1], 16], {0, -0.16}]}]; t3 = Graphics[[{Black, Text[Style[Subscript[f, 1], 16], {-0.08, 0.0}]}]; t4 = Graphics[{Black, Text[Style["u", 16], {0.5, 0.17}]}]; t5 = Graphics[{Black, Text[Style[Subscript[x, 3], 16], {0.44, 0.36}]}]; t6 = Graphics[{Black, Text[Style[Subscript[x, 2], 16], {1, 0.26}]}]; t7 = Graphics[{Black, Text[Style[Subscript[h, 2], 16], {1.08, 0.8}]}]; t8 = Graphics[{Black, Text[Style[Subscript[f, 2], 16], {1.08, 1}]}]; t9 = Graphics[{Black, Text[Style["f(x)", Bold, 16], {0.75, 0.76}]}]; Show[arrow, plot, line, l1, l2, l3, t1, t2, t3, t4, t5, t6, t7, t8, t9, Epilog -> {PointSize[0.02], Purple, Point[{{0, -0.1}, {0, 0}, {0, 0.3}, {0.45, 0.3}, {0.5, 0.3}, {1, 0.8}, {1, 1}, {1, 0.3}}]}]
Ridders’ algorithm.		Mathematica code

When a root is bracketed between x₁ and x₂, Ridders’ method first evaluates the function h(x) at the midpoin \( u = \left( x_1 + x_2 \right) /2 . \) Then for three equidistance x-values x₁, u, and x₂, the following requirement for the straight line connecting points (x₁, f₁) and (x₂, f₂) is met:

where Let \( d = \left( x_2 - x_1 \right) /2 = u - x_1 = x_2 - u \) and the bracketing condition f₁·f₂ < 0 holds, then from the requirement condition we derive Since this is a quadratic equation in e^m d, it can be analytically solved to give where sign(x) stands for the sign of the function’s argument: Finally, application of the false position method to the points h_u and either h₁ or h₂ (depending which pair brackets the root) yields a new guess for the root, x₃ This formula has some very nice properties. First, x₃ is guaranteed to lie in the interval \( \left( x_1 , x_2 \right) , \) so the Ridders’ method never jumps out of its brackets. Second, the convergence of successive applications of the above formula is quadratic. Since each application of Ridders’ formula requires two function evaluations, the actual order of the method is \( \sqrt{2} , \) not 2; but this is still quite respectably superlinear: the number of significant digits in the answer approximately doubles with each two function evaluations. Third, taking out the function’s “bend” via exponential (that is, ratio) factors, rather than via a polynomial technique (e.g., fitting a parabola), turns out to give an extraordinarily robust algorithm. In both reliability and speed, Ridders’ method is generally competitive with the more highly developed and better established (but more complicated) method of Van Wijngaarden, Dekker, and Brent, which we discuss next.

Algorithm:

Ridders’ rule for determination of the root of the equation f(x) = 0

Find proper bracketing points x1 and x2 for the root α of the equation f(x) = 0. It is irrelevant which is positive end and which is negative, but the function must have different signs at these points, i.e., f(x1)×f(x2) < 0. Start loop. Find the midpoint \( u = \left( x_2 - x_1 \right) . \) Find a new approximation for the root \[ x_3 = \frac{x_1 + x_2}{2} + \frac{x_2 - x_1}{2} \cdot \frac{\mbox{sign} \left[ f(x_1 ) - f(x_2 ) \right] f(u )}{\sqrt{f^2 (u ) - f(x_1 )f(x_2 )}} , \qquad\mbox{with} \quad u = \frac{x_1 + x_2}{2} . \] If x3 satisfies the required precision and the convergence condition is reached, then stop. Rebracket the root, i.e., assign new x1 and x2, using old values: one end of the bracket is x3 and f3 = f(x3); the other is whichever of (x1, x2, u) is closer to x3 and provides the proper bracket. Repeat the loop.

Example 4: Ridders’ rule to find a root of the Bessel function

 

To roughly estimate the locations of the roots of the equation \( f(x) =0 \) over the interval [a,b], it make sense to make a uniform partition of the interval by sample points \( \left( x_n , y_n = f(x_n )\right) \) and use the following criteria:

That is, either \( f(x_{n-1}) \quad\mbox{and} \quad f(x_n ) \) have opposite signs or \( |f(x_n )| \) is small and slope of the curve \( y= f(x) \) changes sign near \( \left( x_n , f(x_n )\right) . \) ■

 

 

References

 

Return to Mathematica page
Return to the main page (APMA0330)
Return to the Part 1 (Plotting)
Return to the Part 2 (First Order ODEs)
Return to the Part 3 (Numerical Methods)
Return to the Part 4 (Second and Higher Order ODEs)
Return to the Part 5 (Series and Recurrences)
Return to the Part 6 (Laplace Transform)
Return to the Part 7 (Boundary Value Problems)