# 10.2: Appendix B - Iteration Methods

• • R.L. Cerro, B. G. Higgins, S Whitaker
• Professors (Chemical Engineering) at University of Alabama at Huntsville & University of California at Davis

## B1: Bisection method

Given some function of $$x$$ such as $$H(x)$$, we are interested in the solution of the equation

$H(x) = 0 , \quad x = x^{*} \label{1}$

Here we have used $$x^{*}$$ to represent the solution. For simple functions such as $$H(x) = x-b$$ we obtain a single solution given by $$x^{*} = b$$, while for a more complex function such as $$H(x) = x^{2} -b$$ we obtain more that one solution as indicated by $$x^{*} = \pm \sqrt{b}$$. In many cases there is no explicit solution for Equation \ref{1}. For example, if $$H(x)$$ is given by

$H(x) = a\sin ({x\pi / 2} ) + b\cos (2\pi x)\label{2}$

we need to use iterative methods to determine the solution $$x = x^{*}$$.

The simplest iterative method is the bisection method1 that is illustrated in Figure $$\PageIndex{1}$$. This method Figure $$\PageIndex{1}$$: Illustration of the bisection method

begins by locating $$x_{o}$$ and $$x_{1}$$ such that $$H(x_{o} )$$ and $$H(x_{1} )$$ have different signs. In Figure $$\PageIndex{1}$$ we see that $$x_{o}$$ and $$x_{1}$$ have been chosen so that there is a change of sign for $$H(x)$$, i.e.,

$H(x_{o} )>0 , \quad H(x_{1} )<0 \label{3}$

Thus if $$H(x)$$ is a continuous function we know that a solution $$H(x^{*} ) = 0$$ exists somewhere between $$x_{o}$$ and $$x_{1}$$. We attempt to locate that solution by means of a guess (i.e., the bisection) indicated by

$x_{2} = \frac{x_{o} + x_{1} }{ 2} \label{4}$

As illustrated in Figure $$\PageIndex{1}$$, this guess is closer to the solution, $$x = x^{*}$$, than either $$x_{o}$$ or $$x_{1}$$, and if we repeat this procedure we will eventually find a value of $$x$$ that produces a value of $$H(x)$$ that is arbitrarily close to zero. In terms of the particular graph illustrated in Figure $$\PageIndex{1}$$, it is clear that $$x_{3}$$ will be located between $$x_{1}$$ and $$x_{2}$$; however, this need not be the case. For example, in Figure $$\PageIndex{2}$$ we have represented a slightly different function for which $$x_{3}$$ will be located between $$x_{o}$$ and $$x_{2}$$. The location of the next guess is based on the idea that the function $$H(x)$$ must change sign. In order to determine the location of the Figure $$\PageIndex{2}$$: Alternate choice for the second bi-section

next guess we examine the products $$H(x_{n} )\times H(x_{n-1} )$$ and $$H(x_{n} )\times H(x_{n-2} )$$ in order to make the following decisions:

\begin{align} H(x_{n} )\times H(x_{n-1} )<0 , \quad x_{n+1} = \frac{x_{n} + x_{n-1} }{ 2} \nonumber \\ H(x_{n} )\times H(x_{n-2} )<0 , \quad x_{n+1} = \frac{x_{n} + x_{n-2} }{ 2} \label{5} \end{align}

Since these two choices are mutually exclusive there is no confusion about the next choice of the dependent variable. The use of Eqs. \ref{5} is crucial when the details of $$H(x)$$ are not clear, and a program is written to solve the implicit equation.

## B2: False position method

The false position method is also known as the method of interpolation 2 and it represents a minor variation of the bisection method. Instead of bisecting the distance between $$x_{o}$$ and $$x_{1}$$ in Figure $$\PageIndex{1}$$ in order to locate the point $$x_{2}$$, we use the straight line indicated in Figure $$\PageIndex{3}$$. Sometimes this line is called Figure $$\PageIndex{3}$$: False position construction

the secant line. The definition of the tangent of the angle $$\varphi$$ provides

$\tan \varphi = \frac{0 - H(x_{1} )}{ x_{2} - x_{1} } = \frac{H(x_{o} ) - H(x_{1} )}{ x_{o} - x_{1} }\label{6}$

