####################################################################
#
# File:   Iter
# Author: J. Carlson
# Date:   August 9, 1993
# Last revised: October 2, 1993
#
# "Iter" is a package of simple tools for studying iteration of mappings. 
# The main tools are the functions are "iter" and "iterplot", "iterplot2",
# and "cobweb". 
#
# EXAMPLES:
#
# If we say
#
#    double := x -> 2*x;   # define doubling function
#
# then
#
#    results := iter( double, 4, 1 ); print(results);
#
# yields
#
#    [ 1, 2, 4, 8, 16 ]
#
# Note that the output is [ 1, f(1), f(f(1)), f(f(f(1))), f(f(f(f(1)))) ]
# It is produced by applying f to the initial value 1 a total of 4 times.
# In general, if one says
#
#    y := iter( f, n, a )
# 
# then the result of iteratively applying   f   to the initial value   a
# is stored in   y  .   To view   y    say
# 
#    print(y);
#
# or
#
#    printarray(y);
#
# The function printarray is defined below.  You can, of course,
# say things like   printarray(iter( f, n, a )) , but then the result
# of the computation is not saved in any variable.  A shortcut is
#
#    iterprint( f, n , a );
# 
# To produce graphical output, use
#
#   iterplot( f, n, a );
#
# Thus 
#
#   iterplot( double, 4, 1 );
#
# displays a graph of the data [ 1, 2, 4, 8, 16 ] computed above 
# using iter.  Here is one last example:
#
#    quad := x -> k*x*(1-x);
#    k := 3.0;
#    iterplot( quad, 40, 0.1 );
#    k := 4.0;
#    iterplot( quad, 40, 0.1 );
#
#
# The function "iterplot2" is like iterplot, but iterates
# two functions instead of one, as in the example below:
#
#    iterplot2( quad, quad, 40, 0.1, 0.9 );
#
# The syntax is
#
#    iterplot2 :=  proc(f,g,n,a,b)
#
# where a and b are the initial values used with f and g.
#
#
#
#
# The function "cobweb" is useful for demonstrating graphicaly
# how a map is iterated.   It is also useful for for finding
# fixed points.  The syntax is 
#
#   cobweb( f, n, s, a..b );
#
# where f is the function iterated, n is the number of iterations,
# s is the starting (initial) value, and a..b is the interval
# used to make the graph.
#
# Example:
#
#   k := 2.5; cobweb( quad, 10, 0.1, 0..1 );
#
# This example "cobwebs" quad 10 times on the interval 0..1 with
# an initial value of 0.1.
#

####################################################################


####################################################################
#
# Definitions:
#
# iter( f, n, a) returns (n+1)-element array [a, f(a), f(f(a)), ... ]
#
# iterprint(f, n, a) print the array [a, f(a), f(f(a)), .. ]
#
# iterplot is like iter, but plots the resulting sequence of iterates 
#
# iterplotna is like iter, but plots the resulting sequence of iterates 
#   for a non-autonomous function
#
# iterplot2 is like iterplot, but plots the result of iterating two
# functions.  The syntax is
#
#    iterplot2(f, g, n, a, b)
#
# where a and b are the initial values for the sequences produced
# by f and g, respectively.
#
# iter2 is used by iterplot
#
# iter2na is used by iterplotna
#
# printarray is used to display results of iter and is used in 
# iterprint
#
#####################################################################


printarray := proc( L )
  local i;
  for i from 0 to op(2,op(2,eval(L))) do 
    print(i, L[i])
   od;
end:


iter :=  proc(f,n,a) 
  local x, i;
  x := array(0..n);
  x[0] := a;
  for i from 1 to n do 
    x[i] := f(x[i-1]) 
  od;
  RETURN(x);
end:

iterprint := proc(f,n,a)
  local y, z;
    y := iter(f,n,a);
    printarray(y);
end:
		 
iterplot := proc(f,n,a)
  local y, z;
  y := iter2(f,n,a);
  z := convert(y,list);
  plot(z, 0..n, style=LINE);
end:

iterplotna := proc(f,n,a)
  local y, z;
  y := iter2na(f,n,a);
  z := convert(y,list);
  plot(z, 0..n, style=LINE);
end:

