Numerical Analysis 5

Spring 2005

Chapter 2: Solving x=g(x)

Iteration

You can often find a solution of  x = g (x) by starting with some initial (Subscript[p, 0]) and forming the sequence of iterates, (Subscript[p, n + 1])  =  g ((Subscript[p, n])).  Terminate the iterations when the relative error estimate,   2|(((Subscript[p, n + 1]) -  (Subscript[p, n]))/(small +  ∣ 
                (Subscript[p, n + 1]) ∣ ))|,   in consecutive iterates  (Subscript[p, 0])  and  (Subscript[p, n + 1])  differs by less than a preassigned value called tolerance.  The next Mathematica function is a modification of the function formulated by the author of your textbook.  Following the function form, definitions of function arguments are given 
FixedPointIteration[g0_,X0_,tolerance_,max_,digits_,{leftx_,rightx_},{lowery_,uppery_}]
g0_ is the function g(x),
X0_ is the starting value for the iterations,
tolerance_  is the preassigned tolerance,
max_  is the maximum number of iterations allowed,
digits_ is the number of digits displayed in the output.
leftx_ is the leftmost x-value displayed on the graph
rightx_  is the rightmost x-value displayed on the graph
lowery_ is the minimum y-value displayed on the graph
uppery_ is the maximum y-value displayed on the graph

\(\(FixedPointIteration[g0_, X0_, tolerance_, max_, 
          digits_, {leftx_, rightx_}, {lowery_, uppery_}] := 
	
        Module[{g, k, x0, x1, relerr, plist, doubleplist, gr1, gr2, gr3, 
            gr4}, 
 (*Initial\ Values\ for\ Variables*) 
	
          Set @@ {g[x_], g0}; 
		x0 = 
            N[X0, digits]; 
	\ \ plist = {x0}; 
		k = 
            0; 
		err = 1. ; 
		relerr = 1; 
		small = 
            N[1\/10\^20, 
              digits]; 
Print[\*""Iterates (Subscript[p, k]),  "", "    Relative Error |Δp|"]; 
		Print["p"\_0, " = ", N[x0, digits]]; 
 (*Iteration*) 
		While[
            k < max && 
              tolerance ≤ 
                relerr, 
			Module[{}, 
				x1 = 
                g[x0]; 
				k = k + 1; 
				err = 
                Abs[x1 - x0]; 
				Print["p"\_k, " = ", 
                N[x1, digits], ",   |Δp| ≤ ", 
                N[err, digits]]; 
			relerr = \(2  Abs[x1 - x0]\)\/\(Abs[x1] + small\); 
				x0 = x1; 
              plist = Append[plist, x1]]]; 
		p = 
            x0; 
		Print[" "]; 
 (*Collect\ approximations\ into\ a\ list\ to\ form\ points\ for\ graphing*) 
	
          doubleplist = 
            Flatten[Transpose[{plist, plist}]]; 
	
          pts = Transpose[{Drop[doubleplist, {Length[doubleplist]}], 
                Prepend[Drop[doubleplist, 2], 
                  0]}]; 
 (*Print\ Results*) 
		Print["The function is g[x] = ", g[x]]; 
          Print["Approximation to a fixed point of g[x] is:"]; 
          Print["  p  = ", N[p, digits]]; 
          Print["g[p] at the approximated fixed point= ", 
            N[g[p], digits]]; 
          Print["|Δp| ≤ ", 
            N[err, digits]]; 
 (*Graph\ Results*) 
	
          gr1\  = Plot[g[x], {x, leftx, rightx}, 
              PlotRange -> {{leftx, rightx}, {lowery, uppery}}, 
              PlotStyle -> RGBColor[1, 0, 0], 
              DisplayFunction -> Identity]; 
		gr2\  = 
            Plot[x, {x, leftx, rightx}, 
              PlotRange -> {{leftx, rightx}, {lowery, uppery}}, 
              PlotStyle -> RGBColor[0, 1, 0], 
              DisplayFunction -> Identity]; 
	
          gr3\  = ListPlot[pts, 
              PlotRange -> {{leftx, rightx}, {lowery, uppery}}, 
              PlotStyle -> {PointSize[0.01]}, 
              DisplayFunction -> Identity]; 
	
          gr4 = Show[Graphics[Line[pts]], 
              DisplayFunction -> Identity]; 
		Show[{gr1, gr2, gr3, 
              gr4}, DisplayFunction -> $DisplayFunction];];\)\)

