Preface

---

This tutorial was made solely for the purpose of education and it was designed for students taking Applied Math 0330. It is primarily for students who have very little experience or have never used Mathematica before and would like to learn more of the basics for this computer algebra system. As a friendly reminder, don't forget to clear variables in use and/or the kernel.

Finally, the commands in this tutorial are all written in bold black font, while Mathematica output is in normal font. This means that you can copy and paste all commands into Mathematica, change the parameters and run them. You, as the user, are free to use the scripts for your needs to learn the Mathematica program, and have the right to distribute this tutorial and refer to this tutorial as long as this tutorial is accredited appropriately.

Bracketing Methods

---

In this section we develop 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 [a,b], where f(a) 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 = r that lies somewhere within the interval [a,b]. One of the ways to test a numerical method for solving the equation \( f(x) =0 \) is to check its performance on the Wilkinson polynomials that have the roots \( \{ 1, 2, \ldots , n \} , \) but provide surprisingly tough numerical root-finding problems. For example,

 

I. The Bisection Method of Bolzano

---

The bisection method is a root-finding method that repeatedly bisects an interval and then selects a subinterval in which a root must lie so that the endpoints of the subinterval become closer and closer together until we obtain an interval of arbitrarily small width that brackets the zero. It is a very simple and robust method, 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.

The decision step for this process of interval halving is first to choose the midpoint \( c= (a+b)/2 \) and then to analyze the three possibilities that might arise:

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

If either of 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 [a,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 [a,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, [a,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

Theorem (Bolzano): Suppose that f is a continuous function on the closed interval [a,b] and that there exists a number r from this interval such that \( f(r) =0 . \) If f(a) and f(b) have opposite signs, and \( \left\{ c_n \right\}_{n\ge 0} \) represents the sequence of midpoints generated by the bisection algorithm, then

and the sequence of midpoints converges to the null r. ■

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: Consider the function \( f(x) = e^x\, \cos x - x\, \sin x \) that has three nulls on the interval [0,10], and Mathematica confirms

We bracket the first root by the interval [0,3] because \( f(0)\, f(3) \approx -20.3079 . \) Dividing it in half, we choose the middle point \( c_0 =1.5, \) with \( f(0)\, f(1.5) \approx -1.17922 . \) Since the latter product is negative, our next interval will be \( [a_1 , b_1 ] = [0 , 1.5] . \) Its middle point is \( c_1 = 0.75 . \) Since \( f(0.75)\, f(1.5) \approx -1.22374 , \) we choose \( [a_2 , b_2 ] = [0.75 , 1.5] . \) Its middle point is \( c_2 = 1.125 . \) Since \( f(1.125)\, f(1.5) \approx -0.369174 , \) we choose \( [a_3 , b_3 ] = [1.125 , 1.5] . \) Again, we halve it to obtain the midpoint \( c_3 = 1.3125 . \) Since \( f(1.125)\, f(1.3125) \approx -0.100152 , \) we choose \( [a_4 , b_4 ] = [1.125 , 1.3125] . \) We can continue this process, but we see that convergence to the zero r = 1.22539 of the function f is slow: only fifth iteration gave us the first correct decimal place after the period.

Now we show step by step how it works using Mathematica. First we plot the function to roughly identify the roots.

Then we define the range of x you want to see the function

We break the guess interval [0,3] into 10 uniform parts and calculate the values of the given function at each end point of these subintervals. So we see that the first root is within the interval [6/5, 3/2].

Now we check whether the center of this interval is closer to the root: Therefore, we find that the root is within the interval [6/5,27/20] or [1.2 , 1.35]. Finding the value of the function at the lower and upper guesses and the estimated root: We plot the root bracketing with two iterations:

Now we do the third iteration: We check the interval between which the root lies.

Here the bisection method algorithm is applied to generate the values of the roots, true error, absolute relative true error, approximate error, absolute relative approximate error, and the number of significant digits at least correct in the estimated root as a function of number of iterations.

Absolute true error: Absolute relative true error: Absolute approximate error: Absolute relative approximate error: Significant digits at least correct Graphs

We present some scripts for application of bisection method.

Then we apply the above code to find the null of the function \( f(x) = e^x\, \cos x - x\, \sin x \) using the bisection method: Another code: Another code:

 

II. False Position Method

---

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.

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 end points becomes fixed and the other does not march into the solution.

Example: We reconsider the function \( f(x) = e^x\, \cos x - x\, \sin x \) on interval [0,3]. From the previous example, we know that f has the null r = 1.22539. We ask Mathematica for help:

Using the above code, we conclude that after 100 iterations the approximate value is 1.20382, which is far away from the true one.

 

III. Ridders’ Method

---

A powerful variant on false position is due to Ridders (Ridders, C.J.F. 1979, IEEE Transactions on Circuits and Systems, vol. CAS-26, pp. 979–980). When a root is bracketed between x₁ and x₂. Ridders’ method first evaluates the function at the midpoint \( x_3 = \left( x_1 + x_2 \right) /2 . \) It then factors out that unique exponential function which turns the residual function into a straight line. Specifically, it solves for a factor e^Q that gives This is a quadratic equation in e^Q, which can be solved to give Now the false position method is applied, not to the values \( f(x_1 ), \ f(x_3 ), \ f(x_2 ) , \) but to the values \( f(x_1 ), \ f(x_3 )\, e^Q , \ f(x_2 )\, e^{2Q} , \) yielding a new guess for the root, x₄. The overall updating formula is This formula has some very nice properties. First, x₄ is guaranteed to lie in the interval \( \left( x_1 , x_2 \right) , \) so the 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.

 

IV. Location of Roots

---

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) . \) ■

 

 

 

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)