@samanthasumo,

I recommend you input your equations into a maple format and copy-paste them into the text window of a comment below. I also recommend that you use the worksheet format rather than the document format and that you use the 1D maple input format rather than 2D format. It will be a lot easier for you to copy-paste into mapleprimes. You can still have the output in the prettier 2D format, if you wish. These are options you set in Tools --> Options. Preferably input everything as an editable format, not as an image (sometimes when you copy-paste Maple worksheets you get an image, may happen if you copy-paste 1-line output, not usually if you copy-paste several lines of Maple output.

@acer 

There may be situations in which one would want to overlay transparent backgrounds on top of one another --- you might have a black background for aesthetic purposes and overlay a transparent background with the (white) words "Copyright me me me", or something like that. In such cases, you wouldn't want the system to remove backgrounds automatically without your permission, so you may want to have an option to disable the automatic reshuffling of the order of the background images. Just a thought in passing.

@acer,

that's a very useful one to know,

ImageTools:-View(ImageTools:-Flip(m,'vertical'));

could easily be added as an option too.

The type declarations I added in my code make it unnecessarily complicated to enter the input. For some reason I thought it was a good thing to have these declarations, but not quite so.

So what's with the image being outside of the Array, can that be fixed with a translation into positive values? (I haven't had time to look at it).

At any rate, there is plenty here for anyone to get started with plotting chaotic maps. These images are prettier than most of those retrieved by an internet search. Great stuff, thanks again to both of you.

Out of curiosity (and for no other reason, since I'm supposed to be very busy with something else right now), I tried a few of the well-known chaotic maps to see if the code could handle them.

http://en.wikipedia.org/wiki/List_of_chaotic_maps

To take 3 examples from the page referenced above:

# the "Tent map" on [1,2]
f := (x,r) -> r*min(x,1-x) ;

# a "Gauss map" on [-1,1]
f := (x,r) -> r+exp(-evalhf(6)*x^2) :

# a "Circle map" (if I understand correctly), can be read by the proc if entered as:
`&modulo`:= proc(x,n)x/n-floor(x/n);end proc:
f := (x,r) -> ( x+evalhf(1/3)-r/evalhf(2*Pi)*sin(evalhf(2*Pi)*x) ) &modulo Pi :

I haven't yet managed to get anything like what we're supposed to get with the Gauss map. I get an image in the top left corner only, with much of the image black (using, say, -0.1 as an initial value and plotting 1000 by 1000 points). I will leave this aside for now, as something to get back to in the future.

@acer, thanks a lot for this implementation!

@epostma, thanks a lot for this improvement!

I can see what you've done with the dummy acer, but I have no understanding of why it is needed; this is really beyond me at this stage.

I thought I would add the option of having a black and white or colored version. As I was playing around with the code, I wanted to substitute the 0 of the rtable for something else (the "background" color). I discovered that
subs[inplace]( 0 = something) , img_h); wouldn't work, but
subs[inplace]( 0.0 = something) , img_h); would work, reflecting the hfloat type of the zero (I guess).

Usage of the variant below is:

  BF( f(x,r)
    , 'range' = [2, 4]
    , 'grid' = [1000, 1000]
    , 'iter' = 5000
    , 'init' = 0.5
    , 'bg' = 0.75 #0 = black, 0.25=yellow,0.5=cyan,0.75=magenta,1=red
    , 'hsv' = [1,5,10] #  shouldn't s be between 0 and 1? or am I getting tired?
    , 'color' = true
  )

where 'color' = false returns a black and white image,
where 'bg' = 0 returns a black background,
where 'bg' > 0 returns different colors according to the hue spectrum yellow, cyan, magenta, red.
Different values for 'bg' and the 'v' of 'hsv' can be combined to obtain a fairly wide range of colors.
other options were described earlier.

Other maps are pretty too, e.g.

f2 := (x,r) -> r*x*(1-x)^3 :
m := BF( f2(x,r), range = [0,11], hsv = [1,1,10] );

