(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

Calculus Projects

Here we suggest some project directions that may be interesting to pursue.
For subsequent use, we introduce some new functions.

Chap13.m

points/: points[{a_, b_}, {c_, d_}] := Block[{x, y},
         Table[N[{x, y}], {x, a, b, (b - a)/10}, {y, c, d, (d - c)/10}]]
	 

f/: f[x_] := x^4 - 1


iteration/: iteration[f_][x_] := x - f[x]/Derivative[1][f][x]


attractors/: attractors[f_, {a_, b_}, {c_, d_}, n_] := 
             Block[{x, y}, Table[Nest[iteration[f], N[x + y*I], n],
	     {x, a, b, (b - a)/10}, {y, c, d, (d - c)/10}]]
	     
	     
	     
colour/: colour[i_, j_] := Which[
         attractor[[i, j]] == 1., RGBColor[1, 0, 0],
	 attractor[[i, j]] == 1. I, RGBColor[0, 1, 0],
	 attractor[[i, j]] == -1., RGBColor[0, 0, 1],
	 attractor[[i, j]] == -1. I, RGBColor[1, 1, 0],
	 True, RGBColor[1, 0, 1]]


patch/: patch[f_, {a_, b_}, {c_, d_}, n_] := (
        lattice = points[{a, b}, {c, d}];
	attractor = attractors[f, {a, b}, {c, d}, n];
	Show[Table[Graphics[{colour[i, j], Point[lattice[[i, j]]]}],
	{i, 1, 11}, {j, 1, 11}]])

Newton's Method in the Complex Plane

Recall that the formula for Newton's Method is defined as 
iteration[f]  in the package  Chap4.m above. We need this package called in,
and then proceed to investigate the convergence of Newton's method for the
function f(x) =  x^3 - 3 x^2 - x + 3   which has zeros 1, -1, 3 -3.

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

Solve[f[x] == 0, x]

By making use of  iteration[f], define a function   root, as
follows, to perform  n  iterations of Newton's method with initial choice x:

root[f_, x_, n_] := Nest[iteration[f], N[x], n]

root[f, 4.5, 20]

root[f, -4.5, 20]

root[f, 0.3, 20]

Let us now consider the equation x^4 - 1 = 0, which has two real roots, and 
two imaginary roots: 1,-1, i,-i.

f[x_] := x^4 - 1

Solve[f[x] == 0, x]

This time we tabulate our root outputs:

Table[{i, root[f, i, 30]}, {i, -2, 2, 0.1}]

Except for the singularity at the initial choice of 0 (which is not
unexpected since  f'(0) = 0), in each case of the above table,
Newton's method converges to the nearest root. However, things do not stay
so predictable if we allow complex numbers as initial choices. We
shall not cover any calculus of complex numbers here; we only point out
that all the calculus required in Newton's method works equally well for
complex numbers.

To begin with, let us take an initial choice of 1 + 2i, where i is the
square root of -1. Newton's method then converges to the root i:

root[f, 1 + 2 I, 20]

Instead of investigating things further with tables of calculations, let
us use colours. To do this we shall represent the complex number x + y i
by the point  (x, y) in the plane. If the initial choice x + y i
results in Newton's method converging to the root 1 we shall colour in
the point (x, y) red; if the convergence is to the root i we shall
colour the point (x, y) green. Similarly we shall use blue to colour in
initial points which result in convergence to -1, and yellow for initial
points which result in convergence to -i. If there is no convergence, we
shall colour in the initial point magenta.

In the package for this chapter, you will find all the necessary functions to do 
all of the above colouring. In particular, all you need is the function  patch.

?patch

patch[f, {1, 2}, {1, 2}, 30];

Fixed Point Iteration

In this section we extend examples of Chap 4, Section 4, to investigate functional iteration and chaos. We shall need the function   fixedpoint  from Chap4.m:

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[Cos, 0.0, {0, 1}][30];

Consider next the family of parabolic functions  for k>0

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

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

Which homes in on the fixed point quite quickly. However, no such convergence is
obtained if  k=4, even after 100 iterations:

fixedpoint[g[4], 0.1, {0, 1}][100];  (* this will take about 2 min *)

More exercises, discussion and background references can be found in Chapter 13 of 
our book.

Analysing Time-Series Data

Consider a simple time series generated by sine functions

s[t_]:=(Sin[20 2 Pi t/256]+Sin[40 2 Pi t/256]+
Sin[80 2 Pi t/256])/3

This is evidently a sum of three sine waves with frequencies
20, 40 and 80 cycles per second.  Make a table of values of s[t]

st=Table[ N[s[n]], {n,256}];

Plotting we get

ListPlot[st, PlotJoined->True];

Now suppose that we had actually recorded this plot as data from some experiment; 
the usual technique to employ in discovering any periodicities is the Fourier Transform.
We apply this Mathematica function to our data and it gives us  a function of 
frequency which has amplitude proportional to the contribution of that frequency:

fst=ListPlot[Abs[Fourier[st]],
PlotJoined->True, PlotRange->All];

This correctly reveals that there were frequency components at frequencies f=20, 40 and 
80 cycles per second. (Their ghosts appear again at frequencies 256-f where f=20,40,80.
More discussion of the Fourier transform function is given in the Mathematica book.)

Nesting Functions

The function Nest can quickly generate iterations of funtions. For example,
if we  take a quadratic and an initial value , say 0.1, then we can tabulate 
the first three  iterations as follows:

p[x_] := x^2

ptab=Table[ Nest[p, 0.1, n], {n,3}]

Which correctly gives p[0.1], p[p[0.1]], p[p[p[0.1]]].

The following function

f[x_] := 1 - 2 x^2

looks simple enough. It is bounded by +1 and -1 on the interval [-1,1] 
and achieves its single maximum, 1,  at x=0; it has minima of -1 at the
endpoints of the interval and it has two fixed points, at x=-1 and  x=0.5. 
Now pick any number in the interval [-1,1]  and Nest f. For example, 
with initial value again 0.1

Table[ Nest[f, 0.1, n], {n,3}]

which has no obvious pattern. Try a few more iterations

Table[ Nest[f, 0.1, n], {n,10}]

In fact, try a large number and ListPlot them; so give the Table a name first

ftab=Table[ Nest[f, 0.1, n], {n,256}];

plotftab=ListPlot[ftab, PlotJoined->True];

We take the Fourier transform and plot its absolute value, as before

fpt=ListPlot[Abs[Fourier[ftab]], PlotJoined->True];

This time there is no obvious evidence of periodicity, which suggests that there
may be no order in the sequence of values. However, we know from the form of
the function f and the definition of the Nest iteration operation that there is order
present. For if we ListPlot the successive points in the iteration, each against its 
predecessor, we can see the actual structure:

ftabseq = Partition[ftab,2];
ListPlot[ftabseq];

They fall, as is clear from the definition, on the curve defined by the function f that 
is being iterated; order in a very strict sense. Other initial values for the iteration 
give similar results.