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.

Secant method

---

In this section, we consider a problem of finding root of the equation \( f(x) =0 \) for sufficiently smooth function f(x). If this equation has a solution, it is called a zero/null of the function f. Every rootfinding problem can be transformed into any number of fixed point problems. Consider the following example.

Example: Let \( f(x) = x^5 -5x+3 \) and we try to find its null, so we need to solve the equation \( x^5 -5x+3 =0 . \) Mathematica provides two real positive roots:

Solving for x, we convert the original equation into When you try the starting point \( x_0 = 0.5 , \) the fixed point iteration will provide you a very good approximation to the null x = 0.618034. However, when you start with \( x_0 = 1.5 , \) the fixed point iteration will diverge.

Let us add and subtract x from the equation: \( x^5 -5x+3+x =x . \) Expressing x, we derive another fixed point formula

Whatever the starting point was chosen, the above fixed point iteration diverges. ■

 

The formula involved in the secant method is very close to the one used in regula falsi:

When secant method is applied to find a square root of a positive number A, we get the formula

Example: Let us find a positive square root of 6. To start secant method, we need to pick up two first approximations,which we choose by obvious bracketing: \( x_0 =2, \quad x_1 =3 . \)

which is a number with 16 correct decimal places.

Example: We consider the function \( f(x) = e^x\, \cos x - x\, \sin x \) that has one root within the interval [0,3]. We define the range of x:

and then plot the function:

We make lower initial guess and upper initial guess: and then plot the function along with two guess bounds: We make the first iteration: Plot of the line through two guess points:

It should be noted that these two guesses do not have to bracket the root. We have called the two initial guesses in the way which we have because that will be the format for subsequent iterations. It does not matter which guess is which. Try switching the two initial guesses to see what happens. One should converge faster than the other.

We make the second iteration.

Then we plot bounds.

Third iteration gives.

 

II. Newton’s Method or Newton--Raphson Method

---

The following iterative method, named after the English mathematicians Isaac Newton (1642--1726) and Joseph Raphson (1648--1715), used for solving the equation f(x) = 0: The basic idea behind Newton's method is quite straightforward. Let p_n denote the most recent approximation to zero, p, of the function f. Replace f by its tangent line approximation that goes through the point p_n and takes the x-intercept of the tangent line as the next approximation p_n+1 to the root. ■

Newton's method can be realized with the aid of FixedPoint command:

Example: Suppose we want to find a square root of 4.5. First, we apply NewtonZero command:

which is actually 2.1213203435596424 ■

When Newton's method is applied to find a square root of a positive number A, we get the Babylonian formula Ancient Mesopotamia was a civilization that existed in the area of modern Turkey, Syria, Iraq and Iran, between the Mediterranean Sea and the Persian Gulf. In the period 1900--1600 BC, Babylon was the capital city of Mesopotamia and the mathematics recorded at this time came to be known as Babylonian Mathematics. It turns out that Babylonian mathematicians discovered Pythagorean theorem about 1000 years before the ancient Greek mathematician Pythagoras (570--495 BC).

The Babylonians had an accurate and simple method for finding the square roots of numbers. This method is also known as Heron's method, after the Greek mathematician who lived in the first century AD. Indian mathematicians also used a similar method as early as 800 BC. The Babylonians are credited with having first invented this square root method, possibly as early as 1900 BC.

Example: We use Newton's method to find a positive square root of 6. Starting with one of the two initial positions, we get

As we see it converges pretty fast. In fact, this method converges quadratically---the number of correct decimal places approximately doubles with every step! We can convert the square root problem to a fixed point iteration: which converges to \( \sqrt{6} = 2.449489742783178\ldots , \) but very slowly. However, another suggested iteration diverges.