FixedPointIteration[Exp[-x], 0.5, 0.01, 30, 
        12, {0, 1}, {0, 1}];

Clearly there is a point of intersection, a fixed point.

Fixed Points

A fixed point of function g[x], is a point P with the property that g[P]=P.  A fixed point is a point at the intersection of the graphs of y=g[x] and y=x.  It turns out that if a sequence of iterates of g[x], converges, it must converge to a fixed point.  This is shown by the following theorem.

Theorem 2.1  Assume that g is a continuous function and {(Subscript[p, n])}_(n = 0)^∞ is a sequence generated by (Subscript[p, k + 1])= g[(Subscript[p, k])].  If  (Underscript[limit, n -> ∞]) (Subscript[p, n])=P, then P is a fixed point of g[x]

Proof: If (Underscript[limit, n -> ∞]) (Subscript[p, n])=P, then (Underscript[limit, n -> ∞]) (Subscript[p, n + 1])=P so (using continuity of g)
				 g[P]=g[(Underscript[limit, n -> ∞]) (Subscript[p, n])]=(Underscript[limit, n -> ∞]) g[(Subscript[p, n])]=(Underscript[limit, n -> ∞]) (Subscript[p, n + 1])=P

Conditions for Existence of Fixed Points

In the example of (e^(-x)) shown above there is a fixed point on [0,1], and the iteration scheme converges to this fixed point.  This isn't always the case as shown in the following examples.

FixedPointIteration[1 - x, 0.4, 0.1, 30, 
        12, {0, 1}, {0, 1}];

FixedPointIteration[1 - x^2, 0.5, 0.1, 30, 
        12, {-1, 2}, {-1, 2}];

FixedPointIteration[-1 + 3  x, 0.7, 0.1, 30, 
        12, {-10, 10}, {-10, 10}];

FixedPointIteration[1 - 0.25  x^2, 0.9, 0.01, 30, 
        12, {0, 1}, {0, 1}];