iterplot2 :=  proc(f,g,n,a,b)
  local y, z, y2, z2;
  y := iter2(f,n,a);
  z := convert(y,list);
  y2 := iter2(g,n,b);
  z2 := convert(y2,list);
  plot({ z, z2 }, 0..n, style=LINE);
end:

iterplotP := proc(f,n,a)
  local y, u, z, i, j;
  u := array(0..2*n);
  y := iter2(f,n,a);
  j := 0;
  for i from 0 to n-1 do
    u[j] := [y[i][2],y[i][1]] ;
    j := j+1;
    u[j] := [y[i+1][2],y[i][1]] ;
    j := j+1;
  od;
  u[2*n] := [y[n][2],y[n][1]] ;
  z := convert(u,list);
  plot(z, 0..y[n][2], style=LINE);
end:
		 
iter2 := proc(f,n,a) 
  local x, i;
  x := array(0..n);
  x[0] := [0, a];
  for i from 1 to n do 
    x[i] := [ i, f(x[i-1][2]) ] 
  od;
  RETURN(x);
end:

iter2na := proc(f,n,a) 
  local x, i;
  x := array(0..n);
  x[0] := [0, a];
  for i from 1 to n do 
    x[i] := [ i, f(i,x[i-1][2]) ] 
  od;
  RETURN(x);
end:

makethread := proc( f, n, a )
   local T, i, x, y;
   T := array(0..n);
   x := evalf(a);
   y := 0.0;
   T[0] := [x,y];
   for i from 1 by 2 to n do
      x := T[i-1][1];
      y := evalf( f(x) );
      T[i] := [ x, y ];
      T[i+1] := [ y, y ];
   od;
   convert( T, list );
end;

discretize := proc( f, n, a, b)
  local T, i, x, y, h;
  T := array(0..n);
  x := evalf(a);
  y := f(x);
  h := evalf( (b-a)/n );
  T[0] := [x,y];
  for i from 1 to n do
    x := x + h;
    y := f(x);
    T[i] := [x,y];
  od;
  convert( T, list );
end;


xminoflist := proc(L)
  local x, xm;
    xm:= L[1][1];
    for i from 2 to nops(L) do
       x := L[i][1];
      if (x < xm) then xm := x fi;
    od;
  xm;
end;

xmaxoflist := proc(L)
  local x, xm;
    xm:= L[1][1];
    for i from 2 to nops(L) do
       x := L[i][1];
      if (x > xm) then xm := x fi;
    od;
  xm;
end;

yminoflist := proc(L)
  local y, ym;
    ym:= L[1][2];
    for i from 2 to nops(L) do
       y := L[i][1];
      if (y < ym) then ym := y fi;
    od;
  ym;
end;

ymaxoflist := proc(L)
  local y, ym;
    ym:= L[1][2];
    for i from 2 to nops(L) do
       y := L[i][1];
      if (y > ym) then ym := y fi;
    od;
  ym;
end;

id := x -> x; # identity function # 

cobweb := proc( f, n, a, xrange )
          local thread1, thread2, thread3, epsilon, x1, x2, y1, y2, yrange;
          yrange := xrange;
          epsilon := 0.1;
	          thread1 := makethread( f, 2*n, a );
          x1 := op(1,xrange); x2 := op(2,xrange);
          y1 := op(1,yrange); y2 := op(2,yrange);

          xmin := xminoflist(thread1);
          xmax := xmaxoflist(thread1);
          ymin := yminoflist(thread1);
          ymax := ymaxoflist(thread1);
          if (x1 > xmin) then x1 := xmin fi;
          if (x2 < xmax) then x2 := xmax fi;
          if (y1 > ymin) then y1 := ymin fi;
          if (y2 < ymax) then y2 := ymax fi;
          x1 := x1 - epsilon*(x2-x1);
          x2 := x2 + epsilon*(x2-x1);

          thread2 := discretize(f, 20, x1, x2 );
          thread3 := discretize(id, 1, x1, x2 );
          plot( { thread1, thread2, thread3 },
                x = x1..x2, y=y1..y2, style = LINE );
    end;



## polygon for testing:
#
# square := [ [0,0], [0,1], [1,1], [1,0], [0,0] ];

## functioin for testing:
#
# f := x -> k*x*(1-x); k := 2;

## test:
#
# print(`Fred, try this: cobweb( f, 10, 0.1, 0..1 );`);

######################################################################