yields:

### Version 6. Erik Postma and acer.
### Procedure takes map as input.
### Both color and black and white.
### 2 HSV models currently in competition (fixed V versus V/H=constant)

BF := proc(fexpr
        , { range :: list := [2,4] }
        , { [grid,numpoints] :: list := [1000,1000] }
        , { [iter,iterations] :: posint := 1000 }
        , { [init,initial] :: float := 0.5 }
        , { [bg,background] :: nonnegative := 0 }
        , { [hsv,HSV] :: list := [1,1,1] }
        , { [color,colour] :: boolean := true }
        , $
    )

   description "1). fexpr specifies the map under study, e.g. fexpr:=r*x*(1-x).
   2). range = [ra,rb] specifies the range of values for the bifurcation parameter; optional; default values [2,4].
   3). grid = [xr,yr] specifies the grid of points; optional; default values [1000,1000].
   4). iterations (short form iter) specifies the number of times the map is to be iterated; optional; default value 1000.
   5). initial (short form init) specifies the initial value of the iteration.
   6). background (short form bg) specifies either a black background if bg=0 (default), the hue value of the background if bg>0; optional.
   Rough correspondence for hue bg: 0=black, 0.25=yellow,0.5=cyan,0.75=magenta,1=red. Use in conjunction with v to increase range of colors.
   7). hsv = [h,s,v] controls colors of datapoints,  where h is a SCALE parameter for the hue, s is THE saturation (default range?),
   v is THE value or brightness if bg>0 (between 0 and 1) but a SCALE parameter if bg=0 (any positive value);
   hsv is optional; default values are [1,1,1]
   8). color = false for black and white only; optional; default is true. ";

   local img_h,img_s,img_v,worker,ra,rb,xr,yr,n,h,s,v;
   uses LinearAlgebra, ImageTools, Compiler;

   ra,rb:=op(range): xr,yr:=op(grid): n:=iter:
   h:=HFloat(evalf(hsv[1])): s:=HFloat(evalf(hsv[2])): v:=HFloat(evalf(hsv[3])):

   worker := subs(':-DUMMY'=fexpr,
                  proc(data::Array(datatype=float[8],order=Fortran_order),
                       ra::float[8],rb::float[8],xr::posint,yr::posint,
                       n::posint,init::float[8])
                  local ix::posint,iy::posint,r::float[8],rscale::float[8],
                        x::float[8],y::posint;
                  rscale := (rb - ra) / (xr - 1);
                  for ix to xr do
                     r := ra + rscale * (ix - 1);
                     x := init;
                     for iy to 0 do
                        x := DUMMY;
                     end do;
                     for iy to n do
                        x := DUMMY;
                        y := trunc(x * yr) + 1;
                        if ix>=0 and xr>=ix and y>=0 and yr>=y then
                           data[y, ix] := data[y, ix] + iy^(-0.5); # emphasize early iterations
                        end if;
                     end do;
                  end do;
                  NULL;
               end proc);
   worker := Compile(worker):

   ### Set the "Hue", for h=1 hue varies between 0 and 360
   img_h:=Array(1..xr,1..yr,fill=0,'datatype'='float[8]','order'='Fortran_order');
   try
      worker(img_h,HFloat(evalf(ra)),HFloat(evalf(rb)),xr,yr,n,init);
   catch:
   end try;
   rtable_options(img_h,'subtype'=Matrix);
   MatrixScalarMultiply(img_h,h*360.0/max(img_h),'inplace');
   if bg > 0 then
     subs[inplace]( 0.0 = HFloat(evalf(bg*360)) ,img_h); end if; #took me a while to figure out 0.0
   rtable_options(img_h,'subtype'=Array);

   ### Set the "Saturation"
   img_s:=Array(1..xr,1..yr,fill=s,'datatype'='float[8]','order'='Fortran_order');

   ### Set the "Value", i.e. brightness
   if bg = 0 then
      ### brightness varies with hue (to get black background)
      img_v:=Array(img_h);
      rtable_options(img_v,'subtype'=Matrix);
      MatrixScalarMultiply(img_v,v*1/3.6,'inplace'); #1/3.6 = arbitrary scaling
      rtable_options(img_v,'subtype'=Array);
   else
      ### Fix brightness independently of hue (to get colorful backgrounds)
      img_v:=Array(1..xr,1..yr,fill=v,'datatype'='float[8]','order'='Fortran_order');
   end if;

   ### Combine HSV layers to create colored image or black and white
   if color=true then return HSVtoRGB(CombineLayers(img_h, img_s, img_v)):
   else return RGBtoGray(HSVtoRGB(CombineLayers(img_h, img_s, img_v))): end if:

