Got Questions? Get Answers.
Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Bifurcation Diagram: more matrix operations to make code faster

Subject: Bifurcation Diagram: more matrix operations to make code faster

From: Alexander

Date: 27 Nov, 2009 23:16:04

Message: 1 of 7

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.*(1-yn)

I know no other way than to put this into a loop (here's a simplified version)

for i=1:200
        yn=a.*yn.*(1-yn);
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.*(1-yn);
        if i>ndark
            plot(a,yn,'.','MarkerSize',1)
            drawnow;
        end
    end
end

hold off;

Subject: Bifurcation Diagram: more matrix operations to make code faster

From: Thomas Clark

Date: 27 Nov, 2009 23:54:03

Message: 2 of 7

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.*(1-yn)
>
> I know no other way than to put this into a loop (here's a simplified version)
>
> for i=1:200
> yn=a.*yn.*(1-yn);
> 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.*(1-yn);
> if i>ndark
> plot(a,yn,'.','MarkerSize',1)
> drawnow;
> end
> end
> end
>
> hold off;

Subject: Bifurcation Diagram: more matrix operations to make code faster

From: Alexander

Date: 28 Nov, 2009 08:49:02

Message: 3 of 7

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 linspace-like 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.*(1-yn)
> >
> > I know no other way than to put this into a loop (here's a simplified version)
> >
> > for i=1:200
> > yn=a.*yn.*(1-yn);
> > 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.*(1-yn);
> > if i>ndark
> > plot(a,yn,'.','MarkerSize',1)
> > drawnow;
> > end
> > end
> > end
> >
> > hold off;

Subject: Bifurcation Diagram: more matrix operations to make code faster

From: Thomas Clark

Date: 28 Nov, 2009 14:52:04

Message: 4 of 7

Alexander,

I've rewritten the function in MEX form... but with limited success. It gets quicker, but only where the number of iterations is large - so it could have some use, but possibly not for your immediate application.

I've just submitted it to the file exchange, so when it appears I'll link back to it from here.

For large a, rather than a large number of iterations, the vector operations performed by MATLAB are already highly optimised - MATLAB is only letting us down on the overhead for each iteration. Where this is only ~200 iterations there's very little benefit of using this mex file.

Short of programming on a GPU, I'm running out of ideas (and I'm sorry, I don't know how to do that!).

Kind regards

Tom Clark

Subject: Bifurcation Diagram: more matrix operations to make code faster

From: Rune Allnor

Date: 28 Nov, 2009 15:26:33

Message: 5 of 7

On 28 Nov, 00:16, "Alexander" <alexander.erl...@gmail.com> wrote:

> I know no other way than to put this into a loop (here's a simplified version)
>
> for i=1:200
>         yn=a.*yn.*(1-yn);
> end

There are no other way. The update expression expands
to a 2nd order polynomial in yn, which has to be computed
iteratively.

Rune

Subject: Bifurcation Diagram: more matrix operations to make code faster

From: Lars Barring

Date: 30 Nov, 2009 21:16:20

Message: 6 of 7

"Alexander" <alexander.erlich@gmail.com> wrote in message <heqo5u$6kc$1@fred.mathworks.com>...
> ... 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 linspace-like divide the a's like this: linspace(3.5,4,1000), so the
> whole computation must be repeated).

Try removing the plotting from the loop and store the
result of your calculation instead --- something like this:

---------------------------------
function feigenbaum1(y0, ndark, nsave)
axis([0 4 -0.1 1.1]);

a=[0:0.001:4]';
asav=repmat(a,[1 nsave]);
ysav=repmat(NaN,[length(a) nsave]);
yn=ones(length(a),1).*y0;

for i=1:(ndark+nsave)
  yn=a.*yn.*(1-yn);
  if i>ndark
    ysav(:,i-ndark)=yn;
  end
end
plot(asav(:),ysav(:),'.','MarkerSize',1)
-----------------------------------------------------

Your version times liek this:
>> axis; tic; for k=1:5; feigenbaum(0.5,100,100); drawnow; end; toc
Elapsed time is 38.322195 seconds.

And the version above
>> axis; tic; for k=1:5;feigenbaum1(0.5,100,100);drawnow;end;toc
Elapsed time is 1.942350 seconds.

I.e. something like a 20-fold speedup.

> ... our professor has written a Visual Basic program which is faster
> (and in VB, all you have is loops, matrix operations are not available).

Plotting a large number of individual (simple) objects, like lines and dots,
is rather slow in Matlab compared to VB.


hth
Lars

Subject: Bifurcation Diagram: more matrix operations to make code faster

From: Thomas Clark

Date: 1 Dec, 2009 09:44:19

Message: 7 of 7

Nice one Lars!

Alexander, re. my message above, my file has been rejected from the file exchange (due to having a mex file in it!). It may not be much use to you now, but if you need a quicker way of producing just the end result, drop me an email.
t.clark at cantab.net
... and I'll send it across.

:)

Cheers

Tom

Tags for this Thread

No tags are associated with this thread.

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us