Fixed Point Iteration

Andy Long, 2/16/15
Thanks to: 
1. http://www.math.missouri.edu/~stephen/num-methods/2.2.iterate.nb
2. http://mathworld.wolfram.com/notebooks/DynamicalSystems/LogisticMap.nb

Here we are going to illustrate fixed point iteration using mathematica.  Let us set up a function for which we desire a root:

  f(x) = cos(x).  

We will work in the interval a<=x<=b, where a = 0, and b = 1.  Our starting value is (Subscript[p, 0]) = Pi/8.

Clear[f,g,a,b,p0]

We’re going to find a root of cosine, but in an interesting way!

f[x_]=Cos[x];
g[x_]=f[x]+x;
p0=Pi/8;
a=0;
b=Pi;
precision=50;

Now let us do the iteration:

Clear[p,k]
p[0]=N[p0,precision];
Do[p[k]=N[g[p[k-1]],precision],{k,1,100}]

Let’s write out the values:

TableForm[Table[{k,NumberForm[p[k],precision]},{k,0,11}]]

Now let’s do a bunch of plots, and then make a movie.  What do you think the Mathematica code is doing here?

FixedPointSegments[g_,x0_,its_]:=Module[
{l=NestList[g[#]&,x0,its]},
Partition[Flatten[Transpose[{l,l}]],2,1]
];
MakeMovie[g_,p0_,a_,b_,steps_]:=
Module[{p2=Partition[FixedPointSegments[g,p0,steps],2,1]},
movie=Table[Plot[{g[x],x},{x,a,b},PlotStyle->{{},Dashing[{.02}]},
Prolog->MapIndexed[{Hue[#2[[1]]/Length[p2]],Line[#]}&,Take[p2,i]],
AspectRatio->Automatic,
AxesLabel->TraditionalForm/@{(Subscript[x, n]),(Subscript[x, n+1])}
],{i,Length[p2]}];
ListAnimate[movie]
]

MakeMovie[g,N[p0],a,b,5]

We’re going to find a root of Exp[-x]-.5,still in an interesting way.

f[x_]=Exp[-x]-.5;
g[x_]=f[x]+x;
p0=1.7;
a=0;
b=2;
precision=50;
Clear[p,k]
p[0]=N[p0,precision];
Do[p[k]=N[g[p[k-1]],precision],{k,1,100}]
(* Let’s write out the values: *)
TableForm[Table[{k,NumberForm[p[k],precision]},{k,0,11}]]

MakeMovie[g,N[p0],a,b,10]

Let’s try to find a root of ln(x)....

f[x_]=Log[x]/5;
g[x_]=f[x]+x;
p0=1.1;
a=0.1;
b=2;
precision=50;
Clear[p,k]
p[0]=N[p0,precision];
Do[p[k]=N[g[p[k-1]],precision],{k,1,100}]
(* Let’s write out the values: *)
TableForm[Table[{k,NumberForm[p[k],precision]},{k,0,11}]]

MakeMovie[g,N[p0],a,b,15]

Let’s try another function: this one won’t be constructed in exactly the same way, however: we could say we’re trying to find a solution to  
x=(1/2)+((Sin[x])/2)

Now let’s try other functions.  When does the sequence p[k] converge?

Clear[p,g,p0,a,b]
g[x_]=(1/2)+((Sin[x])/2);
p0=0;
a=0;
b=1;
Clear[p,k]
p[0]=N[p0,precision];
Do[p[k]=N[g[p[k-1]],precision],{k,1,100}]
(* Let’s write out the values: *)
TableForm[Table[{k,NumberForm[p[k],precision]},{k,0,11}]]

MakeMovie[g,p0,a,b,15]

A nice quadratic:

Clear[p,g,p0,a,b,k]
f[x_]=1-(x^2);
g[x_]=f[x]+x;
p0=0.6;
a=0;
b=2;
Clear[p,k]
p[0]=N[p0,precision];
Do[p[k]=N[g[p[k-1]],precision],{k,1,100}]
(* Let’s write out the values: *)
TableForm[Table[{k,NumberForm[p[k],precision]},{k,0,11}]]

MakeMovie[g,p0,a,b,25]

A tricky quadratic:

Clear[p,g,p0,a,b,k]
f[x_]=1-(x^2)-x;
g[x_]=f[x]+x;
p0=.5;
a=-2;
b=2;
Clear[p,k]
p[0]=N[p0,precision];
Do[p[k]=N[g[p[k-1]],precision],{k,1,100}]
(* Let’s write out the values: *)
TableForm[Table[{k,NumberForm[p[k],precision]},{k,0,11}]]

MakeMovie[g,p0,a,b,45]