and we can solve for $$x_{2}$$ to obtain

$x_{2} = x_{1} - \frac{\left(x_{o} - x_{1} \right)H(x_{1} )}{ H(x_{o} ) - H(x_{1} )} \label{7}$

This replaces Equation \ref{3} in the bisection method and it can be generalized to obtain

$x_{n+2} = x_{n+1} - \frac{\left(x_{n} - x_{n+1} \right)H(x_{n+1} )}{ H(x_{n} ) - H(x_{n+1} )} \label{8}$

Application of successive iterations will lead to a value of $$x$$ that approaches $$x^{*}$$ shown in Figure $$\PageIndex{3}$$.

## B3: Newton’s method

Newton’s method (also known as the Newton-Raphson method)3 is named for Sir Isaac Newton and is perhaps the best known method for finding roots of real valued functions. The method is similar to the false position method in that a straight line is used to locate the next estimate of the root of an equation; however, in this case it is a tangent line and not a secant line. This is illustrated in Figure $$\PageIndex{4a}$$ where we Figure $$\PageIndex{4a}$$: First estimate using Newton’s method

have chosen $$x_{o}$$ as our first estimate of the solution to Equation \ref{1} and we have constructed a tangent line to $$H(x)$$ at $$x_{o}$$. The slope of this tangent line is given by

$\left. \frac{dH}{ dx} \right|_{x = x_{o} } = \frac{H(x_{o} ) - 0}{ x_{o} - x_{1} } \label{9}$

and we can solve this equation to produce our next estimate of the root. This new estimate is given by

$x_{1} = x_{o} - \frac{H(x_{o} )}{ \left. {(dH / dx)} \right|_{x = x_{o} } } \label{10}$

and we use this result to determine $$H(x_{1} )$$ as indicated in Figure $$\PageIndex{4a}$$. Given $$H(x_{1} )$$ and $$x_{1}$$ we can construct a second estimate as indicated in Figure $$\PageIndex{4b}$$, and this process can be continued to find the solution given by $$x^{*}$$. The general iterative procedure is indicated by

$x_{n+1} = x_{n} - \frac{H(x_{n} )}{ \left. {(dH / dx)} \right|_{x = x_{n} } } , \quad n = 0,1,2,..., \label{11}$

Newton’s method is certainly an attractive technique for finding solutions to implicit equations; however, it does require that one know both the function and its derivative. For complex functions, calculating the derivative at each step in the iteration may require more effort than that associated with the bisection method or the false position method. In addition, if the derivative of the function is zero in the region of interest, Newton’s method will fail. Figure $$\PageIndex{4b}$$: Second estimate using Newton’s method

## B4: Picard’s method

Picard’s method for solving Equation \ref{1} begins by defining a new function according to

Definition: $f(x) = x + H(x) \label{12}$

Given any value of the dependent variable, $$x_{n}$$, we define a new value, $$x_{n+1}$$, by

Definition: $x_{n+1} = f(x_{n} ) , \quad n = 0,1,2,3,... \label{13}$

This represents Picard’s method or the method of direct substitution or the method of successive substitution. If this procedure converges, we have

$f(x^{*} ) = x^{*} + H(x^{*} ) = x^{*} \label{14}$

In Equation \ref{13} we note that the function $$f(x_{n} )$$, maps the point $$x_{n}$$ to the new point $$x_{n+1}$$. If the function $$f(x)$$ maps the point $$x$$ to itself, i.e., $$f(x) = x$$, then $$x$$ is called the fixed point of $$f(x)$$. In Figure $$\PageIndex{5}$$ we again consider the function represented in Figures B-1, B-3 and B-4, and we illustrate the functions $$f(x)$$, $$y(x)$$ and $$H(x)$$. The graphical representation of the fixed point, $$x^{*}$$, is the intersection of the Figure $$\PageIndex{5}$$: Picard’s method

function of $$f(x)$$ with the line $$y = x$$. Note that not all functions have fixed points. For example if $$f(x)$$ is parallel to the line $$y = x$$ there can be no intersection and no fixed point. Given our first estimate, $$x_{o}$$, we use Equation \ref{13} to compute $$x_{1}$$ according to

$x_{1} = f(x_{o} ) \label{15}$

