Root Finding

Root finding is a fundamental problem in numerical analysis and has many applications in science and engineering such as solving nonlinear equations, optimization problems, and differential equations. Usually, a closed form of the root is not available, and we need to compute the root numerically. In this chapter, we will discuss some of the most common methods for root finding.

Bracket Methods

If \(f\) is a continuous function, and \(f(a)\) and \(f(b)\) have opposite signs, then by the Intermediate Value Theorem, there exists a root of \(f\) on the interval \([a, b]\). The bracket method is based on this fact and iteratively locates the pair of points \(a\) and \(b\) such that \(f(a)\) and \(f(b)\) have opposite signs. The most common bracket methods are the bisection method and the false position method.

Bisection Method

The simplest bracket method is the bisection method. Once \(f(a)f(b) < 0\), one can select the midpoint \(c = \frac{a + b}{2}\) and check the sign of \(f(a) f(c)\).

• If \(f(a)f(c) < 0\), then the root is in the interval \([a, c]\).

• If \(f(a)f(c) > 0\), then the root is in the interval \([c, b]\).

• If \(f(a)f(c) = 0\), then \(c\) is the root.

Iterative vs Recursive

The bisection method can be implemented either iteratively or recursively. The recursive program usually is more compact, but may suffer from inefficiency (more instructions) and the risk of stack overflow.

For a variety of programming languages, recursive implementation can be made more efficient by using the so called tail recursion optimization. This is a compiler feature that allows the recursive program to be executed with the same efficiency as the iterative one. However, neither Python nor MATLAB natively supports tail recursion optimization. : )

For the first two cases, we can repeat the process with the new interval until certain stop criteria are met. Each iteration reduces the size of the interval by half (gaining one bit each iteration), the total number of iterations required to reduce the interval to a certain size is \(\lceil\log_2\left(\frac{b - a}{\epsilon}\right)\rceil\), where \(\epsilon\) is the desired tolerance.

Once the function \(f\) has a sign change over the bracket \([a, b]\), the bisection method is guaranteed to converge to a root, but it is not very efficient. It is usually used to obtain a rough estimate of the root, which is then taken as an initial guess for a more efficient method.

False Position Method

The bisection method only uses \(\text{sgn}(f(a))\) and \(\text{sgn}(f(b))\) instead of the function values. The false position method (Regula falsi in Latin) improves the bisection method by taking the function values into account. Instead of selecting the midpoint \(c = \frac{a + b}{2}\), the false position method selects the point \(c\in[a, b]\) that lies on the line connecting \((a, f(a))\) and \((b, f(b))\), that is

\[c = \frac{f(b) a - f(a) b}{f(b) - f(a)} = a - f(a)\frac{b - a}{f(b) - f(a)}.\]

The false position method is also guaranteed to converge to a root if \(f(a)f(b) < 0\) and the implementation is quite similar to the bisection method. It usually converges faster than the bisection method, but sometimes exceptions occur.

Note

When \(f''\) keeps the same sign over \([a, b]\), it is not hard to show that only one side of the bracket is updating. The bracket size will never decrease to zero, which is different from the bisection method. In the following, we use an example to illustrate this, see Fig. 1. The function \(f(x) = x^2 - 1\) over the initial bracket \([0, 2]\). The left endpoint is updating to the root while the right endpoint is fixed at \(2\).

Fig. 1 False position method

It is easy to improve the false position method by forcing more weight toward the other endpoint. This is called the Illinois method. Once the same side has updated in two consecutive iterations, the Illinois method will adjust \(c\) using a slightly different formula.

\[ {c} = a - f(a)\frac{ b - a }{\lambda f(b) - f(a)} \]

or

\[ {c} = b - f(b) \frac{a - b}{\lambda f(a) - f(b)}, \]

where the weight \(\lambda\) controls the position of \(c\). The weight is initially set to \(1\) which corresponds to the standard false position method. If the same side is about to update twice, the weight of the other side will be halved. If not, the weight on the new bracket point will be reset to \(1\). See Algorithm 2.

The choice of decay factor \(\frac{1}{2}\) is optimal if it has to be a constant (explain later). The factor can be replaced with other variable values. A usual replacement is the Pegasus method.

Essentially, the Pegasus method replaces \(\lambda = 1/2\) with \(\lambda=\frac{f_1}{f_1 + f_2}\).

Algorithm 2 (Illinois Method)

Inputs \(f\), \(a, b\), \(\epsilon\)

Outputs \(c\)