end proc:

I see, wow thanks, I too have got to rush, will be back at the beginning of next week! have a good weekend. :-)

@acer 

thanks a lot for these explanations, I'm thrilled by this new toy. Thanks a lot. Awesome diagrams!

> So one nice improvement would be to rework it to allow the formula as a parameter.

I thought this should be within my ability, but unfortunately not quite. I'm probably doing something silly, but having attempted to debug for 2 hours, it's high time I turn to the experts for help.

I tried to have the map, e.g. f(x,r)=r*x*(1-x), as an argument. I thought f should be entered as an argument and declared as a procedure, but that didn't work, so I omitted any declaration of type but then (if I understand correctly the error message) the default type (float) is assumed, which doesn't work either.

As I was reading the Maple programming guides, I came across little things like multiple spellings and the trace option, so I have thrown them in too, of limited value added but goes to show I did try hard.

### Version 4. Erik Postma and acer, tweaked by me.

BF8 := proc(
        f ### omitted any type declaration, the plan is to pass f to the "worker" procedure
        , { range :: list := [2.5,5] }
        , { [numpoints,grid] :: list := [1000,1000] }   ### first called it numpoints, then thought perhaps grid was better
        , { [iter,iterations] :: posint := 1000 }
        , { [hsv,HSV] :: list := [1,1,1] }
        , $
    )

    option trace;  ### help in tracking the mess I've made

    local img_h, img_s, img_v, img_hsv, ra, rb, xr, yr, h, s, v, n;

    h:=HFloat(evalf(hsv[1])): s:=HFloat(evalf(hsv[2])): v:=HFloat(evalf(hsv[3])):

    ra,rb:=op(range): xr,yr:=op(numpoints): n:=iter:
   
    ### Hue h
    ### h is a scale parameter for the hue
        img_h := Array(1..xr,1..yr,fill=h,'datatype'='float[8]','order'='Fortran_order');
        worker(img_h,f,HFloat(evalf(ra)),HFloat(evalf(rb)),xr,yr,n); ### Note the second argument f
        rtable_options(img_h,'subtype'=Matrix);
        LinearAlgebra:-MatrixScalarMultiply(img_h,360./h/max(img_h),'inplace');
        rtable_options(img_h,'subtype'=Array);
    ### Saturation s
    ### s is the saturation
        img_s:=Array(1..xr,1..yr,fill=s,'datatype'='float[8]','order'='Fortran_order'):
    ### Value v
    ### v is a scale parameter for the value
        img_v:=Array(1..xr,1..yr,fill=v,'datatype'='float[8]','order'='Fortran_order'):
        # img_v:=Array(img_h);
        # rtable_options(img_v,'subtype'=Matrix);
        # LinearAlgebra:-MatrixScalarMultiply(img_v,10/v,'inplace');
        # rtable_options(img_v,'subtype'=Array);
    ### Hue, Saturation, Value (HSV)
        img_hsv:=ImageTools:-CombineLayers(img_h, img_s, img_v):

    return ImageTools:-HSVtoRGB(img_hsv):

end proc:

worker := proc(
    data::Array(datatype=float[8],order=Fortran_order)
    , f    ### couldn't work out a valid type declaration
    , ra::float[8]
    , rb::float[8]
    , xr::posint
    , yr::posint
    , n::posint
)

