Hi Thomas,
thanks for your reply. If I am not mistaken, this kind of diagram was first created by Robert May (created only in parts as he lacked computing power), a theoretical Biologist. He was investigating the logistic map, which is he regarded as a very simplified model for population growth. He also chose a from 1 to 4.
The interesting thing (and the reason why I care about speeding up the algorithm) is that there is a meaning the period doubling (the ratio is Feigenbaum's constant 4,66920...). Feigenbaum found the same ratio even in totally different maps containing e.g. trigonometric functions. I would like to rewrite this program for other maps, too, just to see that indeed the ratio is preserved in other maps. Also a nice way to study is to make a zoom box function using rbbox. But for this purpose, whenever you zoom, it is necessary to redo all the computations with different a's (like if you zoom from a=3.5 to a=4, it is necessary to linspacelike divide the a's like this: linspace(3.5,4,1000), so the whole computation must be repeated). This shouldn't take as long as it does, as our professor has written a Visual Basic program which is faster (and in VB, all you have is loops, matrix operations
are not available).
Regards,
Alexander
"Thomas Clark" <t.clark@remove.spamcantab.net> wrote in message <hepoqr$4kr$1@fred.mathworks.com>...
> Hi Alexander,
>
> That's a really interesting algorithm. Do you mind if I ask what the application is?
>
> In terms of speeding it up, I'm afraid i tried to embed it using emlmex, which usually speeds these things up  but I actually got a factor of two slower!
>
> Possibly writing a short FORTRAN mex file might help. Unfortunately, the algorithm is deterministic  it isn't possible to vectorise the outer loop.
>
> What size of a are we talking about for your practical purposes?
>
> Kind Regards
>
> Tom Clark
>
>
>
>
> "Alexander" <alexander.erlich@gmail.com> wrote in message <hepmjk$jce$1@fred.mathworks.com>...
> > Hello,
> >
> > I've written a small function which creates a bifurcation diagram (like this one: http://en.wikipedia.org/wiki/File:LogisticMap_BifurcationDiagram.png). The program is very short, I've posted it below. It uses the logistic map.
> >
> > In the whole concept, there are only two loops: the outer one varies the a parameter which regulates how "chaotic" the behaviour of the system is. The inner loop calculates a couple of iterations with a given a value and plots them (see program comment for more on this if you like).
> >
> > I managed to vectorize the a variation (putting all a's into a large vector). What remains is this iteration: yn=a.*yn.*(1yn)
> >
> > I know no other way than to put this into a loop (here's a simplified version)
> >
> > for i=1:200
> > yn=a.*yn.*(1yn);
> > end
> >
> > As a is very large, this computation must be crutially slow. I wish there was a way to use fast Matlab routines for this (instead of a loop). Could you think of a way to vectorize this code?
> >
> > The entire program is given below. Here's an example call:
> > Feigenbaum(0.5,100,100);
> >
> > I'd be glad if you could hint me on a way to make this code faster.
> >
> > Regards,
> > Alexander
> >
> > Feigenbaum.m
> > function Feigenbaum(y0, ndark, nsave)
> > %Feigenbaum Draws a Feigenbaum bifurcation diagram for the logistic
> > %equation.
> > % y0: initial value of the logistic equation iteration
> > % ndark: darkroom iterations. These iterations are not plotted. They are
> > % required to get the current iterated value close to the current
> > % attractor.
> > % nsave: Once ndark calculations have been computed, an additional
> > % number (nsave) of calculations is performed. These ones are
> > % actually plotted. Thus, with a given a value, ndark+nsave
> > % calculations are performed.
> > %
> > % example call: Feigenbaum(0.5,100,100);
> >
> > h=plot(0,0); hold on;
> >
> > axis([0 4 0.1 1.1]);
> >
> > a=[0:0.001:4]';
> > yn=ones(length(a),1).*y0;
> >
> > for i=1:(ndark+nsave)
> > if ishandle(h)==1
> > yn=a.*yn.*(1yn);
> > if i>ndark
> > plot(a,yn,'.','MarkerSize',1)
> > drawnow;
> > end
> > end
> > end
> >
> > hold off;