1. \(x_0\gets a\), \(x_1\gets b\), \(f_0\gets f(x_0)\), \(f_1\gets f(x_1)\). \(iter = 0\) //initialization.

2. While True do:

   1. \(x_2 \gets x_0 - f_0 \frac{x_1 - x_0}{f_1 - f_0}\), \(f_2 \gets f(x_2)\), \(iter \gets iter + 1\). //standard false position step, \(x_1\) and \(x_2\) are two latest iterations.

   2. If \( |f_2| < \epsilon\), then return \(x_2\); //check stopping criteria

   3. While \(f_1 f_2 > 0\), then // adjust until the sign changes

      • \((x_0, f_0)\gets (x_0, \lambda f_0)\), where \(\lambda = \frac{1}{2}\)

      • \((x_1, f_1)\gets (x_2, f_2)\) and \(x_2 \gets x_0 - f_0 \frac{x_1 - x_0}{f_1 - f_0}\), \(f_2 \gets f(x_2)\), \(iter \gets iter + 1\).

      • If \( |f_2| < \epsilon\), then return \(x_2\); //check stopping criteria

   4. If \(f_1 f_2 < 0\), then // perform a false position step

      • \((x_0, f_0)\gets (x_1, f_1)\) and \((x_1, f_1)\gets (x_2, f_2)\).

Note

We use the previous example to illustrate the difference between the false position method and the Illinois method. At the second iteration, the Illinois method finds the updating is still on left side, so it modifies right endpoint \(f(b)\) into \(\frac{1}{2} f(b)\) to compute the new \(c\), which makes the selected point \({c}\) closer to the right endpoint than the false position method, see Fig. 2.

Fig. 2 Illinois method

Stopping Criteria

There are several types of stopping criteria to terminate the iteration. Common ones include:

• \(|f(c)| < \texttt{ftol}\), tolerance on the function value.

• \(|c_n - c_{n-1}| < \texttt{atol}\), absolute tolerance.

• \(|c_n - c_{n-1}| < \texttt{rtol} |c_n|\), relative tolerance.

Example 2

Let us try the aforementioned methods to find the root of \(f(x) = x^3 - 2x^2 - 4\) on the interval \([1, 3]\). The root \(x^{\ast}\) can be computed analytically through cubic root formula, which is roughly 2.5943130163548496.

Using the previous methods, we obtain the sequence of selection \(c_n\) and the error \(|c_n - x^{\ast}|\). The tolerance is set to \(|f(c)| < 10^{-6}\). The results are shown in the following table.

Bisection Method

iter  1 | 2.0000000000000000000000000 | 5.94e-01 
iter  2 | 2.5000000000000000000000000 | 9.43e-02 
iter  3 | 2.7500000000000000000000000 | 1.56e-01 
iter  4 | 2.6250000000000000000000000 | 3.07e-02 
iter  5 | 2.5625000000000000000000000 | 3.18e-02 
iter  6 | 2.5937500000000000000000000 | 5.63e-04 
iter  7 | 2.6093750000000000000000000 | 1.51e-02 
iter  8 | 2.6015625000000000000000000 | 7.25e-03 
iter  9 | 2.5976562500000000000000000 | 3.34e-03 
iter 10 | 2.5957031250000000000000000 | 1.39e-03 
iter 11 | 2.5947265625000000000000000 | 4.14e-04 
iter 12 | 2.5942382812500000000000000 | 7.47e-05 
iter 13 | 2.5944824218750000000000000 | 1.69e-04 
iter 14 | 2.5943603515625000000000000 | 4.73e-05 
iter 15 | 2.5942993164062500000000000 | 1.37e-05 
iter 16 | 2.5943298339843750000000000 | 1.68e-05 
iter 17 | 2.5943145751953125000000000 | 1.56e-06 
iter 18 | 2.5943069458007812500000000 | 6.07e-06 
iter 19 | 2.5943107604980468750000000 | 2.26e-06 
iter 20 | 2.5943126678466796875000000 | 3.49e-07 
iter 21 | 2.5943136215209960937500000 | 6.05e-07 
iter 22 | 2.5943131446838378906250000 | 1.28e-07 
iter 23 | 2.5943129062652587890625000 | 1.10e-07 
iter 24 | 2.5943130254745483398437500 | 9.12e-09 

False Position Method

