Program 14 - Density for 2-Dimensional Systems

Program Code

Off[General::"spell1", General::"spell"];

Needs["Graphics`Graphics3D`"];

Options[Density] = {
	Type -> DensityPlot
};

RealQ[r_] := NumericQ[r] && Im[r] == 0;

Density::badvars = "Invalid variable list `1` in position 2.";
Density::badlist = "`1` is not a list of pairs of real numbers.";
Density::badtype = "Unknown plot type `1`.";
Density::badbin = "Number of bins `2` for variable `1` is not a positive integer.";
Density::badits = "Number of iterations `1` is not a positive integer.";

DrawDensityPlot[points_, a_, b_, kind_, args___] := 
		Module[{xs, ys, xmin, xmax, ymin, ymax, dx, dy, bins},
	{xs, ys} = Transpose[points];
	
	xmin = Min[xs];
	xmax = Max[xs];
	ymin = Min[ys];
	ymax = Max[ys];
	If[xmin == xmax, xmin -= 1; xmax += 1];
	If[ymin == ymax, ymin -= 1; ymax += 1];
	
	dx=(xmax - xmin) / (a - 1);
	dy=(ymax - ymin) / (b - 1);
	
	bins = Compile[{{xs, _Real, 1}, {ys, _Real, 1}, xmin, ymin, dx, dy, nx, ny, count}, Module[{i, j, k, results},
		results = Table[0, {nx}, {ny}];
		For[k = 1, k <= count, k++,
			i = Round[(xs[[k]] - xmin) / dx + 1];
			j = Round[(ys[[k]] - ymin) / dy + 1];
			If[i <= nx && i >= 1 && j <= ny && j >= 1, results[[i, j]]++];
		];
		results
	]] [xs, ys, xmin, ymin, dx, dy, a, b, Length[points]];
	
	Return[Switch[kind,
		SurfacePlot, ListPlot3D[bins, args, MeshRange->{{xmin,xmax},{ymin,ymax}}],
		BarPlot, BarChart3D[bins, Sequence@@Flatten[{args}]],
		DensityPlot, ListDensityPlot[bins, args, MeshRange->{{xmin,xmax},{ymin,ymax}}]
	]];
];

Density[f_List, vars_List, p0_, n_, a_, b_, args___] := Module[{kind, opts, x, y, g, points},
	If[!(VectorQ[vars, MatchQ[#, _Symbol]&] && Length[vars] == 2),
		Message[Density::badvars, vars]; Return[]];
	If[!IntegerQ[a] || a < 2, Message[Density::badbin, vars[[1]], a]; Return[]];
	If[!IntegerQ[b] || b < 2, Message[Density::badbin, vars[[2]], b]; Return[]];
	If[!IntegerQ[n] || n < 1, Message[Density::badits, n]; Return[]];
	kind = Type /. {args} /. Options[Density];
	If[!MemberQ[{SurfacePlot, BarPlot, DensityPlot}, kind],
		Message[Density::badtype, kind]; Return[]];
	opts = DeleteCases[{args}, Type->_];
	
	g[{x_, y_}] := f /. {vars[[1]] -> x, vars[[2]] -> y};
	
	points = NestList[g, N[p0], n];
	
	Return[DrawDensityPlot[points, a, b, kind, opts]];
];

Density[points_List, a_, b_, args___] := Module[{kind, opts},
	If[!VectorQ[points, MatchQ[#, {_?RealQ, _?RealQ}]&],
		Message[Density::badlist, points]; Return[]];
	If[!IntegerQ[a] || a < 2, Message[Density::badbin, "x", a]; Return[]];
	If[!IntegerQ[b] || b < 2, Message[Density::badbin, "y", b]; Return[]];
	kind = Type /. {args} /. Options[Density];
	If[!MemberQ[{SurfacePlot, BarPlot, DensityPlot}, kind],
		Message[Density::badtype, kind]; Return[]];
	opts = DeleteCases[{args}, Type->_];
	
	Return[DrawDensityPlot[points, a, b, kind, opts]];
];

On[General::"spell1", General::"spell"];

Examples

To use the program, input the functions in a list {f(x,y), g(x,y)} in position 1, then the variables used in position 2, the initial values of these variables in position 3, followed by the number of iterations to use, and then two numbers giving the number of bins to use along the two axes.

Density[{3.6  x (1 - x) - x y, 3.8  y (1 - y) - 0.1  x y}, {
        x, y}, {0.5, 0.1}, 2000, 50, 50, Mesh -> False];

Density[{3.7  x (1 - x), 3.8  y (1 - y)}, {x, y}, {0.5, 0.1}, 
      2000, 50, 50, Mesh -> False];

Density[{3.7  x (1 - x), 3.8  y (1 - y)}, {x, y}, {0.5, 0.1}, 
      2000, 20, 20, Type -> SurfacePlot, Mesh -> False];

Density[{3.7  x (1 - x), 3.8  y (1 - y)}, {x, y}, {0.5, 0.1}, 
      2000, 20, 20, Type -> BarPlot];