What aspect of a function determines whether the sequence of iterates converage to a fixed point?.  As you can see from the examples above, some of the functions have fixed points but the iteration doesn't converge to a fixed point--what additional or other conditions are necessary to have a successful iteration?  Perhaps the initial point needs to be taken 'close enough' to the fixed point--but then, since you may not know the fixed point, this might be difficult.
Theorem 2.2 gives conditions under which g[x] has a fixed point (but doesn't say whether the fixed point iteration scheme will proceed to some fixed point), while Theorem 2.3 gives conditions for the iteration to be successful.

Theorem 2.2:  Assume that g is continuous on [a,b]
 (a) If the range of y=g[x] satisfies yϵ[a,b] for all x ϵ [a,b], then g has a fixed point in [a,b]
 (b) If g'[x] is defined on (a,b) and a positive constant K<1 exists with |g'[x]|≤K<1 for all x ϵ (a,b), then g has a unique fixed point P on [a,b]

Proof: (a) If g[a]=a or g[b]=b, the assertion in (a) is true and you are done.  Otherwise a<g[a] and b>g[b] by the assumption that the range of g[x] lies on [a,b].  The f[x]=x-g[x] satisfies f[a]=a-g[a]<0 and f[b]=b-g[b]>0, so by the intermediate value theorem, there is a point P between a and b such that f[P]=0 or P-g[P]=0, so g[P]=P.
(b) Suppose there are two fixed points, (Subscript[P, 1])and (Subscript[P, 2]),  Then by the mean value theorem, there is a number d ϵ (a,b) such that g'(d)=((g((Subscript[P, 2])) - g((Subscript[P, 1])))/((Subscript[P, 2]) - (Subscript[P, 1])))=(((Subscript[P, 2]) - (Subscript[P, 1]))/((Subscript[P, 2]) - (Subscript[P, 1])))=1.  But this contradicts that hypothesis that |g'(x)|<1 for all x ϵ (a,b).

Theorem 2.3:Assume that g and g' are continuous on[a,b], that K is a positive constant, that the initial point (Subscript[p, 0]) ϵ (a,b), and g[x] ϵ [a,b] for all x ϵ [a,b]
(a) 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
(b) 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: By the mean value theorem, there is a number (Subscript[c, 0])ϵ (a,b) such that |P-(Subscript[p, 1])|=|g(P)-g((Subscript[p, 0]))|=|g'((Subscript[c, 0]))(P-(Subscript[p, 0]))|=|g'((Subscript[c, 0]))| |(P-(Subscript[p, 0]))|≤K|(P-(Subscript[p, 0]))|<|(P-(Subscript[p, 0]))|.  Therefore (Subscript[p, 1]) is no farther from P than (Subscript[p, 0]) was, and so (Subscript[p, 1])ϵ (a,b).  (Look at the examples above where the sequence of approximations is converging to the fixed point to see that this is true).  In general,  |P-(Subscript[p, n])|=|g(P)-g((Subscript[p, n - 1]))|=|g'((Subscript[c, n - 1]))(P-(Subscript[p, n - 1]))|=|g'((Subscript[c, n - 1]))| |(P-(Subscript[p, n - 1]))|≤K|(P-(Subscript[p, n - 1]))|<|(P-(Subscript[p, n - 1]))|.  So, by this last statement and an induction argument, all points (Subscript[p, n]) for n=0,1,2,..., lie in (a,b).  
Next it can be shown by induction, that  |P-(Subscript[p, n])|≤(K^n)|(P-(Subscript[p, 0]))|.  This will imply that as n->∞, the sequence of (Subscript[p, n])'s approaches P, since, by assumption K<1.
In the induction, the case of n=1 is established in the line  |P-(Subscript[p, 1])|=|g(P)-g((Subscript[p, 0]))|=|g'((Subscript[c, 0]))(P-(Subscript[p, 0]))|=|g'((Subscript[c, 0]))| |(P-(Subscript[p, 0]))|≤K|(P-(Subscript[p, 0]))| from above.
So assume |P-(Subscript[p, n - 1])|≤(K^(n - 1))|(P-(Subscript[p, 0]))|.  Then  |P-(Subscript[p, n])|≤K|(P-(Subscript[p, n - 1]))|≤K(K^(n - 1))|(P-(Subscript[p, 0]))|=(K^n)|(P-(Subscript[p, 0]))|.

Error Considerations

How do you stop the iterative process?  As you can see from the FixedPointIteration function defined above, the process stops when 2|(((Subscript[p, n + 1]) -  (Subscript[p, n]))/(small +  ∣ 
                (Subscript[p, n + 1]) ∣ ))|<tolerance, and if this never happens, after the maximum number of iterations allowed.  The author has chosen 2|(((Subscript[p, n + 1]) -  (Subscript[p, n]))/(small +  ∣ 
                (Subscript[p, n + 1]) ∣ ))|as the criterion for smallness but other choices would be |(((Subscript[p, n + 1]) -  (Subscript[p, n]))/(small +  ∣ 
              (Subscript[p, n + 1]) ∣ ))| or |(Subscript[p, n + 1]) -  (Subscript[p, n])|.

Mathematica Solutions of g(x)=x

The command for solving g(x)=x, or other equations is, from the Mathematica Help System:
 Solve[eqns, vars] attempts to solve an equation or set of equations for the variables vars
 You need to be aware that when specifying an equation in Mathematica you must use == instead of equal.  Here is the Mathematica command for finding a fixed point solution for the function (e^(-x)).

Solve[(E^(-x)) == x, x]

Notice a couple of things about the solution: (1) warnings that there are other possible solutions are given--in this case, in which you are dealing with real numbers, you can disregard them.  (2) Mathematica has found an exact solution for the problem--to get a numerical approximation to the solution, see the next statement.  This shows that our approximate answer is close to this approximation.  If you want to see the Mathematica answer with more decimal places use the second of the following statments.

N[ProductLog[1]]

N[ProductLog[1], 20]

Mathematica has another function, NSolve, that immediately provides a numeric approximation to the solution of an equation.

NSolve[(E^(-x)) == x, x]

Another advantage of the Mathematica solve function for g(x)=x, is that it will work in situations where our iterative scheme fails.

Solve[1 - x == x, x]