Module for Fixed Point Iteration

Numerical Methods using Mathematica  
Complementary software supplements for   
Numerical Methods textbooks  
John H. Mathews and Kurtis Fink, 2004

A fundamental principle in computer science is iteration.  As the name suggests, a process is repeated until an answer is achieved. Iterative techniques are used to find roots of equations, solutions of linear and nonlinear systems of equations, and solutions of differential equations.

A rule or function g(x) for computing successive terms is needed, together with a starting value (Subscript[p, 0]).  Then a sequence of values {(Subscript[p, k])} is obtained using the iterative rule (Subscript[p, k+1])=g((Subscript[p, k])).   The sequence has the pattern
	   (Subscript[p, 0])			(starting value)
	   (Subscript[p, 1])=g((Subscript[p, 0]))
	   (Subscript[p, 2])=g((Subscript[p, 1]))
	        ⋮  
	   (Subscript[p, k])=g((Subscript[p, k-1]))
	(Subscript[p, k+1])=g((Subscript[p, k]))
	        ⋮

What can we learn from an unending sequence of numbers?  If the numbers tend to a limit, we suspect that it is the answer.

Finding Fixed Points

Definition (FixedPoint).  A fixed point of a function g(x) is a number P such that P=g(P).  

Caution.  A fixed point is not a root of the equation  0=g(x),  it is a solution of the equation  x=g(x).

Geometrically, the fixed points of a function  g(x)  are the point(s) of intersection of the curve  y=g(x)  and the line  y=x.

Definition (Fixed Point Iteration).  The iteration (Subscript[p, n])=g((Subscript[p, n-1])) for n=0,1,... is called fixed point iteration.

Theorem (For a converging sequence).  Assume that g(x) is a continuous function and that (Subscript[{(Subscript[p, n])}, n=0]^∞) is a sequence generated by fixed point iteration.

If  (Underscript[ lim , n->∞])(Subscript[p, n])=P,  then P is a fixed point of g(x).

Proof  Fixed Point Iteration  Fixed Point Iteration

The following two theorems establish conditions for the existence of a fixed point and the convergence of the fixed-point iteration process to a fixed point.

Theorem (First Fixed Point Theorem).  Assume that g∈C[a,b], i. e.  g(x)  is continuous on  [a,b].  
Then we have the following conclusions. 
(i).	If the range of the mapping y=g(x) satisfies y∈[a,b] for all x∈[a,b], then  g has a fixed point in [a,b]. 
(ii).	Furthermore, suppose that g'(x) is defined over (a,b) and that a positive constant K<1 exists with 
	|g'(x)|≤ K  for all  x∈(a,b),  then g has a unique fixed point P in [a,b].

Proof  Fixed Point Iteration  Fixed Point Iteration

Theorem (Second Fixed Point Theorem).  Assume that the following hypothesis hold true.  
(a)	P is a fixed point of a function g,
(b)	g,g'∈C[a,b], 
(c)	K is a positive constant,
(d)	(Subscript[p, 0])∈(a,b), and 
(e)	g(x)∈[a,b]  for all  x∈[a,b].  
Then we have the following conclusions. 
(i).	If |g'(x)|≤ K<1  for all   x∈[a,b],  then the iteration  (Subscript[p, n])=g((Subscript[p, n-1]))  will converge to the 
	unique fixed point P∈(a,b).  In this case, P is said to be an attractive fixed point.  
(ii).	If |g'(x)|> 1  for all   x∈[a,b],  then the iteration  (Subscript[p, n])=g((Subscript[p, n-1]))  will not converge to P.  
	In this case, P is said to be a repelling fixed point and the iteration exhibits local divergence.

Proof  Fixed Point Iteration  Fixed Point Iteration

Remark 1.  It is assumed that (Subscript[p, 0])≠P in statement (ii). 

Remark 2.  Because g  is continuous on an interval containing P, it is permissible to use the simpler criterion |g'(P)|≤ K<1  and  |g'(P)|> 1 in (i) and (ii), respectively.

Corollary.  Assume that g  satisfies hypothesis (a)-(e)  of the previous theorem.  Bounds for the error involved when using  (Subscript[p, n])  to approximate  P  are given by         

	|P - (Subscript[p, n])| ≤ (K^n) |P - (Subscript[p, 0])|          	for  n≥1,  
and
	|P - (Subscript[p, n])| ≤ ((K^n)/(1-K)) |(Subscript[p, 1]) - (Subscript[p, 0])|  	for  n≥1.

