Fixed Point Iteration

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:

f[x_]:=Cos[x]
Plot[{x,f[x]},{x,0,Pi/2}]

f[.5]
Cos[0.8]

f[%]

f[x_]:=-Sin[x]+.8+x
Plot[{x,f[x]},{x,0,Pi/2}]

f[.5]
ArcSin[0.8]

f[%]

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]

x0=0.7;
f[x_]:=-Sin[x]+x0+x
results=fixedPoint[f,.1,maxits->100,verbose->True,tolerance->10^-8];
points=Table[{i,results[[2]][[i]]},{i,1,Length[results[[2]]]}];
ListPlot[points,PlotStyle->{Thick,PointSize[Large]}]
ArcSin[x0]

This example illustrates chaos:

al=3.2;
x0=0.1;
f[x_]:=al* x (1-x)
results=fixedPoint[f,x0,maxits->100,verbose->True,tolerance->10^-8];
points=Table[{i,results[[2]][[i]]},{i,1,Length[results[[2]]]}];
ListPlot[points,PlotStyle->{Thick,PointSize[Large]}]

Danny suggested a variation on my fixed point function -- scaling it by a parameter. That can be beneficial (but not at the tricky endpoint...).

f[x_,a_,b_]:=a(-Sin[x]+b)+x
fp[x_,a_]:=-a Cos[x]+1
Manipulate[
Module[{},
Plot[{x,f[x,a,b]},{x,0,Pi/2},PlotLabel->fp[Pi/2,a]]
],
{a,0.01,2},
{b,0.01,1}
]