Fixed Point Iteration Acceleration

Andy Long, Spring, 2024

Here’s our author’s picture of two scenarios:  on the left, convergence; on the right, failure to converge.

Let’s try to use fixed point iteration to solve for Arcsin:

fixedPoint code:

(* These are options that one can give at game time: 
-- epsilon: the value we consider "zero" (0 by default)
  -- tolerance: the width of the interval we consider "close enough" to the root
		-- verbose: print the details of the convergence
-- debug: print out a little more for the coder....
*)
Options[fixedPoint]={maxits->100,epsilon->0,tolerance->10^(-6),verbose->True,debug->False} ;

(* Here's the function itself: *)
fixedPoint[f_,a0_,OptionsPattern[]] :=

Module[{a=a0,c,fa=f[a0],fc,slope,width,output,i=1,its=OptionValue[maxits],eps=OptionValue[epsilon],tol=OptionValue[tolerance],verbose=OptionValue[verbose],debug=OptionValue[debug]},
width=10000;
(* Update and continue: let's start looping: *)
While[i<=its && width>tol,
(* c is the function value at a: *)
c=fa;
fc=f[c];
width=Abs[c-a];
slope=(fc-fa)/width;

(* Otherwise, store output and do it again: *)
If[i==1,
(* Let's collect the output, since we've gone to the trouble of computing it: *)
output={{a,c,fc,slope}},
output=Append[output,{a,c,fc,slope}]
]; 

(* Did we stomp on a root by good fortune? *)
If[width<=tol, 
If[verbose,
Print["Root at c."];printoutFixedPoint[output,a,c,fc,width,slope];Return[Transpose[output]],
Return[fc] (* if we're converging, this is better than c *)
]
];

(* update: *)
a=c;
fa=fc;
i=i+1;
]; 

(* We've popped out, so print some nice formatted output: *)
If[verbose,
printoutFixedPoint[output,a,c,fc,width,slope];
(* Return the output, so that it can be plotted, etc. *)
Return[Transpose[output]],
Return[fc]
];
]

This is the printout command for verbose results:

printoutFixedPoint[output_,a_,c_,fc_,width_,slope_]:=
Module[{a0=output[[1]][[1]]},
Print[a0];
Print[NumberForm[TableForm[output,
TableHeadings->{None,{"(Subscript[a, k])","(Subscript[c, k])","f[(Subscript[c, k])]","slope"}}],16]]; 
Print["  c  = ",NumberForm[c,16] ]; 
Print[" Δc  = ±",width]; 
Print["f[c] = ",NumberForm[fc,16] ]; 
Print["f'[c] = ",NumberForm[slope,16] ]; 
]

We  illustrate that the method works (at least on these problems):

f[x_]:=Cos[x]
results=fixedPoint[f,0.5,maxits->100,verbose->False];
results

f[x_]:=-Sin[x]+.7+x
results=fixedPoint[f,.1,maxits->100,verbose->False,tolerance->10^-8]
results
ArcSin[.7]

Let’s accelerate the convergence:

f[x_]:=Cos[x]
results=fixedPoint[f,0.5,maxits->100,verbose->True];
iterates=results[[1]]

aitkensStep[a_,b_,c_]:=a-(b-a)^2/(c-2b+a)
aitkens[iterates_]:=
Table[aitkensStep[iterates[[i-2]],iterates[[i-1]],iterates[[i]]],{i,3,Length[iterates]}]
improved=aitkens[iterates]
ListPlot[Transpose[{improved,Rest[Rest[iterates]]}]]
ListPlot[{improved,Rest[Rest[iterates]]}]

andysStep[a_,b_,c_]:=c+(c-a)(c-b)^2/(a^2-2b(a-c)-c^2)
andysStep[a_,b_,c_]:=c-(c-b)^2/(a-2b+c)
andys[iterates_]:=
Table[andysStep[iterates[[i-2]],iterates[[i-1]],iterates[[i]]],{i,3,Length[iterates]}]
better=andys[iterates]
ListPlot[Transpose[{better,Rest[Rest[iterates]]}]]
ListPlot[{better,Rest[Rest[iterates]]}]

ListPlot[{better-improved}]

x0=0.99
f[x_]:=-Sin[x]+x0+x
results=fixedPoint[f,.1,maxits->100,verbose->True,tolerance->10^-8];
iterates=results[[1]];
ArcSin[x0]

better=andys[iterates];
ListPlot[Transpose[{better,Rest[Rest[iterates]]}]]
ListPlot[{better,Rest[Rest[iterates]]}]

steffensen[g_,x0_,tol_,n_]:=Module[{i=0,a,x1,x2,q0,y1,y2,y0,denom},
q0=x0;
x1=N[g[q0]];
x2=g[x1];
y1=x1;
y2=x2;
denom=(x2-2*x1+q0);
While[Abs[denom]>tol && i<n,
q0=q0-(x1-q0)^2/denom;
a=y2-(y1-y2)^2/(y2-2*y1+y0);
y0=y2;
y1=g[y2];
y2=g[y1];
Print[{y2,a,q0}];
x1=g[q0];
x2=g[x1];
q0=x2;
denom=(x2-2*x1+q0);
i=i+1
];
q0
]
steffensen[Cos,0.999,10^-8,30]