Clearly $$x_{1}$$ is further from the solution, $$x^{*}$$, than $$x_{o}$$ and we can see from the graphical representation in Figure $$\PageIndex{5}$$ that Picard’s method diverges for this case. If $$x_{o}$$ were chosen to be less than the solution, $$x^{*}$$, we would also find that the iterative procedure diverges. If the slope of $$f(x)$$ were less than the slope of $$y(x)$$, we would find that Picard’s method converges. This suggests that the method is useful for “weak” functions of $$x$$, i.e., $$\left|{df / dx} \right|<1$$ and this is confirmed in Sec. B6.

## B5: Wegstein’s method

In Figure $$\PageIndex{6}$$ we have illustrated the same function, $$f(x)$$, that appears in Figure $$\PageIndex{5}$$. For some point $$x_{ 1}$$ in the neighborhood of $$x_{o}$$ we can approximate the derivative of $$f(x)$$ according to Figure $$\PageIndex{6}$$: Wegstein’s method

$\frac{df}{ dx} \approx \frac{f(x_{1} ) - f(x_{o} )}{ x_{1} - x_{o} } \approx {slope} = S \label{16}$

and we can use this result to obtain an approximation for the function $$f(x_{1} )$$.

$f(x_{1} )\approx f(x_{o} ) + S\left(x_{1} - x_{o} \right) \label{17}$

At this point we recall Equation \ref{14} in the form

$f(x^{*} ) = x^{*} \label{18}$

and note that if $$x_{ 1}$$ is in the neighborhood of $$x^{*}$$ we obtain the approximation

$f(x_{1} )\approx x_{1} \label{19}$

We use this result in Equation \ref{17} to produce an equation

$x_{1} = f(x_{o} ) + S\left(x_{1} - x_{o} \right) \label{20}$

in which $$S$$ is an adjustable parameter that is used to determine the next step in the iterative procedure. It is traditional, but not necessary, to define a new adjustable parameter according to

$q = \frac{S}{ S-1} \label{21}$

Use of this representation in Equation \ref{20} leads to

$x_{1} = (1-q)f(x_{o} ) + q x_{o} \label{22}$

and we can generalize this result to Wegstein’s method given by

$x_{n+1} = (1-q)f(x_{n} ) + q x_{n} , \quad n = 0, 1, 2, 3, ...{etc} \label{23}$

When the adjustable parameter is equal to zero, $$q = 0$$, we obtain Picard’s method described in Sec. B4. When the adjustable parameter greater than zero and less than one, $$0<q<1$$, we obtain a damped successive substitution process that improves stability for nonlinear systems. When the adjustable parameter is negative, $$q<0$$, we obtain an accelerated successive substitution that may lead to an unstable procedure.

## B6: Stability of iteration methods

In this section we consider the linear stability characteristics of Newton’s method, Picard’s method, and Wegstein’s method that have been used to solve the implicit equation given by

$H(x) = 0 , \quad x = x^{*} \label{24}$

The constraint associated with the linear analysis will be listed below and it must be kept in mind when interpreting results such as those presented in Chapter 7. We begin by recalling the three iterative methods as

Newton’s method: $x_{n+1} = x_{n} - \frac{H(x_{n} )}{ \left. {(dH / dx)} \right|_{x = x_{n} } } , \quad n = 0,1,2,..., \label{25}$

Picard’s method: $x_{n+1} = f(x_{n} ) , \quad n = 0,1,2,.... \label{26}$

Wegstein’s method: $x_{n+1} = (1-q)f(x_{n} ) + q x_{n} , \quad n = 0, 1, 2, ... \label{27}$

in which the auxiliary function, $$f(x)$$, is defined by

Definition: $f(x) = x + H(x) \label{28}$

The general form of these three iterative methods is given by

$x_{n+1} = G(x_{n} ) , \quad n = 0,1,2,..... \label{29}$

and for each of the three methods on seeks to find the fixed point $$x^{*}$$ of $$G(x)$$ such that

$x^{*} = G(x^{*} ). \label{30}$

Our stability analysis of Eqs. \ref{25} through \ref{27} is based on linearizing $$G(x)$$ about the fixed point $$x^{*}$$. We let $$\Delta x_{n}$$ and $$\Delta x_{n+1}$$ be small perturbations from the fixed point as indicated by

