# 10.2: Appendix B - Iteration Methods

- Page ID
- 44514

## 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 method^{1} that is illustrated in Figure \(\PageIndex{1}\). This 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

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

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

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.

## 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

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

\[\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 chaos^{4} ^{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}\]

This leads to the condition

\[\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}\]

which leads to

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

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