iter  1 | 2.0000000000000000000000000 | 5.94e-01 
iter  2 | 2.4444444444444446418174266 | 1.50e-01 
iter  3 | 2.5621621621621617492792211 | 3.22e-02 
iter  4 | 2.5876913365185605364615640 | 6.62e-03 
iter  5 | 2.5929610854818996301673906 | 1.35e-03 
iter  6 | 2.5940374914642010395482430 | 2.76e-04 
iter  7 | 2.5942568846837747997824408 | 5.61e-05 
iter  8 | 2.5943015817106331866170876 | 1.14e-05 
iter  9 | 2.5943106870264029950590157 | 2.33e-06 
iter 10 | 2.5943125418534931370118102 | 4.75e-07 
iter 11 | 2.5943129196954899384763849 | 9.67e-08 

Illinois Method

iter  1 | 2.0000000000000000000000000 | 5.94e-01 
iter  2 | 2.6153846153846154187760931 | 2.11e-02 
iter  3 | 2.5847750865051901669744439 | 9.54e-03 
iter  4 | 2.5941951587569969106539247 | 1.18e-04 
iter  5 | 2.5944267005726100450146987 | 1.14e-04 
iter  6 | 2.5943130084597889606357057 | 7.90e-09 

More Digits

Using bigfloat package, it is possible to extract more digits to observe the dynamics of the convergence. For instance, by setting the precision to 900 bits, the Illinois method’s convergence is shown in the following table.

iter  1 | 2.0000000000000000000000000 | 5.94e-01 
iter  2 | 2.6153846153846154187760931 | 2.11e-02 
iter  3 | 2.5847750865051901669744439 | 9.54e-03 
iter  4 | 2.5941951587569973547431346 | 1.18e-04 
iter  5 | 2.5944267005726091568362790 | 1.14e-04 
iter  6 | 2.5943130084597889606357057 | 7.90e-09 
iter  7 | 2.5943130163543197674869134 | 5.29e-13 
iter  8 | 2.5943130163553775879847763 | 5.29e-13 
iter  9 | 2.5943130163548486777358448 | 1.65e-25 
iter 10 | 2.5943130163548486777358448 | 5.13e-38 
iter 11 | 2.5943130163548486777358448 | 5.13e-38 
iter 12 | 2.5943130163548486777358448 | 1.55e-75 
iter 13 | 2.5943130163548486777358448 | 4.70e-113 
iter 14 | 2.5943130163548486777358448 | 4.70e-113 
iter 15 | 2.5943130163548486777358448 | 1.30e-225 

Remark 2

The bracket methods need to first locate an interval \([a, b]\) such that \(f(a)f(b) < 0\). A common approach is to sample a few equidistant points in a large interval and then use the sign of the function values to identify the bracket. This is a simple and robust approach, but it may require a large number of function evaluations.

Order of Convergence

The order of convergence quantifies how fast the sequence approximates the limiting value.

More About Convergence

In practice, the limit or even upper bound may not exist for \(\frac{|x_{n+1} - x^{\ast}|}{|x_n - x^{\ast}|^p}\), it is possible to consider the convergence rate in weaker sense. For instance, for sufficiently large \(n\) that the inequality

\[ \lim_{k\to\infty}\sqrt[k]{\frac{|x_{n+k} - x^{\ast}|}{|x_n - x^{\ast}|^{p^k}}}= \rho \]

holds for certain \(p > 0\) and \(\rho > 0\), then the sequence has a mean order of convergence is \(p\) and a mean convergence rate \(\rho\).

On average, each iteration contributes an order of \(p\) and a rate of \(\rho\).

Definition 2

The order of convergence of a sequence \(\{x_n\}\) is \(p > 0\) if

\[ \lim_{n\to\infty}\frac{|x_{n+1} - x^{\ast}|}{|x_n - x^{\ast}|^p} = \rho, \]

The constant \(\rho\) is the rate of convergence. If \(p = 1\), the sequence is said to have linear convergence. If \(p = 2\), the sequence is said to have quadratic convergence.

If the limit does not exist while the upper bound exists for sufficiently large \(n\):

\[ \frac{|x_{n+1} - x^{\ast}|}{|x_n - x^{\ast}|^p} \le \rho, \]

then the order of convergence is at least \(p\) and the rate of convergence is at most \(\rho\).

Theorem 4

The bisection method has a linear convergence rate.

Proof. Without loss of generality, we may assume the initial bracket is on \([0, 1]\). Let the root \(x^{\ast} = 0.b_1 b_2\cdots\) be the binary representation, then the sequence of bisection method can be written as

\[ x_n = 0.b_1 b_2\cdots b_{n-1} 1, \]

where \(b_i\) is the \(i\)-th bit of the binary representation. The error at the \(n\)-th iteration is

