# Iteration of mappings

# Jim Carlson 
# http://www.math.utah.edu/~carlson)
# August 9, 1993
# Last revised: October 2, 1993


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=POINT);
end:

iterplotna := proc(f,n,a)
  local y, z;
  y := iter2na(f,n,a);
  z := convert(y,list);
  plot(z, 0..n, style=POINT);
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=POINT);
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 i,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 i,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 i,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 i,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,
                 xmin,xmax,ymin,ymax;
                 
          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: