Prev
Unit 2 · Nonlinear Equations
Fixed-Point Iteration and Convergence Order
60m100 XP
Next
🔍Nonlinear Equations

Fixed-Point Iteration and Convergence Order

x = g(x), the contraction mapping theorem, linear rate, Aitken and Steffensen acceleration

Any calculator, any day
Take a scientific calculator. Set it to radians. Type 1, and press cos. Press it again. And again. 0.5403, 0.8576, 0.6543, 0.7935, 0.7014, 0.7640, 0.7221… The numbers swing above and below some value, the swings shrinking each time. After about sixty presses the display stops changing: 0.7390851332. Start from 100 instead. Or from −3. The same number appears. Press √ repeatedly on any positive number and you get 1. Press x2x^2 on 0.9 and you get 0; on 1.1, the calculator overflows. Nothing here is a 'method' anybody designed. It is just a function applied to its own output. Yet it solves an equation — cos⁡x=x\cos x = x — that has no closed-form answer. Some of these loops settle and some explode. The difference turns out to explain every iterative method in this course.
A fixed point of a function gg is a value with g(x∗)=x∗g(x^*) = x^*. The calculator was running the fixed-point iteration
xk+1=g(xk)x_{k+1} = g(x_k)
with g=cos⁡g = \cos. If it converges, it converges to a fixed point (take the limit of both sides, assuming gg is continuous). Any root-finding problem can be put in this form: f(x)=0f(x) = 0 is equivalent to x=x−f(x)x = x - f(x), or x=x+c f(x)x = x + c\,f(x) for any c≠0c \ne 0, or many other rearrangements. The question is which rearrangement converges.
💡The contraction mapping theorem
Suppose gg maps an interval II into itself and is a contraction there: ∣g(x)−g(y)∣≤L∣x−y∣|g(x) - g(y)| \le L|x - y| for some L<1L < 1 — for a differentiable gg, ∣g′(x)∣≤L<1|g'(x)| \le L < 1 on II. Then:
  1. gg has exactly one fixed point x∗x^* in II;
  2. the iteration converges to it from every starting point in II;
  3. the error shrinks at least by LL each step: ∣ek+1∣≤L∣ek∣|e_{k+1}| \le L|e_k|.
Stefan Banach proved the general version in 1922. It is the reason the calculator converges: ∣cos⁡′(x)∣=∣sin⁡x∣≤sin⁡1=0.84|\cos'(x)| = |\sin x| \le \sin 1 = 0.84 for x∈[0,1]x \in [0, 1], and cos⁡\cos maps [0,1][0, 1] into itself.
Derivation
Why a contraction must converge
Consecutive steps shrink. Apply the contraction inequality to two consecutive iterates:
∣xk+1−xk∣=∣g(xk)−g(xk−1)∣≤L ∣xk−xk−1∣≤⋯≤Lk ∣x1−x0∣.|x_{k+1} - x_k| = |g(x_k) - g(x_{k-1})| \le L\,|x_k - x_{k-1}| \le \dots \le L^k\,|x_1 - x_0|.
So the iterates settle. The distance from xkx_k to any later iterate is at most the sum of the steps in between, a geometric series:
∣xk+m−xk∣≤(Lk+Lk+1+… ) ∣x1−x0∣=Lk1−L ∣x1−x0∣  ⟶  0.|x_{k+m} - x_k| \le (L^k + L^{k+1} + \dots)\,|x_1 - x_0| = \frac{L^k}{1 - L}\,|x_1 - x_0| \;\longrightarrow\; 0.
The iterates crowd together, so they converge to some limit x∗x^* in II; and since gg is continuous, taking the limit of xk+1=g(xk)x_{k+1} = g(x_k) gives x∗=g(x∗)x^* = g(x^*). That is existence.
Only one fixed point. If x∗x^* and y∗y^* were both fixed points, ∣x∗−y∗∣=∣g(x∗)−g(y∗)∣≤L ∣x∗−y∗∣|x^* - y^*| = |g(x^*) - g(y^*)| \le L\,|x^* - y^*|, which with L<1L < 1 forces ∣x∗−y∗∣=0|x^* - y^*| = 0.
How far from the answer. Letting m→∞m \to \infty in the second line, started one step earlier, gives a bound you can compute while iterating: ∣xk−x∗∣≤L1−L ∣xk−xk−1∣|x_k - x^*| \le \dfrac{L}{1 - L}\,|x_k - x_{k-1}|. A small step does not mean a small error when LL is close to 1 — a warning that returns in lesson 4-2.
Near the fixed point, the mean value theorem gives ek+1=g(xk)−g(x∗)=g′(ξk) eke_{k+1} = g(x_k) - g(x^*) = g'(\xi_k)\,e_k, so asymptotically
ek+1≈g′(x∗) eke_{k+1} \approx g'(x^*)\,e_k