\[\begin{split}|x_n - x^{\ast}| = |2^{-n} - \sum_{j\ge n} 2^{-j} b_j|=\begin{cases}\sum_{j>n} 2^{-j}b_j & \text{if } b_n=1\\\sum_{j>n} 2^{-j}(1-b_j) &\text{if }b_n=0\end{cases}\end{split}\]

For each \(n\) such that \(b_n = 1\) (similar argument holds for \(b_n=0\)), we can find \(s\in\N\) such that \(b_{s+n} = 1\), otherwise we arrive at the exact solution.

\[ \frac{|x_{n+k} - x^{\ast}|}{|x_n - x^{\ast}|} = \frac{\sum_{j>n+k} 2^{-j} b_j }{\sum_{j>n} 2^{-j}b_j} \le \frac{2^{-(n+k)}}{2^{-(n+s)}} =2^{s - k} \]

Therefore, the geometric mean of the convergence rate is bounded by \(\frac{1}{2}\).

\[ \rho =\lim_{k\to\infty} \sqrt[k]{\frac{|x_{n+k} - x^{\ast}|}{|x_n - x^{\ast}|}} \le \lim_{k\to\infty}2^{s/k}\frac{1}{2} =\frac{1}{2}. \]

A common technique to study the order of convergence is to use the Taylor expansion. Let us use the Illinois method as an example.

Theorem 5

The order of convergence of the Illinois method is \(\sqrt[3]{3}\).

Big \(O\) Notation

The big \(O\) notation is used to describe the upper bound of a function. For instance, \(f(x) = \cO(g(x))\) means that there exists a constant \(C\) such that \(|f(x)| \le C|g(x)|\) for sufficiently large \(x\).

Another notation is the big \(\Theta\) notation, which is used to describe the upper and lower bounds of a function. For instance, \(f(x) = \Theta(g(x))\) means that there exists constants \(C_1\) and \(C_2\) such that \(C_1|g(x)| \le |f(x)| \le C_2|g(x)|\) for sufficiently large \(x\).

Selection of Decay Factor

In theory, any decay factor smaller than one would produce the same order of convergence. The decay factor \(\frac{1}{2}\) is optimal in the sense that any smaller factor will make the adjustment iteration move away from the root. While having a larger factor will not be efficient to reduce the bracket size.

