Convergence of fixed point iteration¶

We revisit Fixed point iteration and investigate the observed convergence more closely. Recall that above we calculated $g'(r)\approx-0.42$ at the convergent fixed point.

using FundamentalsNumericalComputation

p = Polynomial([3.5,-4,1])
r = roots(p)
@show rmin,rmax = sort(r);

(rmin, rmax) = sort(r) = [1.2928932188134525, 2.7071067811865475]

Here is the fixed point iteration. This time we keep track of the whole sequence of approximations.

g = x -> x - p(x);
x = 2.1; 
for k = 1:12
    x = [x;g(x[k])]
end
x

      
    

13-element Array{Float64,1}:
 2.1
 2.59
 2.7419000000000002
 2.69148439
 2.713333728386328
 2.7044887203327885
 2.7081843632566587
 2.7066592708954196
 2.7072919457529734
 2.7070300492259465
 2.707138558717502
 2.707093617492436
 2.7071122335938966

It’s easiest to construct and plot the sequence of errors.

err = @. abs(x - rmax)
plot(0:12,err,m=:o,
    leg=:none,xaxis=("iteration number"),yaxis=("error",:log10),title="Convergence of fixed point iteration")

It’s quite clear that the convergence quickly settles into a linear rate. We could estimate this rate by doing a least-squares fit to a straight line. Keep in mind that the values for small $k$ should be left out of the computation, as they don’t represent the linear trend.

y = log.(err[5:12])
p = fit(5.0:12.0,y,1)

-0.6680573888302028 - 0.8807181589734455∙x

We can exponentiate the slope to get the convergence constant $\sigma$ .

sigma = exp(p.coeffs[2])

0.41448513854854724

The numerical values of the error should decrease by a factor of $\sigma$ at each iteration. We can check this easily with an elementwise division.

@. err[9:12] / err[8:11]

4-element Array{Float64,1}:
 0.4137660520817109
 0.4143987269383
 0.4141368304124451
 0.4142453399049934

The methods for finding $\sigma$ agree well.

Fundamentals of Numerical Computation

Convergence of fixed point iteration¶