Fixed-point iteration converges linearly with rate |g′(x*)|. The sign of g′ decides whether the iterates approach monotonically (g′ > 0) or alternate above and below (g′ < 0).

📐Why the calculator needs 58 presses
At x∗=0.7390851332x^* = 0.7390851332, ∣g′(x∗)∣=sin⁡x∗=0.6736|g'(x^*)| = \sin x^* = 0.6736. Each press multiplies the error by 0.67 and flips its sign — that is the swinging you saw. One decimal digit costs log⁡10/log⁡(1/0.6736)≈5.8\log 10/\log(1/0.6736) \approx 5.8 presses, so ten digits cost about 58. The measured error ratios in the lab are 0.6755, 0.6723, 0.6745, 0.6730, 0.6740. For x\sqrt{x}, the fixed point is 1 and g′(1)=1/2g'(1) = 1/2: one binary digit per press, like bisection. For x2x^2 at the fixed point 1, g′(1)=2g'(1) = 2: errors double, which is why 1.1 explodes and 0.9 collapses to the other fixed point, 0, where g′=0g' = 0.
Same equation, five iterations. x3+4x2−10=0x^3 + 4x^2 - 10 = 0 has a root at x∗=1.36523x^* = 1.36523. Here are five ways to write it as x=g(x)x = g(x), all started from 1.5. The last one is Newton's method. That is not a coincidence.
g(x)g′(x*)Behaviour
x−x3−4x2+10x - x^3 - 4x^2 + 10−15.5-15.5diverges in 7 steps
10/x−4x\sqrt{10/x - 4x}—square root of a negative number at step 3
1210−x3\tfrac12\sqrt{10 - x^3}−0.51-0.51converges in 34 steps
10/(4+x)\sqrt{10/(4 + x)}−0.13-0.13converges in 12 steps
x−x3+4x2−103x2+8xx - \dfrac{x^3 + 4x^2 - 10}{3x^2 + 8x}00converges in 5 steps

Five rearrangements of x³ + 4x² − 10 = 0, iterated from x₀ = 1.5 to a tolerance of 10⁻¹⁰