Proof. Let \(f(x)\in C^2[a, b]\) and \(x^{\ast}\) be the root. For simplicity, we assume the root is simple. Let \(\ell(x)\) be the linear interpolation of \(f\) at the bracket \([a, b]\), the slope is \(f'(\zeta)=\frac{f(b) - f(a)}{b -a}\), then

\[ f(x) - \ell(x) = \frac{f''(\xi)}{2}(x - a)(x - b), \]

therefore, the root satisfies \(-\ell(x^{\ast}) = \frac{f''(\xi)}{2}(x^{\ast} - a)(x^{\ast} - b)\). Let \(c\) be the next iteration, then

\[ -\ell(x^{\ast}) + \ell(c) = -(x^{\ast} - c) f'(\zeta) = \frac{f''(\xi)}{2}(x^{\ast} - a)(x^{\ast} - b), \]

which implies \((x^{\ast} - c) = -\frac{f''(\xi)}{2 f'(\zeta)} (x^{\ast} - a)(x^{\ast} - b)\). When \(a\) and \(b\) are already near the root, the right-hand side has a fixed sign, which implies that using the unadjusted false position method, the new point \(c\) will always fall into a fixed side of the root.

Now, we consider the adjustment of the Illinois method. Without loss of generality, we assume that \(c\) falls into the left side of the root. Then the next iteration will be using the bracket \([c, b]\) with adjusted weight for \(b\). The new point \(c'\) satisfies

\[\begin{split} \begin{aligned} c' - x^{\ast} &= \frac{\frac{1}{2} f(b) (c - x^{\ast}) - f(c) (b - x^{\ast})}{\frac{1}{2} f(b) - f(c)} = (c - x^{\ast}) - \frac{f(c)(b-x^{\ast})}{\frac{1}{2}f(b) - f(c)} + \frac{f(c)(c - x^{\ast})}{\frac{1}{2}f(b) - f(c)}\\ &= (c - x^{\ast}) - \frac{2f(c)}{f(b)}(b - x^{\ast}) \frac{1}{1 - \frac{2f(c)}{f(b)}} + \frac{2f(c)}{f(b)}(c - x^{\ast}) \frac{1}{1 - \frac{2f(c)}{f(b)}} \end{aligned} \end{split}\]

Using Taylor expansion at \(x^{\ast}\), since \(f'\) is bounded from below near the root,

\[\begin{split} \begin{aligned} \frac{f(c)}{f(b)} &= \frac{(c - x^{\ast}) + \cO(|c - x^{\ast}|^2)}{(b - x^{\ast}) + \cO(|b - x^{\ast}|^2)} = \frac{(c - x^{\ast})}{(b - x^{\ast})}\left(1 + \frac{\cO(|b-x^{\ast}|) + \cO(|c - x^{\ast}|)}{1 + \cO(|b-x^{\ast}|)}\right) \\&= \frac{(c - x^{\ast})}{(b - x^{\ast})}\left(1 + \cO(|b-x^{\ast}|)\right) = \frac{f''(\xi)}{2 f'(\zeta)} (a - x^{\ast}) \left(1 + \cO(|b-x^{\ast}|)\right) . \end{aligned} \end{split}\]

Therefore, we obtain

\[\begin{split} \begin{aligned} c' - x^{\ast} &= (c- x^{\ast}) + 2(x^{\ast} - c)\left(1 + \cO(|a - x^{\ast}|)\right) + \cO(|a - x^{\ast}||c - x^{\ast}|) \\ &= (x^{\ast} - c)(1 + \cO(|a - x^{\ast}|)), \end{aligned} \end{split}\]

which moves the new point to the other side of the root by almost reflection. The next bracket becomes \([c, c']\) such that \(|c - x^{\ast}|\approx |c' - x^{\ast}|\). Therefore, if denote the last two points as \(c = x_{n}\) and \(c'=x_{n+1}\), then \(|x_{n} - x^{\ast}|\approx |x_{n+1} - x^{\ast}|\), currently \(c'\) is on the right side of the root. The next three iterations satisfy

• \(|x_{n+2} - x^{\ast}| = \Theta(|x_n - x^{\ast}| |x_{n+1} - x^{\ast}|) = \Theta(|x_{n+1} - x^{\ast}|^2) \), the iteration \(x_{n+2}\) is on left side.

• \(|x_{n+3} - x^{\ast}| = \Theta(|x_{n+2} - x^{\ast}||x_{n+1} - x^{\ast}|) = \Theta(|x_{n+1} - x^{\ast}|^3)\), the iteration \(x_{n+3}\) is on left side, now two consecutive iterations on left side

• \(|x_{n+4} - x^{\ast}| = \Theta(|x_{n+3} - x^{\ast}|) = \Theta(|x_{n+1} - x^{\ast}|^3)\), adjusted by Illinois method, the iteration \(x_{n+4}\) is on right side, which completes a cycle.

By the previous definition of order of convergence, the Illinois method has an order of convergence \(\sqrt[3]{3}\).

Remark 3

A quick way of proof uses the asymptotic analysis.

If we denote \(\epsilon_{n} =x_n - x^{\ast}\), then if \(\epsilon_{i} > 0\) and \(\epsilon_{i-1} < 0\), then \(\epsilon_{i+1} \simeq C \epsilon_{i}\epsilon_{i-1}\), where \(C\approx \frac{f''(x^{\ast})}{2f'(x^{\ast})}\), here we may assume \(C > 0\), then \(\epsilon_{i+1} < 0\).

The next iteration applies the same rule, which is \(\epsilon_{i+2} \simeq C \epsilon_{i+1}\epsilon_{i} = C^2 \epsilon_i^2 \epsilon_{i-1} < 0\).

Now, we need an adjustment step, which gives \(\epsilon_{i+3} = -C^2 \epsilon_i^2 \epsilon_{i-1} > 0\), which completes a cycle. Every three iterations make a cycle, to derive the order of convergence, we continue to derive \(\epsilon_{i+6}\), which equals to \(-C^8\epsilon_{i}^6\epsilon_{i-1}^3 = C^2\epsilon_{i+3}^3\).

Therefore, the order of convergence is \(\sqrt[3]{3}\).

By taking the decay factor as a constant, the order of convergence is at most \(\sqrt[3]{3}\). However, if the decay factor can be chosen as a variable one, the order of convergence can be improved. Next, we apply the same technique to analyze the Pegasus method, which is a variant of the Illinois method. It is interesting that the Pegasus method was initially discovered in a subroutine for the Ferranti Pegasus computer, but no author information is included. The convergence analysis was given by [Dowell and Jarratt, 1972] after its discovery.

Theorem 6 (Pegasus Method)

The Pegasus method is a variant of Illinois method, whose updating scheme for the decay factors is replaced by

\[\lambda=\frac{f_1}{f_1+f_2}\]

in Algorithm 2. Then the order of convergence is \(\sqrt[4]{\frac{7+\sqrt{57}}{2}}\).

Proof. Similar to the previous Remark 3, we assume \(\epsilon_{i-1} < 0\) and \(\epsilon_{i} > 0\) such that \(|\epsilon_i|\ll |\epsilon_{i-1}|\), the constant \(C = \frac{f''}{2f'}|_{x=x^{\ast}} > 0\). Similar to the previous analysis, we can derive more terms in the asymptotic form

\[ \epsilon_{i+1}\simeq C \epsilon_{i}\epsilon_{i-1} + D \epsilon_{i} \epsilon_{i-1}(\epsilon_{i} + \epsilon_{i-1})< 0 \]

where \(D = -C^2 + \frac{f'''}{6f'}|_{x=x^{\ast}}\). Then due to a different sign for \(\epsilon_{i}\) and \(\epsilon_{i+1}\), we can derive

\[\epsilon_{i+2}\simeq C \epsilon_{i+1}\epsilon_{i} + D \epsilon_{i+1} \epsilon_{i}(\epsilon_{i+1} + \epsilon_{i}) < 0.\]

Now use the adjustment step, we obtain (needs some calculation)

\[\begin{split} \begin{aligned} \epsilon_{i+3} &= \frac{\epsilon_{i+2}\lambda f(x_i) - \epsilon_i f(x_{i+2})}{\lambda f(x_{i}) - f(x_{i+2})},\quad \lambda = \frac{f(x_{i+1})}{f(x_{i+2}) + f(x_{i+1})}\\ &\approx C^2 \epsilon_{i}^2\epsilon_{i+2} - D \epsilon_{i} \epsilon_{i+1} \epsilon_{i+2}, \end{aligned} \end{split}\]

it implies \(\epsilon_{i+3} < 0\) as well. Therefore, another adjustment step is needed and (after some more calculations, see Remark 4)

\[\begin{split} \begin{aligned} \epsilon_{i+4} &= \frac{\epsilon_{i+3}\lambda f(x_i) - \epsilon_i f(x_{i+3})}{\lambda f(x_{i}) - f(x_{i+3})},\quad \lambda = \frac{f(x_{i+2})}{f(x_{i+3}) + f(x_{i+2})}\frac{f(x_{i+1})}{f(x_{i+2}) + f(x_{i+1})}\\ &\approx C^5 \epsilon_{i}^4\epsilon_{i+1}^2 > 0. \end{aligned} \end{split}\]

Therefore, \(\epsilon_{i+4}\approx C \epsilon_{i+3} \epsilon_{i+2}\). A full cycle consists of 4 iterations and

\[ \epsilon_{i+8} \simeq C^7 \epsilon_{i+4}^6 \epsilon_{i+3}^2 \simeq C^{8}\epsilon_{i+4}^7 \epsilon_i^2. \]

The order of convergence \(p\) solves \(p^2 - 7p - 2 = 0\), which gives \(p = \sqrt[4]{\frac{7+\sqrt{57}}{2}}\).

Remark 4

Actually, the asymptotic expansion of \(\epsilon_{i+4}\) is (expanded in \(\epsilon_{i+1}\) first )

\[ \epsilon_{i+4} = C^5 \epsilon_{i}^4\epsilon_{i+1}^2 + C^7 \epsilon_i^7\epsilon_{i+1} + \cO(\epsilon_{i}^8 \epsilon_{i+1} + \epsilon_{i}^5 \epsilon_{i+1}^2 + \epsilon_{i}^3 \epsilon_{i+1}^3)\]

It is at first not clear why we can retain the first term and drop the rest, because it requires the following inequalities to hold

\[\epsilon_{i}^3 \ll \epsilon_{i+1}\ll \epsilon_i.\]

The latter one is correct because of the relation \(\epsilon_{i+1} \simeq C \epsilon_{i}\epsilon_{i-1} \ll \epsilon_i\) once the iterations are close to the root. The former one is equivalent to \(\epsilon_{i}^2 \ll \epsilon_{i-1}\). This can be made into an assumption because \(\epsilon_{i-1}\) and \(\epsilon_{i}\) are the last two iterations in the previous cycle, which correspond to \(\epsilon_{i+3}\) and \(\epsilon_{i+4}\) in the current cycle, we can use the above asymptotic estimate to find \(\epsilon_{i+3}^3 \approx C^6 \epsilon_{i}^9 \epsilon_{i+1}^3 \ll \epsilon_{i+4}\).

Therefore, if the inequality does not hold, one can use the current cycle as the starting point to perform the same analysis.

In the Pegasus method, two consecutive standard false position steps and two adjustment steps are performed in a full cycle. Although the decay over a full cycle is significant, such advantage will be gone if the cycle is long.

In order to make the cycle shorter, we need to drop at least one step. According to the previous analysis (see Illinois method), the standard false position step does not change the sign of the error asymptotically, thus it is preferred to drop one of the standard false position steps.

Let us finish this section with a brief discussion on the improved Pegasus method [King, 1973], which takes the advantage of the symmetry in the false position method to avoid consecutive standard false position steps. The similar technique can be also applied to other methods such as Anderson-Bjorck method [Anderson and Björck, 1973].

The Improved Pegasus method can also use the decay factor of Anderson-Bjorck method as

\[\begin{split}\lambda = \begin{cases} 1 - \frac{f_2}{f_1}&\text{ if positive},\\ \frac{1}{2}&\text{ otherwise}.\end{cases} \end{split}\]

Algorithm 3 (Improved Pegasus Method)

Inputs \(f\), \(a, b\), \(\epsilon\)

Outputs \(c\)

1. \(x_0\gets a\), \(x_1\gets b\), \(f_0\gets f(x_0)\), \(f_1\gets f(x_1)\), \(iter = 0\). //initialization.

2. While True do:

   1. \(x_2 \gets x_0 - f_0 \frac{x_1 - x_0}{f_1 - f_0}\), \(f_2 \gets f(x_2)\), \(iter \gets iter + 1\). //standard false position step, \(x_1\) and \(x_2\) are two latest iterations.

   2. If \( |f_2| < \epsilon\), then return \(x_2\); //check stopping criteria

   3. If \(f_1 f_2 < 0\), then swap \((x_0, f_0)\) and \((x_1, f_1)\). // avoid false position step

   4. While \(f_1 f_2 > 0\), then // adjust until the sign changes

      • \((x_0, f_0)\gets (x_0, \lambda f_0)\), where \(\lambda =\frac{f_1}{f_1 + f_2}\).

      • \((x_1, f_1)\gets (x_2, f_2)\) and \(x_2 \gets x_0 - f_0 \frac{x_1 - x_0}{f_1 - f_0}\), \(f_2 \gets f(x_2)\), \(iter\gets iter + 1\)

      • If \( |f_2| < \epsilon\), then return \(x_2\); //check stopping criteria

   5. If \(f_1 f_2 < 0\), then // perform a false position step

      • \((x_0, f_0)\gets (x_1, f_1)\) and \((x_1, f_1)\gets (x_2, f_2)\).

Theorem 7

The improved Pegasus method has an order of convergence at least \(\sqrt[3]{5}\).

Proof. In the same setting as the previous theorem, we will find the first iteration is the same as the Pegasus method (false position) that

\[\epsilon_{i+1}\simeq C \epsilon_i \epsilon_{i-1} + D \epsilon_i \epsilon_{i-1}(\epsilon_{i} + \epsilon_{i-1}) < 0.\]

For the second iteration, although \(x_{i+1}\) and \(x_{i}\) are on different sides, but a false position step has been performed in the previous step, thus it will perform adjustment step with \(\lambda = \frac{f(x_{i}) - f(x_{i+1})}{f(x_{i})}\), then

\[\epsilon_{i+2} \simeq -D \epsilon_{i+1}\epsilon_{i}\epsilon_{i-1},\]

note the leading term is different from the Pegasus method. There are two options.

• If \(D < 0\), then \(\epsilon_{i+2} > 0\) which completes a cycle with two iterations, which is very compact. In this case, we find

  \[\epsilon_{i+4} \simeq C^{2} \epsilon_{i+2}^3,\]

  which implies an order of convergence at \(\sqrt{3}\).

• If \(D > 0\), then it will perform another adjustment step,

  \[\epsilon_{i+3} = \frac{\epsilon_{i+2}\lambda f(x_i) - \epsilon_i f(x_{i+2})}{\lambda f(x_{i}) - f(x_{i+2})},\quad \lambda = \frac{f(x_{i+1})-f(x_{i+2})}{f(x_{i+1})}\frac{f(x_{i})-f(x_{i+1}) }{ f(x_{i})}\]

  which will make \(\epsilon_{i+3} \simeq -C^3 D \epsilon_{i-1}^3 \epsilon_i^3 > 0\), which completes a cycle with three iterations. In this case, we find

  \[\epsilon_{i+6} \simeq D^{2} \epsilon_{i+3}^5,\]

  which implies an order of convergence at \(\sqrt[3]{5}\).

Iterative Methods

In theory, the bracket methods are also iterative. However, one important phenomenon of bracket methods is that the current iteration may not be generated by the latest iterations. For instance, the false position method will have a stalled endpoint once the iterations are near the root.

In the following, we will introduce some iterative methods, which generate new iterations based on the latest iterations.

Newton-Raphson Method

Let \(f\in C^2\), if the current iteration \(x_{n}\) is close to a root, then the Newton-Raphson method computes the next iteration \(x_{n+1}\) by

\[ x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}. \]

Geometrically, \(x_{n+1}\) is the \(x\)-intercept of the tangent line of \(f\) at \((x_n, f(x_n))\). Asymptotically, the Newton-Raphson method has a quadratic convergence rate if the root is simple.

Fig. 3 Newton-Raphson method

Limitation

The analysis can only provide convergence for the initial guess close to the root. If the initial guess is far from the root, then it is possible that the sequence diverges.

Theorem 8

Suppose \(f\in C^2\) and the root is simple, then the order of convergence for Newton-Raphson method is 2.

Proof. The Taylor expansion of \(f\) at \(x_n\) is

\[f(x) = f(x_n) + f'(x_n)(x - x_n) + \frac{f''(\zeta)}{2}(x - x_n)^2\]

where \(\zeta\in (x, x_n)\). The root is \(x^{\ast}\), then \(f(x^{\ast}) = 0\), thus

\[0 = f(x_n) + f'(x_n)(x^{\ast} - x_n) + \frac{f''(\zeta)}{2}(x^{\ast} - x_n)^2.\]

Therefore,

\[ -\frac{f''(\zeta)}{f'(x_n)} (x^{\ast} - x_n)^2= x^{\ast} - x_n + \frac{f(x_n)}{f'(x_n)} = x^{\ast} - x_{n+1}. \]

If we denote the error \(\epsilon_{n} = x_n - x^{\ast}\), then \(\epsilon_{n+1} = \frac{f''(\zeta)}{f'(x_n)} \epsilon_{n}^2\), which implies the order of convergence is 2.

Secant Method

The secant method is similar to the false position method, but it always uses the latest two iterations to perform the next iteration. The secant method computes the next iteration \(x_{n+1}\) by

\[ x_{n+1} = x_n - f(x_n)\frac{x_n - x_{n-1}}{f(x_n) - f(x_{n-1})}. \]

Using the same derivation for the Illinois method, the error \(\epsilon_{n} = x_n - x^{\ast}\) satisfies \(\epsilon_{n+1} \simeq C \epsilon_{n}\epsilon_{n-1}\), see Fig. 4.

Fig. 4 Secant method

The order of convergence \(p\) solves the equation \(p^2 - p - 1 = 0\), which gives \(p = \frac{1 + \sqrt{5}}{2}\).

Applications in Optimization

Exercises

Theoretical Part

Problem 1

Briefly explain why the formula \(x_2 = x_0 - f_0 \frac{x_1 - x_0}{f_1 - f_0}\) in bracket methods will not suffer from catastrophic cancellation.

Problem 2

The only difference between Algorithm 3 and the usual Pegasus method is at the 3rd step in the while loop, which eliminates consecutive false position steps. This change is very simple, but it leads to an improvement in the order of convergence.

Explain why the Illinois method \(\lambda = \frac{1}{2}\) cannot be faster by the same technique.

Problem 3

Suppose \(f\in C^2\) and the root \(x^{\ast}\) has a multiplicity of \(m>1\), then the order of convergence for the Newton-Raphson method is \(1\). The modified Newton-Raphson method

\[ x_{n+1} = x_n - m\frac{f(x_n)}{f'(x_n)} \]

has an order of convergence \(2\).

Computational Part

Extended Reading

[ABjorck73]

Ned Anderson and Åke Björck. A new high order method of regula falsi type for computing a root of an equation. BIT Numerical Mathematics, 13(3):253–264, 1973.

[DJ72]

M Dowell and P Jarratt. The “pegasus” method for computing the root of an equation. BIT Numerical Mathematics, 12:503–508, 1972.

[Kin73]

Richard F King. An improved pegasus method for root finding. BIT Numerical Mathematics, 13:423–427, 1973.