Graphical Interpretation of Fixed-point Iteration

Since we seek a fixed point P to g(x),  it is necessary that the graph of the curve  y=g(x)  and the line  y=x  intersect at the point (P,P).  
The following animations illustrate two types iteration: monotone and oscillating.

Algorithm (Fixed Point Iteration).  To find a solution to the equation  x=g(x) by starting with  (Subscript[p, 0]) and iterating  (Subscript[p, n])=g((Subscript[p, n-1])).

Computer Programs  Fixed Point Iteration  Fixed Point Iteration

Mathematica Subroutine (Fixed Point Iteration).

FixedPointIteration[x0_,max_]:=
Module[{},
(Subscript[p, 0]) = N[x0]; 
k = 0; 
Print[(Subscript["  p", 0])," = ",PaddedForm[(Subscript[p, 0]),{15,15}]]; 
While[k<max,
Module[{},
(Subscript[p, 1]) = g[(Subscript[p, 0])]; 
k = k+1; 
Print[(Subscript["  p", k])," = ",PaddedForm[(Subscript[p, 1]),{15,15}] ]; 
(Subscript[p, 0]) = (Subscript[p, 1]); ]; ]; 
p = (Subscript[p, 0]); 
Print[" "]; 
Print["The function is g[x] = ",g[x]]; 
Print["  p  = ",PaddedForm[p,{15,15}]]; 
Print["g[p] = ",PaddedForm[g[p],{15,15}]];  ];

Example 1.  Use fixed point iteration to find the fixed point(s) for the function   g(x)=1+x-((x^2)/3).

Solution 1.

Plot the function and determine graphically that there are two solutions to the equation  x=g(x).

(g[x_]=1+x-((x^2)/3);) (Needs["Graphics`Colors`"];) (gr1=Plot[{x,g[x]},{x,-2,4},PlotStyle->{Green,Magenta}]) (Print["The function is g[x] = ",g[x]];)

Will iteration prove successful as a computational tool to numerically find these two fixed points.