$x_{n} = x^{*} + \Delta x_{n} , \quad x_{n+1} = x^{*} + \Delta x_{n+1} \label{31}$

This allows us to express Equation \ref{29} as

$x^{*} + \Delta x_{n+1} = G(x^{*} +\Delta x_{n} ) , \quad n = 0,1,2,..... \label{32}$

and a Taylor series expansion (See Problems 5-30 and 5-31) leads to

$G(x^{*} +\Delta x_{n} ) = G(x^{*} ) + \Delta x_{n} \left(\frac{dG}{ dx} \right)_{x^{*} } + { \frac{1}{ 2}} \Delta x_{n}^{2} \left(\frac{d^{2} G}{ dx^{2} } \right)_{x^{*} } + ..... \label{33}$

On the basis of Equation \ref{30} this infinite series simplifies to

$G(x^{*} +\Delta x_{n} ) = x^{*} + \Delta x_{n} \left(\frac{dG}{ dx} \right)_{x^{*} } + { \frac{1}{ 2}} \Delta x_{n}^{2} \left(\frac{d^{2} G}{ dx^{2} } \right)_{x^{*} } + ..... \label{34}$

and we can use Equation \ref{32} to represent the left hand side in a simpler form to obtain

$x^{*} + \Delta x_{n+1} = x^{*} + \Delta x_{n} \left(\frac{dG}{ dx} \right)_{x^{*} } + { \frac{1}{ 2}} \Delta x_{n}^{2} \left(\frac{d^{2} G}{ dx^{2} } \right)_{x^{*} } + ..... \label{35}$

At this point we impose a constraint on the higher order terms expressed as

Constraint: $\Delta x_{n} \left(\frac{dG}{ dx} \right)_{x^{*} } >>{ \frac{1}{ 2}} \Delta x_{n}^{2} \left(\frac{d^{2} G}{ dx^{2} } \right)_{x^{*} } + ..... \label{36}$

so that Equation \ref{35} takes the form

$\Delta x_{n+1} = \Delta x_{n} \left(\frac{dG}{ dx} \right)_{x^{*} } , \quad n = 0,1,2,3,4,..... \label{37}$

If we write a few of these equations explicitly as

$\Delta x_{1} = \Delta x_{o} \left(\frac{dG}{ dx} \right)_{x^{*} } \label{38a}$

$\Delta x_{2} = \Delta x_{1} \left(\frac{dG}{ dx} \right)_{x^{*} }\label{38b}$

$\Delta x_{3} = \Delta x_{2} \left(\frac{dG}{ dx} \right)_{x^{*} }\label{38c}$

........................

$\Delta x_{n} = \Delta x_{n-1} \left(\frac{dG}{ dx} \right)_{x^{*} }\label{38d}$

it becomes clear that they can be used to provide a general representation given by

$\Delta x_{n} = \Delta x_{o} \left[\left({dG / dx} \right)_{x^{*} } \right]^{n} , \quad n = 0,1,2,..... \label{39}$

At this point we see that $$\Delta x_{n} \to 0$$ when $$n\to \infty$$ provided that

$\left|\left({dG / dx} \right)_{x^{*} } \right|<1 \label{40}$

When $$\Delta x_{n} \to 0$$ as $$n\to \infty$$ the system converges and one says that the fixed point $$x^{*}$$ is attracting. The three special cases represented by Equation \ref{39} can be expressed as

I. $\left|\left({dG / dx} \right)_{x^{*} } \right|<1 , \quad \text{the fixed point }x^{*} \text{ is }attracting\text{ and the iteration }converges \label{41}$

II. $\left|\left({dG / dx} \right)_{x^{*} } \right|>1 , \quad \text{the fixed point }x^{*} \text{ is }repelling\text{ and the iteration }diverges \label{42}$

III. $\left|\left({dG / dx} \right)_{x^{*} } \right| = 1 , \quad \text{the fixed point }x^{*} \text{ is neither }attracting\text{ nor }repelling \label{43}$