local
    ix::posint
    , iy::posint
    , r::float[8]
    , rscale::float[8]
    , x::float[8]
    , y::posint
;
 
rscale := (rb - ra) / (xr - 1);

for ix to xr do
    r := ra + rscale * (ix - 1);
    x := HFloat(0.5);
    for iy to 3 do  # Warmup:
        x := f(x,r);
    end do;
    for iy to n do
        x := f(x,r);
        y := trunc(x * yr) + 1;
        data[y, ix] := data[y, ix] + 1;
    end do;
end do;

NULL;

end proc:

worker := Compiler:-Compile(worker):

# Test functions
f := proc(x,r) r*x*(1-x); end proc;
f := (x,r) -> r*x*(1-x);

m := CodeTools:-Usage( BF8( f(x,r), range = [2.5, 5], numpoints = [1200, 1200], iter = 50000, hsv = [0.008,0.8,0.4] ) ):
here := cat(currentdir(),kernelopts(dirsep),"BifurcationDiagramInColorV4.tiff"):
ImageTools:-Write(here, m, 'format' = 'TIFF') : # 3 valid formats: TIFF, JPEG, BMP
{--> enter BF8, args = r*x*(1-x), range = [2.5, 5], numpoints = [1200, 1200], iter = 50000, hsv = [0.8e-2, .8, .4]
Error, (in worker) invalid input: expecting the 2nd argument to be of type float, but received r*x*(1-x)
Error, invalid input: ImageTools:-Write expects its 2nd argument, image, to be of type ImageTools:-Image, but received m

@Markiyan Hirnyk > The notion *N is used in nonstandard analysis.

thanks, that's not something I know, have followed up on the link :-)

in my recollection, the notation would be (a, b) \in N* X N, with a cross between N* and N, rather than the dot the OP used. Anyways, I think the OP is happy now, thanks!

@Markiyan Hirnyk > The notion *N is used in nonstandard analysis.

thanks, that's not something I know, have followed up on the link :-)

in my recollection, the notation would be (a, b) \in N* X N, with a cross between N* and N, rather than the dot the OP used. Anyways, I think the OP is happy now, thanks!

@acer and @epostma

This is excellent, I just spent several hours playing around with this, trying to understand the code and how to manipulate colours. I found it convenient to add an hsv option to the main procedure, to control the colours. Below is a sample of some of the prettiest (?) pictures I got.

Several questions popped into my head, not particularly pressing, consigned here for reference, in random order:
1. what is the "warm-up" for?
2. I wasn't entirely sure of where the scaling of Hue and Value came from, namely 360.0/max(img_h) and 1/3.6
3. I tried to shorten the code by removing the transformation to a Matrix and the call to LinearAlgebra, for instance by writing:
map[evalhf,inplace](x -> 360.0/HFloat(evalf(max(img_h)))*x, img_h);
but this returns the following error:
Error, (in BF8) lexical scoping is not supported in evalhf
Is that why the transformation to a Matrix is necessary?

In the following I have three new parameters, h, s, v to control the colours.

BF8 := proc(ra,rb,xr,yr,n,{hsv::list:=[1,1,1]})
    local img_h, img_s, img_v, img_hsv, h, s, v;
    h:=HFloat(evalf(hsv[1])): s:=HFloat(evalf(hsv[2])): v:=HFloat(evalf(hsv[3])):

    ### Hue h
    ### h is a scale parameter for the hue
        img_h := Array(1..xr,1..yr,fill=h,'datatype'='float[8]','order'='Fortran_order');
        worker(img_h,HFloat(evalf(ra)), HFloat(evalf(rb)),xr,yr,n);
        rtable_options(img_h,'subtype'=Matrix);
        LinearAlgebra:-MatrixScalarMultiply(img_h,360/h/max(img_h),'inplace');
        rtable_options(img_h,'subtype'=Array);
    ### Saturation s
    ### s is the saturation
        img_s:=Array(1..xr,1..yr,fill=s,'datatype'='float[8]','order'='Fortran_order'):
    ### Value v
    ### v is the value
        img_v:=Array(1..xr,1..yr,fill=v,'datatype'='float[8]','order'='Fortran_order'):
    ### Hue, Saturation, Value (HSV)
        img_hsv:=ImageTools:-CombineLayers(img_h, img_s, img_v):

    return ImageTools:-HSVtoRGB(img_hsv):

end proc:

[rest of the code unchanged]

and then I would use it as follows:

m := CodeTools:-Usage( BF8(2.6, 4, 1200, 1200, 50000, hsv = [10,0.9,0.9] ) ):
here := cat(currentdir(),kernelopts(dirsep),"BifurcationDiagramInColorV3b.tiff"):
ImageTools:-Write(here, m, 'format' = 'TIFF') : # 3 formats: TIFF, JPEG, BMP

where the list hsv = [10,0.9,0.9] helps control the colours somewhat.

Actually, if wikipedia may be trusted,

There is no universal agreement about whether to include zero in the setof natural numbers: some define the natural numbers to be the positive integers{1, 2, 3, ...}, while for others the term designates the non-negative integers {0, 1, 2, 3, ...}.

Some authors use the term "natural number" to exclude zero and "whole number" to include it; others use "whole number" in a way that excludes zero, or in a way that includes both zero and the negative integers.

further down the page, there is confirmation of the notation N*=N\{0}.

http://en.wikipedia.org/wiki/Natural_number

Actually, if wikipedia may be trusted,

There is no universal agreement about whether to include zero in the setof natural numbers: some define the natural numbers to be the positive integers{1, 2, 3, ...}, while for others the term designates the non-negative integers {0, 1, 2, 3, ...}.

Some authors use the term "natural number" to exclude zero and "whole number" to include it; others use "whole number" in a way that excludes zero, or in a way that includes both zero and the negative integers.

further down the page, there is confirmation of the notation N*=N\{0}.

http://en.wikipedia.org/wiki/Natural_number

I don't have Maple right now to check, but the problem may be that the bvp solver is less flexible than the ivp solve, you should check the maple help files (http://www.maplesoft.com/support/help/), perhaps the problem is that the bvp solver cannot take a parameter as an option. If you look inside Preben's procedure, you'll see that while the proc defines the solution in terms of the unknown parameter, the computation is done after a value is assigned to the parameter; so depending on your particular purpose it may be equally easy to simply set a:=1; (or whatever) at the start of your worksheet, or like I did use subs to replace the parameter by a value of your choice (this could be done in a loop with the seq command if you wanted to see how the solution depends on the parameter); or you could use Preben's procedure which does it in a more elegant and efficient way. What I meant about the boundary conditions is that if you actually know the slope of the solution at some point (for instance, the solution may approach 0 along the horizontal axis, with a slope of zero) then you may be able to transform the bvp into an ivp, which is easier to solve.

I don't have Maple right now to check, but the problem may be that the bvp solver is less flexible than the ivp solve, you should check the maple help files (http://www.maplesoft.com/support/help/), perhaps the problem is that the bvp solver cannot take a parameter as an option. If you look inside Preben's procedure, you'll see that while the proc defines the solution in terms of the unknown parameter, the computation is done after a value is assigned to the parameter; so depending on your particular purpose it may be equally easy to simply set a:=1; (or whatever) at the start of your worksheet, or like I did use subs to replace the parameter by a value of your choice (this could be done in a loop with the seq command if you wanted to see how the solution depends on the parameter); or you could use Preben's procedure which does it in a more elegant and efficient way. What I meant about the boundary conditions is that if you actually know the slope of the solution at some point (for instance, the solution may approach 0 along the horizontal axis, with a slope of zero) then you may be able to transform the bvp into an ivp, which is easier to solve.

@epostma

wow, so elegant and, as acer points out, so fast. A big thank you.

P.S. I'll need some time to fully understand your code, right now I just marvel.