#
#	file:   Chaos.mpl
#	a package developed for One-Dimensional Dynamical Systems Lab
#	author: Hinke Osinga, University of Bristol
#	date:   November 2003
#
# 
# 

printf("\n ---  Copyright by Hinke Osinga, 2003  ---\n");
printf(" ---       University of Bristol       ---\n");
printf("Chaos.mpl is a package designed for:\n");
printf("      `One-Dimensional Dynamical Systems'\n");
printf("Course material can be found at:\n");
printf("       http://www.enm.bris.ac.uk/staff/hinke/courses/Chaos/.\n\n");
printf("This package produces pictures of orbits, either by plotting ");
printf("iterates versus time, or by graphical iteration. ");
printf("It can also produce bifurcation diagrams.\n\n\n");


printf("The Logistic Map: ");
printf("f(x) = lambda x (1 - x)\n");

f:= x -> evalf(lambda * x * (1 - x)):
DynSys := x -> evalf(lambda * x * (1 - x));

x0 := 0.94:
lambda := 2:
xrange := 0..1:
lrange := [0, 4]:
critical := 0.5:

printf("lambda = %f, x0 = %f\n\n", lambda, x0);


with(plots):


Iterates := proc(its)
option `Copyright (c) 2003 by Hinke Osinga; all rights reserved.`;
description `Iterates(its,lastits) computes the orbit of x0 up to the n-th iterate, where n = its, and returns a list that contains the last lastits iterates. Iterates(its) uses lastits=10, unless its is smaller`;
      local q, i, fst, lst, orb;
      global x0, DynSys;
      q := x0:
      if (nargs > 1) then lst := args[2]: else lst:= 10: end if;
      if (its < lst) then lst := its: end if;
      fst := its - lst:
      if (fst > 0) then for i to fst do q := DynSys(q) od: end if;
      lst := lst + 1:
      orb := []:
      for i from 1 to lst do orb := [op(orb), q]: q := DynSys(q) od:
      evalm(orb);
end;


CritOrbit := proc(p, m, n)
option `Copyright (c) 2003 by Hinke Osinga; all rights reserved.`;
      local q, i, lst, orb;
      global DynSys;
      q := p:
      for i from 1 to m do q := DynSys(q) od:
      orb := []:
      lst := n-m+1:
      for i from 1 to lst do orb := [op(orb), q]: q := DynSys(q) od:
      evalm(orb);
end;


HigherIts := proc(n, its)
option `Copyright (c) 2003 by Hinke Osinga; all rights reserved.`;
description `HigherIts(n, its) computes the n-th iterate of x0 its times, and returns a list that contains these iterates.`;
      local q, i, orb;
      global x0, DynSys;
      q := x0:
      orb := []:
      for i to its do orb := [op(orb), q]: q := (DynSys@@n)(q) od:
      evalm(orb);
end;


PlotIterates := proc(its)
option `Copyright (c) 2003 by Hinke Osinga; all rights reserved.`;
description `PlotIterates(its,lastits) computes the orbit of x0 up to the n-th iterate, where n = its, and then plots the last lastits iterates versus time. PlotIterates(its) uses lastits=its.`;
      local i, lastits, strt, orb, pts;
      global Iterates, xrange;
      if (nargs > 1) then lastits := args[2]: else lastits:= its: end if:
      if (lastits > its) then lastits:= its: end if:
      orb := Iterates(its,lastits):
      strt := its - lastits - 1:
      pts := [[strt+i, orb[i]] $i=1..lastits+1]:
      plot([pts, pts], x=strt+1..its, y=xrange, color=[red,blue], 
		 style=[point,line], symbol=circle, 
		 tickmarks=[lastits,default],
		 labels=["time","x"], title=`Iterates versus time`);
end;


PlotGraphical := proc(its)
option `Copyright (c) 2003 by Hinke Osinga; all rights reserved.`;
description `PlotGraphical(its) first computes and then shows the iteration of the orbit in a so-called cobweb diagram. The last few iterates are shown in red. The graph of the function (in blue) and the diagonal line y = x are drawn as well.`;
      local orb, pts, lpts, fst, k, gr1, gr2, gr3;
      global Iterates, xrange;
      orb := Iterates(its,its):
      pts := [[orb[1], 0], [orb[1], orb[2]]]:
      lpts := []:
      fst := floor(0.6*its):
      if (fst > 1) then
	 for k from 2 to fst do
	     pts := [op(pts), [orb[k], orb[k]], [orb[k], orb[k+1]]]:
	 od: 
	 fst := fst + 1:
      else
	 fst := 2:
      end if:
      for k from fst to its+1 do
	  lpts := [op(lpts), [orb[k-1], orb[k]], [orb[k], orb[k]]]:
      od:
      gr1 := plot(pts, x=xrange, y=xrange, 
		  color=COLOR(RGB,0.5,0.5,0.5), thickness=1,
	          labels=["x","x"], title=`Cobweb Diagram`):
      gr2 := plot(lpts, x=xrange, y=xrange, color=red, thickness=1,
	  	  labels=["x","x"], title=`Cobweb Diagram`): 
      gr3 := plot([x, 'DynSys'(x)], x=xrange, y=xrange, 
	  	   color=[black,blue],
		   thickness=[1,2], labels=["x","x"], 
		   title=`Cobweb Diagram`):
      display([gr2, gr1, gr3]);