💡Newton as a fixed-point iteration
Newton's method is the fixed-point iteration for g(x)=x−f(x)/f′(x)g(x) = x - f(x)/f'(x). Differentiate: g′(x)=f(x) f′′(x)f′(x)2,g'(x) = \dfrac{f(x)\,f''(x)}{f'(x)^2}, which is zero at a simple root, where f(x∗)=0f(x^*) = 0. The linear term in the error vanishes, leaving the quadratic one — that is where quadratic convergence comes from. In general, if g′(x∗)=⋯=g(p−1)(x∗)=0g'(x^*) = \dots = g^{(p-1)}(x^*) = 0 and g(p)(x∗)≠0g^{(p)}(x^*) \ne 0, the iteration has order pp. And the failures of Newton from the last lesson are failures of contraction: far from the root, ∣g′∣>1|g'| > 1 and the map pushes points apart.
Derivation
Order of convergence from the derivatives of gg, and what a double root does
The derivative of Newton's map. With g(x)=x−f/f′g(x) = x - f/f', the quotient rule gives
g′(x)=1−f′ f′−f f′′f′2=f f′′f′2,g'(x) = 1 - \frac{f'\,f' - f\,f''}{f'^2} = \frac{f\,f''}{f'^2},
which is zero wherever f=0f = 0 and f′≠0f' \ne 0.
Why a zero derivative raises the order. Expand gg about the fixed point and use g(x∗)=x∗g(x^*) = x^*:
ek+1=g(xk)−g(x∗)=g′(x∗) ek+12g′′(x∗) ek2+16g′′′(x∗) ek3+…e_{k+1} = g(x_k) - g(x^*) = g'(x^*)\,e_k + \tfrac12 g''(x^*)\,e_k^2 + \tfrac16 g'''(x^*)\,e_k^3 + \dots
If g′(x∗)≠0g'(x^*) \ne 0 the first term dominates: linear convergence with rate ∣g′(x∗)∣|g'(x^*)|. If g′(x∗)=0g'(x^*) = 0 the second term leads: quadratic. In general, the first non-vanishing derivative g(p)(x∗)g^{(p)}(x^*) sets the order pp.
A root of multiplicity mm. Write f(x)=(x−x∗)mh(x)f(x) = (x - x^*)^m h(x) with h(x∗)≠0h(x^*) \ne 0. Then f′(x)=(x−x∗)m−1(m h+(x−x∗)h′)f'(x) = (x - x^*)^{m-1}\big(m\,h + (x - x^*)h'\big), and
g(x)=x−(x−x∗) hm h+(x−x∗)h′⟹g′(x∗)=1−1m.g(x) = x - \frac{(x - x^*)\,h}{m\,h + (x - x^*)h'} \qquad\Longrightarrow\qquad g'(x^*) = 1 - \frac{1}{m}.
For a double root g′(x∗)=12g'(x^*) = \tfrac12: Newton becomes a linear iteration that halves the error each step — no faster than bisection. If mm is known, the modified step xk+1=xk−mf/f′x_{k+1} = x_k - m f/f' makes g′(x∗)=0g'(x^*) = 0 and restores quadratic convergence.
Mathematics is the art of giving the same name to different things.
– Henri Poincaré, Science and Method (1908)
Accelerating a slow sequence. If the error really behaves like ek+1≈λeke_{k+1} \approx \lambda e_k, then three consecutive iterates xk,xk+1,xk+2x_k, x_{k+1}, x_{k+2} are enough to solve for x∗x^* and λ\lambda. The result is Aitken's Δ² process (Alexander Aitken, 1926):
x^k=xk−(xk+1−xk)2xk+2−2xk+1+xk\hat x_k = x_k - \frac{(x_{k+1} - x_k)^2}{x_{k+2} - 2x_{k+1} + x_k}

Aitken extrapolation. Exact for a sequence whose error is exactly geometric; for a linearly convergent iteration it converges much faster than the sequence itself.

Derivation
Solving for the limit from three iterates
Assume the error is exactly geometric for three consecutive iterates: xk+1−x∗=λ(xk−x∗)x_{k+1} - x^* = \lambda(x_k - x^*) and xk+2−x∗=λ(xk+1−x∗)x_{k+2} - x^* = \lambda(x_{k+1} - x^*). Dividing one by the other removes the unknown λ\lambda:
xk+2−x∗xk+1−x∗=xk+1−x∗xk−x∗⟹(xk+2−x∗)(xk−x∗)=(xk+1−x∗)2.\frac{x_{k+2} - x^*}{x_{k+1} - x^*} = \frac{x_{k+1} - x^*}{x_k - x^*} \quad\Longrightarrow\quad (x_{k+2} - x^*)(x_k - x^*) = (x_{k+1} - x^*)^2.
Multiply out. The x∗2x^{*2} terms cancel, leaving an equation linear in x∗x^*:
x∗ (xk−2xk+1+xk+2)=xkxk+2−xk+12⟹x∗=xkxk+2−xk+12xk+2−2xk+1+xk.x^*\,(x_k - 2x_{k+1} + x_{k+2}) = x_kx_{k+2} - x_{k+1}^2 \quad\Longrightarrow\quad x^* = \frac{x_kx_{k+2} - x_{k+1}^2}{x_{k+2} - 2x_{k+1} + x_k}.
Adding and subtracting xkx_k times the denominator in the numerator turns this into the form above, xk−(xk+1−xk)2/(xk+2−2xk+1+xk)x_k - (x_{k+1} - x_k)^2/(x_{k+2} - 2x_{k+1} + x_k), which is better in floating point: it computes a small correction to xkx_k instead of a ratio of two nearly equal large products (lesson 1-2). The denominator is a second difference, Δ2xk\Delta^2 x_k — hence the name.
Applied to the twelve cos iterates from 1.0, the error drops from 3.5×10−33.5\times10^{-3} to 8.8×10−68.8\times10^{-6} with no extra function evaluations.
Steffensen's method (Johan Steffensen, 1933) goes one step further: take two fixed-point steps, apply Aitken, and restart the iteration from the extrapolated value. The result converges quadratically — like Newton — without any derivative. On cos⁡x=x\cos x = x the errors go 1.1×10−21.1\times10^{-2}, 1.8×10−51.8\times10^{-5}, 4.9×10−114.9\times10^{-11}, 0.
🕰Polynomials: roots as eigenvalues
For polynomials there is a cleaner route than iterating on each root. The roots of xn+cn−1xn−1+⋯+c0x^n + c_{n-1}x^{n-1} + \dots + c_0 are exactly the eigenvalues of its n×nn \times n companion matrix. numpy.roots builds that matrix and calls an eigenvalue solver — it finds all roots at once, complex ones included. You met its sensitivity in Wilkinson's polynomial; Unit 6 shows how the eigenvalue solver itself works.
🛰Kepler's equation
Every orbit propagator solves E−esin⁡E=ME - e\sin E = M for the eccentric anomaly EE, millions of times. The obvious fixed-point form E←M+esin⁡EE \leftarrow M + e\sin E is a contraction with rate ecos⁡Ee\cos E. For a near-circular orbit that is tiny and a few steps suffice. For e=0.9e = 0.9 near periapsis (EE small, cos⁡E≈1\cos E \approx 1) the rate is 0.9 and the iteration needs 135 steps for 10−1210^{-12}; Newton needs 8. That is why flight software uses Newton — with a carefully chosen starting guess.
💡Where you will use this
Geotechnical Engineering: Bishop's method for slope stability has the factor of safety on both sides of its equation and is solved by exactly this iteration. Orbital Mechanics: Kepler's equation. Structures and CFD: every 'iterate until converged' loop — Picard iteration for nonlinear materials, pressure–velocity coupling — is a fixed-point iteration, and its convergence rate is ∣g′∣|g'|, measured from the residual history.
📚Further reading
  • Endre Süli and David Mayers, An Introduction to Numerical Analysis (Cambridge University Press, 2003), chapter 1, 'Solution of equations by iteration' — the contraction mapping theorem and convergence orders, with full proofs.
  • Michael Heath, Scientific Computing: An Introductory Survey (2nd ed.; SIAM Classics 2018), chapter 5 — fixed-point iteration, convergence rates and polynomial roots.
  • Alan Edelman and H. Murakami, 'Polynomial roots from companion matrix eigenvalues', Mathematics of Computation 64, 763–776 (1995) — why numpy.roots is reliable.
↓But there is a problem
Every method in this unit has chased one unknown at a time. On 9 November 1965, at 5:16 in the evening, a protective relay at a hydroelectric station in Ontario tripped a single transmission line. Power surged onto the neighbouring lines, they tripped too, and within minutes about 30 million people in Ontario and the north-eastern United States were in the dark, some for up to 13 hours. A power grid is thousands of buses, each with a voltage magnitude and angle, every one coupled to its neighbours through equations with sines and cosines in them. Knowing how power will flow through it means solving all of those equations at once.