(C) Copyright 1994 D.C.M. Burbulla and  C.T.J. Dodson
       Notebook to support the book
     Self-Tutor for Computer  Calculus

Using Mathematica 2

by
     D. C. M. Burbulla and C. T. J. Dodson, University of Toronto
 Published by:    
    Prentice-Hall Canada  (English)   ISBN 0130152803     and 
     Prentice-Hall Japan  (Japanese) ISBN4-8101-8558-3 C3055 P2880E

Methods of Approximation

Chap4.m

bisect/: bisect[f_][{a_, b_}] :=
         If[N[f[a]*f[b]] > 0, {a, b},
	 If[N[f[a]*f[(a + b)/2]] <= 0, N[{a, (a + b)/2}], N[{(a + b)/2, b}]]] 
	 
bisect::usage = "bisect[f][{a, b}] determines which half of the interval {a, b} contains 
a solution to the equation f[x] = 0. It is assumed that f is continuous on {a, b}. 
If f[a] and f[b] have the same sign, then the output is the original interval."


tangent/: tangent[f_, a_][x_] := f'[a]*(x - a) + f[a]


normal/: normal[f_, a_][x_] := (a - x)/f'[a] + f[a]


tangentplot/: tangentplot[f_, a_, {b_, c_}] :=
              Plot[{f[x], tangent[f, a][x]}, {x, b, c}] 
	      

normalplot/: normalplot[f_, a_, {b_, c_}] :=
             Plot[{f[x], normal[f, a][x]}, {x, b, c}, AspectRatio -> Automatic] 


errorplot/: errorplot[f_, a_, {b_, c_}] :=
            Plot[Abs[f[x] - tangent[f, a][x]], {x, b, c}, PlotRange -> All] 
	    
errorplot::usage = "errorplot[f, a, {b, c}] plots the absolute value of the difference of 
f and the tangent to f at x = a, for each x in the interval {b, c}."


iteration/: iteration[f_][x_] := x - f[x]/f'[x]


point/: point[f_, a_, n_] := Nest[iteration[f], N[a], n]


newtonsmethod/: newtonsmethod[f_, a_, {b_, c_}][n_] := (one = 
                Block[{x, i}, Plot[{f[x], tangent[f, point[f, a, n]][x]},
		{x, First[Union[{b}, Table[point[f, a, i], {i, 0, n + 1}]]],
		Last[Union[{c}, Table[point[f, a, i], {i, 0, n + 1}]]]},
	        DisplayFunction -> Identity]];
		two = Block[{i}, ListPlot[Table[{point[f, a, i], 0},
		{i, 0, n + 1}], DisplayFunction -> Identity]];
		three = Block[{i}, ListPlot[Table[{point[f, a, i], 
		f[point[f, a, i]]}, {i, 0, n}], DisplayFunction -> Identity]];
		Show[one, two, three, PlotLabel -> 
         If[n == 0, If[a == 0, "0 = initial guess", If[a == 1,
        "1 = initial guess", a " = initial guess"]],
           {(n + 1) " iterations completed",
            N[Nest[iteration[f], a, n + 1]] " = approximation"} ],
              DisplayFunction -> $DisplayFunction])

		
newtonsmethod::usage = "newtonsmethod[f, a, {b, c}][n] shows Newton's method being 
used to generate the (n + 1)st approximation to the solution of the equation f[x] = 0 from 
the nth approximation, starting with initial guess x0 = a. All previous approximations are 
shown along the x-axis, and the value of the current approximation is indicated in the 
plot label. The interval {b, c} is the initial x-interval used for the plot; however the 
x-interval will adjust automatically, if necessary, to include all approximations 
x0, x1, ... , xn, x(n + 1)."


fixedpoint/: fixedpoint[g_, a_, {b_, c_}][n_] := (one = 
             Block[{x, i}, Plot[{g[x], x}, {x, b, c}, 
	     AspectRatio -> Automatic, PlotRange -> {{b, c}, {b, c}}, 
	     DisplayFunction -> Identity]];
	     two = Block[{i}, ListPlot[Table[{Nest[g, N[a], i], 0}, 
	     {i, 0, n + 1}], DisplayFunction -> Identity]];
	     three = Block[{i}, ListPlot[Table[{Nest[g, N[a], i], 
	     Nest[g, N[a], i]}, {i, 0, n + 1}], DisplayFunction -> Identity]];
	     four = Block[{i}, ListPlot[Table[{Nest[g, N[a], i], 
	     Nest[g, N[a], i + 1]}, {i, 0, n + 1}],
	     DisplayFunction -> Identity]];
	     five = Block[{i}, ListPlot[Table[{Nest[g, N[a], Floor[i/2]],
	     Mod[i + 1, 2]*Nest[g, N[a], Floor[i/2]] + 
	     Mod[i, 2]*g[Nest[g, N[a], Floor[i/2]]]}, {i, 0, 2*n + 3}],
	     PlotJoined -> True, DisplayFunction -> Identity]];
	     Show[one, two, three, four, five,
        PlotLabel -> If[n == 0,
             {"beginning fixed point iteration",
        If[a == 0, "0 = initial guess", If[a == 1,
        "1 = initial guess", a " = initial guess"]]},
   {(n + 1) " iterations completed", N[Nest[g, a, n + 1]] 
" = approximation"}],
                DisplayFunction -> $DisplayFunction])

	     

fixedpoint::usage = "fixedpoint[f, a, {b, c}][n] performs (n + 1) iterations 
of the Fixed Point iteration method for approximating solutions to the equation 
f[x] = x, starting with x0 = a. The output plots f and the identity function, y = x, 
on the interval {b, c}; shows all successive approximations x0, x1, ... , xn, x(n + 1) 
along the x-axis; indicates the value of the current approximation in the PlotLabel; 
and shows the `cobweb' of points homing in on the fixed point."

Linear Approximation

tangent[Sqrt, 4][x]

Simplify[%]

tangent[Sqrt, 4][{7.2, 2.5, 4.8, 3.95}] // N

Sqrt[{7.2, 2.5, 4.8, 3.95}] // N

We can see the error in tangent approximations using the new function errorplot:

errorplot[Sqrt, 4, {0, 8}];

errorplot[Sqrt, 4, {2, 6}];

The Bisection Method

This cubic cannot easily be factored so finding exact solutions is incovenient, but 
you can see from the graph that there are three solutions to f(x)=0:

Clear[f]; f[x_] := x^3 + 4 x^2 - 6

Plot[f[x], {x, -4, 3}];

Plot[f[x], {x, -2, -1}];

Plot[f[x], {x, -2, -1.5}];

Plot[f[x], {x, -1.75, -1.5}];

Plot[f[x], {x, -1.625, -1.5}];

Now try the new bisect function:

?bisect

bisect[f][{-2, -1}] (* f as before *)

bisect[f][%]

bisect[f][%]

bisect[f][%]

More neatly, we can use Nest:

?Nest

Nest[bisect[f], {-2, -1}, 4]

Nest[bisect[f], {-2, -1}, 18]

Nest[bisect[f], {-4, -3}, 19]

Nest[bisect[f], {1, 2}, 18]

You can try also NRoots:

?NRoots

NRoots[x^3 + 4 x^2 - 6 == 0, x]

Newton's Method

Newton's method requires an initial choice, and the sequence may not
converge. However, when Newton's method does work, it works so well 
that the convergence is usually much quicker than that of the bisection method.

Example 1

Clear[f]; f[x_] := x^3 + 4 x^2 - 6

Plot[f[x], {x, -4, 3}];

?newtonsmethod

newtonsmethod[f, 0.4, {0, 3}][0];

newtonsmethod[f, 0.4, {0, 3}][1];

newtonsmethod[f, 0.4, {0, 3}][2];

newtonsmethod[f, 0.4, {0, 3}][3];

newtonsmethod[f,  0.4,  {0, 3}][4];

Note that each graph produced with  newtonsmethod  keeps track of all
the previous approximations, and gives you the value of the current
approximation. Of course, if the initial guess were closer to the root
then the convergence would be quicker.

Example 2

?iteration

Nest[iteration[f], 0.4, 6]  (* f as before  *)

Nest[iteration[f], -1.5, 4]

Nest[iteration[f], -3.5, 5]

?FindRoot

FindRoot[f[x] == 0, {x, 1.5}]

FindRoot[f[x] == 0, {x, -1.5}]

FindRoot[f[x] == 0, {x, -3.5}]

FindRoot[Sin[x] == x, {x, 12}]

Functional Iteration and Fixed Points

Finding a fixed point of a function g(x) is  equivalent to solving the equation
g(x)  = x. The Brouwer Fixed Point Theorem says that every continuous 
function from a closed n-ball to itself (in particular, from a a closed interval
to itself)  has a fixed point. Finding it is an interesting exercise in successive
approximations. The Mathematica function FixedPoint provides one method,
for example, we estimate a fixed point of the cosine function starting form an 
initial guess of 0.0.

FixedPoint[Cos,0.0,20]

Cos[%]

After a few more iterations the procedure homes in on the fixed point at  0.739085.
In our book we  defined a function, fixedpoint, to illustrate graphically the progress of
functional iteration to find a fixed point; we look at this and some other functions.

Example 1

FindRoot[Cos[x] == x, {x, 0.0}]

Nest[Cos, 0.0, 35]

?fixedpoint

fixedpoint[Cos, 0.0, {-0.5, 2}][0];

fixedpoint[Cos, 0.0, {-0.5, 2}][1];

fixedpoint[Cos, 0.0, {-0.5, 2}][2];

fixedpoint[Cos, 0.0, {-0.5, 2}][6];

?FixedPoint

FixedPoint[Cos, 0.0]

?FixedPointList

FixedPointList[Cos, 0.0, 40]

ListPlot[%];

Example 2 
Consider next the family of parabolic functions  for k>0:

Clear[g]; g[k_][x_] := k x (1 - x)

FixedPoint[g[2.6], 0.1]

fixedpoint[g[2.6], 0.1, {0, 1}][20];

ListPlot[ FixedPointList[g[2.6], 0.1, 20] ];

Example 3

fixedpoint[g[3.2], 0.1, {0, 1}][50];

ListPlot[ FixedPointList[g[3.2], 0.1, 50] ];

fixedpoint[g[3.561], 0.1, {0, 1}][50];

ListPlot[ FixedPointList[g[3.561], 0.1, 50] ];