end;



ShowGraph := proc(n, colour)
option `Copyright (c) 2003 by Hinke Osinga; all rights reserved.`;
description `ShowGraph(n, colour) draws the graph of the n-th iterate in the specified colour together with the diagonal line. Choose from colours: red, blue, green, yellow, orange, etcetera, or specify using COLOR(RGB, r, g, b) with r, g, b in [0,1].`;
      global DynSys, xrange;
      plot({x, '(DynSys@@n)'(x)}, x = xrange, xrange,
      color=[black,colour], thickness=[1,2], labels=["x", "x"]);
end;



AddGraph := proc(fig, n, colour)
option `Copyright (c) 2003 by Hinke Osinga; all rights reserved.`;
description `AddGraph(picture, n, colour) draws the graph of the n-th iterate in the specified colour inside the given picture. Choose from colours: red, blue, green, yellow, orange, etcetera, or specify using COLOR(RGB, r, g, b) with r, g, b in [0,1].`;
      local graph;
      global DynSys, xrange;
      graph := plot('(DynSys@@n)'(x), x = xrange, xrange,
      color=colour, thickness=2, labels=["x","x"]):
      display([fig,graph]);
end;


GraphSlope := proc(xstar, colgraph, colslope)
option `Copyright (c) 2003 by Hinke Osinga; all rights reserved.`;
description `GraphSlope(x*, colgraph, colslope) plots the graph of f in the colour "colgraph" and draws the line tangent to this graph at the point (x*, f(x*)) in the colour "colslope". Choose from colours: red, blue, green, yellow, orange, etcetera, or specify using COLOR(RGB, r, g, b) with r, g, b in [0,1].`;
      local graph, grsl, slope, ystar;
      global DynSys, xrange;
      slope := D(DynSys)(xstar):
      ystar := DynSys(xstar):
      printf("Slope at (%f, %f) is %f\n\n", xstar, ystar, slope);
      graph := ShowGraph(1,colgraph):
      grsl := plot(ystar + slope * (x - xstar), x = xrange, xrange,
      color=colslope, thickness=2, labels=["x","x"]):
      display([graph,grsl]);
end;


BifDiagram := proc(transient, actual)
option `Copyright (c) 2003 by Hinke Osinga; all rights reserved.`;
description `BifDiagram(transient, actual) draws an orbit diagram by computing the orbit of an arbitrary point up to 'actual' iterates, but only starts displaying after 'transient' iterates. You can specify other ranges for lambda by calling BifDiagram(transient, actual, [ll, lr]), where the range will then be ll..lr. You can also specify another range for x by calling BifDiagram(transient, actual, [ll, lr], [xl, xr]), with lambda range ll..lr and x range xl..xr.`;
      local lst, orb, pts, ptlst, param, step, stp, i, storelambda, xiv, liv;
      global lambda, critical, CritOrbit;
      if transient < actual then
            storelambda := lambda:
	    xiv := xrange:
	    liv := lrange:
	    if (nargs > 3) then xiv := args[4,1]..args[4,2]: end if;
	    if (nargs > 2) then liv := [args[3,1], args[3,2]]: end if;
	    lst := 200:
	    step := evalf(liv[2] - liv[1])/lst:
	    stp := abs(step) / 2:
	    ptlst := []:
	    param := liv[1]:
	    i := 0;
	    while abs(liv[2] - param) > stp do
		  if mods(i, 10) = 0 then 
		     printf("%d orbits to go\n", lst - i) 
		  end if;
		  lambda := param:
		  orb := CritOrbit(critical, transient,actual);
		  pts := [[param, orb[k]] $k=1..actual-transient+1]:
		  param := param + step:
		  i := i+1:
		  ptlst := [op(ptlst), pts]:
	    od:
	    lambda := storelambda:
	    
	    plot(ptlst, x=liv[1]..liv[2], y=xiv, style=point,
		 color=green, symbol=CROSS, axes=BOXED, 
		 labels=["lambda","x"], title=`Bifurcation Diagram`);
      else
	    printf("You should use a value for transient (= %d) < %d!\n", 
		   transient, actual);
      end if;
end;