Theorem: Let f be twice continuously differentiable function on the interval [a,b] with \( p \in \left( a, b \right) \quad\mbox{and} \quad f(p) =0 . \) Suppose that \( f' (p) \ne 0. \) Then there exists a positive number δ such that for any \( p_0 \in \left[ p - \delta , p+\delta \right] , \) the sequence \( \left\{ p_n \right\} \) generated by Newton's algorithm converges to p. ■

Sometimes Newton’s method does not converge; the above theorem guarantees that δ exists under certain conditions, but it could be very small. Newton's method is a good way of approximating solutions, but applying it requires some intelligence. You must beware of getting an unexpected result or no result at all. So why would Newton’s method fail? Well, the derivative may be zero at the root (so when the function at one of the iterated points will have zero slope); the function may fail to be continuously differentiable; one of the iterated points x_n is a local minimum/maximum of f; and you may have chosen a bad starting point, one that lies outside the range of guaranteed convergence. It may happen that in the Newton--Raphson method, an initial guess close to one root can jump to a location several roots away. Depending on the context, each one of these may be more or less likely. Degenerate roots (those where the derivative is 0) are "rare" in general and we do not consider this case.

Finally, there is a chance that Newton's method will cycle back and forth between two value and never converge at all. This failure is illustrated in Figure below.

For example, take a function \( f(x) = x^2 -2x +1 \) and let start with x₀ =0. The first iteration produces 1 and the second iteration returns back to 0, so the sequence will oscillate between the two without converging to a root.

Example: Consider the function \( f(x) = x\,e^{-x^2} \) that obviously has a root at x = 0. However, if your starting point exceeds the root of the equation \( f'(x) = 0 , \) which is \( x_0 = 1/\sqrt{2} \approx 0.707107 , \) than Newton's algorithm diverges.

Newton's method diverges for the cube root, which is continuous and infinitely differentiable, except for x = 0, where its derivative is undefined:

Example: Consider the cubic function \( f(x) = x^3 - 0.926\,x^2 + 0.0371\,x + 0.043 . \) First we plot the function, and then define the range of 'x' you want to see its null.

Taking initial guess we start Newton's iteration. Because this method uses a line tangent to the function at the initial guess, we must calculate the derivative of the function to find the slope of the line at this point. Here we will define the derivative of the function f(x) as g(x). We begin with first Newton's iteration. Then we plot tangent line along with given function and bounds for the root. Now we make the second iteration. Then we compare its value with the actual root obtained with standard Mathematica command. which is an approximation to more accurate root: 0.291095502695723. Indeed, This is zero within Mathmatica accuracy. So we see that two iterations are enough to get at least six accurate decimal places. For instance, the third iteration yields

Example: We reconsider the function \( f(x) = e^x\, \cos x - x\, \sin x \) that has one root within the interval [0,3]. We define the range of x:

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]. Let us take the middle of the interval as out first guess: which is 1.35. We set the maximum number of iterations to be 4: Then we define the derivative of the function f(x) which we denote by g(x). Now we find the actual solution using Mathematica (this is an approximation as well, but with good accuracy) We apply the Newton--Raphson method 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. We generate the sequence of root approximations. Absolute true error: Absolute relative true error: Absolute approximate error: Absolute relative approximate error: Significant digits at least correct. Graphs:

Example: This example shows that Newton's method may converge slowly due to an inflection point occurring in the vicinity of the root. Let us consider the function \( f(x) = (x-0.5)^3 . \) First we define the function and its derivative:

Then we plot the function:

Next we define the initial guess and maximum number of iterations. First iteration.

We start the second iteration. Then we plot tangent line.

Third iteration, and then fourth iteration. So we see that convergence is very slow. ■

 

III. Chebyshev Iteration Scheme

---

The Chebyshev iteration is an iterative method for determining the solutions of a system of linear equations \( f(x) =0. \) The method is named after Russian mathematician Pafnuty Chebyshev (1821--1894), who discoved it in 1838 as a student project (but not published until 1951). When Chebyshev iteration scheme is applied to solve a quadratic equation \( x^2 -A =0, \) it leads to

Example: Let us reconsider the problem of determination of a positive square root of 6 using Chebyshev iteration scheme:

So we see that the third iteration gives an approximation with 16 correct decimal places.

 

IV. Halley's Method

---

The English astronomer, geophysicist, mathematician, meteorologist, and physicist Edmund/Edmond Halley (1656--1741) discovered the following iteration algorithm in 1694. He is well known for first predicting the orbit of the comet that bears his name. Usually Halley's formula makes the convergence of the process of iteration extremely rapid when the first approximation is fairly close to the null. For example, application of Halley's method to solving the quadratic equation \( x^2 -A =0 \) leads to the recurrence

Example: Let us find a few iterations for determination of a square root of 6:

So we see that the third iteration gives an approximation with 16 correct decimal places.

 

V. Dekker's Method

---

Brent's method is a root-finding algorithm combining the bisection method, the secant method and inverse quadratic interpolation. It has the reliability of bisection but it can be as quick as some of the less reliable methods. The algorithm tries to use the potentially fast-converging secant method or inverse quadratic interpolation if possible, but it falls back to the more robust bisection method if necessary. Brent's method was discovered in 1973 by Richard Peirce Brent (born in 1946, Melbourne). He is an Australian mathematician and computer scientist. Richard is an emeritus professor at the Australian National University and a conjoint professor at the University of Newcastle (Australia). Brent's method is based on an earlier algorithm of 1969 proposed by Theodorus Jozef Dekker (born in 1927). Theodorus is a Dutch mathematician who completed his Ph.D degree from the University of Amsterdam in 1958. Consequently, the method is also known as the Brent--Dekker method.

The idea to combine the bisection method with the secant method goes back to Dekker (1969). Suppose that we want to solve the equation f(x) = 0

. As with the bisection method, we need to initialize Dekker's method with two points, say a₀ and b₀, such that \( f \left( a_0 \right) \quad\mbox{and} \quad f \left( b_0 \right) \) have opposite signs. If f is continuous on \( \left[ a_0 , b_0 \right] , \) the intermediate value theorem guarantees the existence of a solution between a₀ and b₀.

Three points are involved in every iteration:

• b_k is the current iterate, i.e., the current guess for the root of f;
• a_k is the "contrapoint," i.e., a point such that \( f \left( a_k \right) \quad\mbox{and} \quad f \left( b_k \right) \) have opposite signs, so the interval \( \left[ a_k , b_k \right] \) contains the solution. Furthermore, \( \left\vert f \left( b_k \right) \right\vert \) should be less than or equal to \( \left\vert f \left( a_k \right) \right\vert , \) so that b_k is a better guess for the unknown solution than a_k.
• b_k-1 is the previous iterate (for the first iteration, we set b_k-1 = a₀).

Two provisional values for the next iterate are computed. The first one is given by linear interpolation, also known as the secant method: and the second one is given by the bisection method If the result of the secant method, s, lies strictly between b_k and m, then it becomes the next iterate (b_k+1 = s), otherwise the midpoint is used (b_k+1 = m).

Then, the value of the new contrapoint is chosen such that \( f \left( a_{k+1} \right) \quad\mbox{and} \quad f \left( b_{k+1} \right) \) have opposite signs. If \( f \left( a_k \right) \quad\mbox{and} \quad f \left( b_{k+1} \right) \) have opposite signs, then the contrapoint remains the same: a_k+1 = a_k. Otherwise, \( f \left( b_{k+1} \right) \quad\mbox{and} \quad f \left( b_{k} \right) \) have opposite signs, so the new contrapoint becomes a_k+1 = b_k.

Finally, if \( \left\vert f \left( a_{k+1} \right) \right\vert < \left\vert f \left( b_{k+1} \right) \right\vert , \) then a_k+1 is probably a better guess for the solution than b_k+1, and hence the values of a_k+1 and b_k+1 are exchanged.

This ends the description of a single iteration of Dekker's method. ■

Dekker's method performs well if the function f is reasonably well-behaved. However, there are circumstances in which every iteration employs the secant method, but the iterates b_k converge very slowly (in particular, \( \left\vert b_k - b_{k-1} \right\vert \) may be arbitrarily small). Dekker's method requires far more iterations than the bisection method in this case.

 

VI. Brent's Method

---

Brent (1973) proposed a small modification to avoid this problem. He inserted an additional test which must be satisfied before the result of the secant method is accepted as the next iterate. Therefore, Brent's method is a hybrid method which combines the reliability of bracketing method and the speed of open methods.

Two inequalities must be simultaneously satisfied: Given a specific numerical tolerance δ if the previous step used the bisection method, the inequality

must hold to perform interpolation, otherwise the bisection method is performed and its result used for the next iteration.

If the previous step performed interpolation, then the inequality

is used instead to perform the next action (to choose) interpolation (when inequality is true) or bisection method (when inequality is not true).

Also, if the previous step used the bisection method, the inequality

must hold, otherwise the bisection method is performed and its result used for the next iteration. If the previous step performed interpolation, then the inequality is used instead.

Brent proved that his method requires at most N² iterations, where N denotes the number of iterations for the bisection method. If the function f is well-behaved, then Brent's method will usually proceed by either inverse quadratic or linear interpolation, in which case it will converge superlinearly.

Furthermore, Brent's method uses inverse quadratic interpolation instead of linear interpolation (as used by the secant method). If \( f(b_k ), \ f(a_k ) , \mbox{ and } f(b_{k-1}) \) are distinct, it slightly increases the efficiency. As a consequence, the condition for accepting s (the value proposed by either linear interpolation or inverse quadratic interpolation) has to be changed: s has to lie between \( \left( 3 a_k + b_k \right) /4 \) and b_k.