It is extremely important to note that the stability analysis leading to these three results is based on the linear approximation associated with Equation \ref{36}. In this development, the word attracting is used for a system that converges since $$x_{n}$$ moves toward $$x^{*}$$ as $$n$$ increases, while the word repelling is used for a system that diverges since $$x_{n}$$ moves away from $$x^{*}$$ as $$n$$ increases. The case in which the fixed point is neither attracting nor repelling can lead to chaos4 5.

At this point we are ready to return to Eqs. \ref{25}, \ref{26} and \ref{27} in order to determine the linear stability characteristics of placeCityNewton’s method, Picard’s method, and Wegstein’s method.

### Newton’s method

In this case we have

$G(x) = x - \frac{H(x)}{ \left({dH / dx} \right)} \label{44}$

and the derivative that is required to determine the stability is given by

$\left({dG / dx} \right) = \frac{H(x)}{ \left({dH / dx} \right)^{2} } \frac{d^{2} H}{ dx^{2} } \label{45}$

Evaluation of this derivative at the fixed point where $$H(x^{*} ) = 0$$ leads to

$\left({dG / dx} \right)_{x^{*} } = 0 \label{46}$

This indicates that Newton’s method will converge provided that $$\left({dH / dx} \right) \neq 0$$ and provided that the initial guess, $$x_{o}$$, is close enough to $$x^{*}$$ so that Equation \ref{36} is satisfied. If Equation \ref{36} is not satisfied, the linear stability analysis leading to Eqs. \ref{41} through \ref{43} is not valid.

### Picard’s method

In this case Eqs. \ref{26} and \ref{29} provide $$G(x) = f(x)$$ and

$\left({dG / dx} \right)_{x^{*} } = \left({df / dx} \right)_{x^{*} } \label{47}$

and from Equation\ref{41} we conclude that Picard’s method is stable when

$\left|\left({df / dx} \right)_{x^{*} } \right|<1 \label{48}$

In Example $$7.5.1$$ we have the fixed point iteration given by Equation 17 that can be expressed as

$x_{n+1} = f(x_{n} ) = \left(1 - {C}\right)\left(1 + x_{n} \right) , \quad {C} = 0.30 , \quad n = 0.1,2,3,..... \label{49}$

$\left|\left({df / dx} \right)\right| = \left|\left(1 - {C}\right)\right|<1 \label{50}$

that produces the stable iteration illustrated in Table $$7.5.1a$$. In Example $$7.5.2$$ we find another example of Picard’s method given by Equation 40 that we repeat here as

$x_{n+1} = f(x_{n} ) = x_{n} + \left[37.6190 - \frac{15.1678 x_{n} }{ 0.8571 - x_{n} } \right] , \quad n = 0,1,2,3,... \label{51}$

The solution is given by $$x^{*} = 0.6108$$ and this leads to

$\left({df / dx} \right)_{x^{*} } = -213.1733 \label{52}$

indicating that Picard’s method is unstable for this particular fixed point iteration. This result is consistent with the entries in Table $$7.5.2a$$.

### Wegstein’s method

In this case Eqs. \ref{27} and \ref{29} provide

$G(x) = (1-q)f(x) + q x \label{53}$

$G(x) = (1-q)f(x) + q x \label{54}$

From this we have

$\frac{dG}{ dx} = (1-q)\frac{df}{ dx} + q \label{55}$

and the stability condition given by Equation \ref{41} indicates that Wegstein’s method will converge provided that

$(1-q)\left|\frac{df}{ dx} \right| + q<1 \label{56}$

Here one can see that the adjustable parameter $$q$$ can often be chosen so that this inequality is satisfied and Wegstein’s method will converge as illustrated in Examples $$7.5.1$$ and $$7.5.2$$.

1. Corliss, G. 1977, Which root does the bisection method find?, SIAM Review 19, 325-327.↩
2. Wylie, C.R., Jr. 1951, Advanced Engineering Mathematics, McGraw-Hill Book Co., Inc., New York.↩
3. Ypma, T.J. 1995, Historical development of the Newton-Raphson method, SIAM Review 37, 531-551.↩
4. Gleick, J. 1988 Chaos: Making a New Science, Penguin Books, New York.↩
5. Peitgen, H-O., Jürgens, H., and Saupe, D., 1992, Chaos and Fractals. New Frontiers of Science, Springer-Verlag.↩