(g[x_]=1+x-((x^2)/3);) (Needs["Graphics`Colors`"];) (gr1=Plot[{x,g[x]},{x,0,4.82`},PlotRange->{{0,4},{0,3}},AspectRatio->(3/4),PlotStyle->{Green,Magenta}]) (Print["The function is g[x] = ",g[x]];)

Use fixed point iteration to find a numerical approximation.
First, do the iteration one step at a time.  Type each of the following commands in a separate cell and execute them one at a time.

(Subscript[p, 0])=3.0

(Subscript[p, 1])=g[(Subscript[p, 0])]

(Subscript[p, 2])=g[(Subscript[p, 1])]

(Subscript[p, 3])=g[(Subscript[p, 2])]

(Subscript[p, 4])=g[(Subscript[p, 3])]

(Subscript[p, 5])=g[(Subscript[p, 4])]

Now use the FixedPointIteration subroutine to perform the computations.

FixedPointIteration[3.0,7];

Are 7 iterations enough to locate the fixed point? 

How many decimal places agreement do you see to the equation  (Subscript[p, 7])=g[(Subscript[p, 6])] ?

If it is not sufficient for your purposes you must request more iterations!

FixedPointIteration[3.0,20];

Comparison 1. Compare with Mathematica's built in "FindRoot" subroutine for numerically finding solutions.

eqn=g[x]==x;  
solset = FindRoot[eqn,{x,0.5}];  
p=(Subscript[solset, [[1,2]]]);  
Print["The function is g[x] = ",g[x]];  
Print[""];  
Print[eqn];  
Print[solset];  
Print[""];  
Print["   p = ",NumberForm[p,16] ];  
Print[""];  
Print["g[p] = ",NumberForm[g[p],16] ];

Did Mathematica get all of the digits?  Why?
There are options to every built in Mathematica subroutine. For FindRoot it is:
Options[FindRoot]
{AccuracyGoal->Automatic,Compiled->True,DampingFactor->1, Jacobian->Automatic,MaxIterations->15,WorkingPrecision->16}
It is the purpose of numerical analysis to study the last two items "MaxIterations,"  "WorkingPrecision."

Comparison 2. Compare with Mathematica's built in "Solve" subroutine for symbolically finding solutions.

eqn=g[x]==x;  
solset = Solve[g[x]==x,x];  
p=(Subscript[solset, [[2,1,2]]]);  
Print["The function is g[x] = ",g[x]];  
Print[""];  
Print[eqn];  
Print[solset];  
Print[""];  
Print["   p = ",p];  
Print[""];  
Print["g[p] = ",g[p]];

Mathematica's answers are "mathematically perfect."

Now investigate fixed point iteration near the other fixed point,  p=-(Sqrt[3]).

(gr2=Plot[{x,g[x]},{x,-7,0},PlotRange->{{-7,0},{-7,0}},AspectRatio->1,PlotStyle->{Green,Magenta}]) (Print["The function is g[x] = ",g[x]];)

Now use the FixedPointIteration subroutine to perform the computations.

FixedPointIteration[-2.0,5];

Remark. You need to look carefully at the above output and determine what it means !  Did it converge ?

The distinguishing property for determining convergence is the size of  |g'[p]|.

If  (Subscript[p, 0])  is near the fixed point  p  and    |g'[p]|< 1  then the iteration will converge to  p.

If  (Subscript[p, 0])  is near the fixed point  p  and    |g'[p]|> 1  then the iteration will not converge to  p.

Print["The  function  is  g[x] = ",g[x]];  
Print["The derivative is g'[x] = ",g'[x]];  
Print[""];  
p=(Sqrt[3]);  
Print["  g'[",p,"] = ",g'[p] ];  
Print["|g'[",p,"]| = ",Abs[g'[p]] ];  
Print["|g'[",p,"]| = ",N[Abs[g'[p]]] ];  
Print[""];  
p=-(Sqrt[3]);  
Print["  g'[",p,"] = ",g'[p] ];  
Print["|g'[",p,"]| = ",Abs[g'[p]] ];  
Print["|g'[",p,"]| = ",N[Abs[g'[p]]] ];

The distinguishing property for determining convergence is the size of  |g'[p]|.

If  (Subscript[p, 0])  is near the fixed point  p=(Sqrt[3])  and  |g'[x]|< 1  in the neighborhood, then the iteration will converge to  p=(Sqrt[3]).

If  (Subscript[p, 0])  is near the fixed point  p=-(Sqrt[3])  and  |g'[x]|> 1  in the neighborhood, then the iteration will not converge to  p=-(Sqrt[3]).

Example 2.  Use Mathematica's built in subroutine "FixedPointList" and experiment finding the fixed point(s) for the function   g(x)=1+x-((x^2)/4).

Solution 2.

Investigate iteration near the "attractive fixed point."

g[x_]=1+x-((x^2)/3);  
Print["The function is g[x] = ",g[x]];

P=FixedPointList[g,3.0,5]

How could we generate this list ?  Execute the following cells one cell at a time. 
It will start with a list of one number and "Append" numbers to the list.

P={3.0};  
For[ i=2,i≤8,i++,
P=Append[P,g[(Subscript[P, [[-1]]])] ] ];  
P

Thus we have constructed the essential ideas for Mathematica's  "FixedPointList" subroutine.

Aside. It has become "pop math" to see all those crazy figures about "chaos", etc. They all depend on drawing line segments between points. This graphical excursion is not the usual academic treatment of numerical analysis.  It is just for fun!

pts={{(Subscript[P, [[1]]]),0}};  
For[ i=1,i≤Length[P]-1,i++,
pts=Append[pts,{(Subscript[P, [[i]]]),(Subscript[P, [[i+1]]])} ];  
pts=Append[pts,{(Subscript[P, [[i+1]]]),(Subscript[P, [[i+1]]])} ];  
 ];  
pts

Consider the following graph consisting of line segments joining the points given above.

lin1=ListPlot[pts,Joined->True,PlotRange->{{0,4.82`},{0,2.1`}},AspectRatio->((2.1`)/(4.82`)),PlotStyle->Blue]

Show this graph and the one in Example 1 on the same plot.

Show[gr1,lin1,PlotRange->{{0,4},{0,3}},AspectRatio->(3/4)] (Print["The function is g[x] = ",g[x]];) (Print["The sequence of fixed point iterates is"];) (Print["P = ",P];)

Example 3.  Use Mathematica's built in subroutine "FixedPointList" and continue investigation of the "repulsive fixed point" for the function   g(x)=1+x-((x^2)/3).

Solution 3.

g[x_]=1+x-((x^2)/3);  
Print["The function is g[x] = ",g[x]];

P=FixedPointList[g,-1.8,10]

Let's generate a list of points and plot that nifty graph.

pts={{(Subscript[P, [[1]]]),0}}; 
For[i=1,i≤Length[P]-1,i++,
pts=Append[pts,{(Subscript[P, [[i]]]),(Subscript[P, [[i+1]]])} ]; 
pts=Append[pts,{(Subscript[P, [[i+1]]]),(Subscript[P, [[i+1]]])} ];  
 ];  
pts

Consider the following graph consisting of line segments joining the points given above.

lin2=ListPlot[pts,Joined->True,PlotRange->{{-7,0},{-7,0}},AspectRatio->1,PlotStyle->Blue]

Show this graph and the one in Example 1 on the same plot.

Show[gr2,lin2,PlotRange->{{-7,0},{-7,0}},AspectRatio->1] (Print["The function is g[x] = ",g[x]];) (Print["The sequence of fixed point iterates is"];) (Print["P = ",P];)

Various Scenarios and Animations for Fixed Point Iteration.

Example 4.  Convergence: Monotone Increasing  Find the solution to  x = g[x] =(Sqrt[2x]).  Use the starting approximation  (Subscript[p, 0])=0.1

Solution 4.

(Needs["Graphics`Colors`"];) (g[x_]=(Sqrt[2 x]);) (p0=0.1`;) (FixedPointIteration[p0,20];) Plot[{g[x],x},{x,0.`,4.05`},PlotStyle->{Magenta,Green}]

Print[" g[2] = ",g[2]]; 
Print["g'[2] = ",g'[2]];

Since |g'(x)|< 1  for all   x∈((1/2),∞),  the iteration  (Subscript[p, n])=g((Subscript[p, n-1]))  will converge to the fixed point  p=2.  In this case,  p=2  is said to be an attractive fixed point.

Needs["Algebra`InequalitySolve`"]; 
Print["g'[x] = ",g'[x]]; 
Print["Solve  -1 < g'[x] < 1"]; 
solset=InequalitySolve[-1<g'[x]<1,x];
Print["Get    ",solset];

Animation 1.  ( Fixed Point Iteration  Fixed Point Iteration ).  Internet hyperlink.

Example 5.  Convergence: Monotone Decreasing  Find the solution to  x = g[x] =(Sqrt[2x]).  Use the starting approximation  (Subscript[p, 0])=4.0

Solution 5.

(Needs["Graphics`Colors`"];) (g[x_]=(Sqrt[2 x]);) (p0=4.`;) (FixedPointIteration[p0,20];) Plot[{g[x],x},{x,0.`,4.05`},PlotStyle->{Magenta,Green}]

Print[" g[2] = ",g[2]]; 
Print["g'[2] = ",g'[2]];

Since |g'(x)|< 1  for all   x∈((1/2),∞),  the iteration  (Subscript[p, n])=g((Subscript[p, n-1]))  will converge to the fixed point  p=2.  In this case,  p=2  is said to be an attractive fixed point.

Needs["Algebra`InequalitySolve`"]; 
Print["g'[x] = ",g'[x]]; 
Print["Solve  -1 < g'[x] < 1"]; 
solset=InequalitySolve[-1<g'[x]<1,x];
Print["Get    ",solset];

Animation 2.  ( Fixed Point Iteration  Fixed Point Iteration ).  Internet hyperlink.

Example 6. Convergence: Spiral  Find the solution to  x = g[x] = (8/(x+2)).  Use the starting approximation  (Subscript[p, 0])=4.0

Solution 6.

(Needs["Graphics`Colors`"];) (g[x_]=(8/(x+2));) (p0=4.`;) (FixedPointIteration[p0,20];) Plot[{g[x],x},{x,0.`,4.05`},PlotStyle->{Magenta,Green}]

Print[" g[2] = ",g[2]]; 
Print["g'[2] = ",g'[2]];

Since |g'(x)|< 1  for all   x∈(2 (-1+(Sqrt[2])),∞)=(0.828427,∞),  the iteration  (Subscript[p, n])=g((Subscript[p, n-1]))  will converge to the fixed point  p=2.  In this case,  p=2  is said to be an attractive fixed point.

Needs["Algebra`InequalitySolve`"]; 
Print["g'[x] = ",g'[x]]; 
Print["Solve  -1 < g'[x] < 1"]; 
solset=InequalitySolve[-1<g'[x]<1,x];
Print["Get    ",solset];

Animation 3.  ( Fixed Point Iteration  Fixed Point Iteration ).  Internet hyperlink.

Example 7. Divergence: Monotone Increasing  Find the solution to  x = g[x] = ((x^2)/4)+(x/2).  Use the starting approximation  (Subscript[p, 0])=2.03

Solution 7.

(Needs["Graphics`Colors`"];) (g[x_]=((x^2)/4)+(x/2);) (p0=2.03`;) (FixedPointIteration[p0,15];) Plot[{g[x],x},{x,0.`,4.05`},PlotStyle->{Magenta,Green}]

Print[" g[2] = ",g[2]]; 
Print["g'[2] = ",g'[2]];

Since|g'(x)|> 1  for all   x∈(1,∞),  then the iteration  (Subscript[p, n])=g((Subscript[p, n-1]))  will not converge to  p=2 .  In this case,   p=2 is said to be a repelling fixed point and the iteration exhibits local divergence.

Needs["Algebra`InequalitySolve`"]; 
Print["g'[x] = ",g'[x]]; 
Print["Solve  1 < g'[x]"]; 
solset=InequalitySolve[1<g'[x],x];
Print["Get    ",solset];

Animation 4.  ( Fixed Point Iteration  Fixed Point Iteration ).  Internet hyperlink.

Example 8. Convergence: Monotone Decreasing  Find the solution to  x = g[x] = ((x^2)/4)+(x/2).  Use the starting approximation  (Subscript[p, 0])=1.95

Solution 8.

(Needs["Graphics`Colors`"];) (g[x_]=((x^2)/4)+(x/2);) (p0=1.95`;) (FixedPointIteration[p0,30];) Plot[{g[x],x},{x,0.`,4.05`},PlotStyle->{Magenta,Green}]

Print[" g[0] = ",g[0]]; 
Print["g'[0] = ",g'[0]];

Since |g'(x)|< 1  for all   x∈(-3,1),  the iteration  (Subscript[p, n])=g((Subscript[p, n-1]))  will converge to the fixed point  p=0.  In this case,  p=0  is said to be an attractive fixed point.

Needs["Algebra`InequalitySolve`"]; 
Print["g'[x] = ",g'[x]]; 
Print["Solve  -1 < g'[x] < 1"]; 
solset=InequalitySolve[-1<g'[x]<1,x];
Print["Get    ",solset];

Print[" g[2] = ",g[2]]; 
Print["g'[2] = ",g'[2]];

Since|g'(x)|> 1  for all   x∈(1,∞),  then the iteration  (Subscript[p, n])=g((Subscript[p, n-1]))  will not converge to  p=2 .  In this case,   p=2 is said to be a repelling fixed point and the iteration exhibits local divergence.

Needs["Algebra`InequalitySolve`"]; 
Print["g'[x] = ",g'[x]]; 
Print["Solve  g'[x] > 1"]; 
solset=InequalitySolve[1<g'[x],x];
Print["Get    ",solset];

Animation 5.  ( Fixed Point Iteration  Fixed Point Iteration ).  Internet hyperlink.

Example 9. Divergence: Spiral  Find the solution to  x = g[x] = -((x^2)/4)-(x/2)+4.  Use the starting approximation  (Subscript[p, 0])=1.97

Solution 9.

(Needs["Graphics`Colors`"];) (g[x_]=-((x^2)/4)-(x/2)+4;) (p0=1.97`;) (FixedPointIteration[p0,20];) Plot[{g[x],x},{x,0.`,4.05`},PlotStyle->{Magenta,Green}]

Print[" g[2] = ",g[2]]; 
Print["g'[2] = ",g'[2]];

Since|g'(x)|> 1  for all   x∈(1,∞),  then the iteration  (Subscript[p, n])=g((Subscript[p, n-1]))  will not converge to  p=2 .  In this case,   p=2 is said to be a repelling fixed point and the iteration exhibits local divergence.

Needs["Algebra`InequalitySolve`"]; 
Print["g'[x] = ",g'[x]]; 
Print["Solve  g'[x] < -1"]; 
solset=InequalitySolve[g'[x]<-1,x];
Print["Get    ",solset];

Animation 6.  ( Fixed Point Iteration  Fixed Point Iteration ).  Internet hyperlink.

Animations (Fixed Point Iteration  Fixed Point Iteration).  Internet hyperlinks to animations.

Research Experience for Undergraduates

Fixed Point Iteration  Fixed Point Iteration  Internet hyperlinks to web sites and a bibliography of articles.