1.7 Bifurcations

September 16, 1998

How to draw a bifurcation diagram

We shall consider the one-parameter family of differential equations dy/dt = y^6-2*y^4+alpha .  The main curve in a bifurcation diagram is the curve consisting of equilibrium points for different values of alpha .  In other words, we simply want to plot the curve y^6-2*y^4+alpha = 0 .  In many simple cases, this is easy to solve for the parameter, as it is here : alpha = -y^6+2*y^4 . This can then be plotted.

> plot(-y^6+2*y^4,y=-5..5);

[Plot]

One of the difficulties with this is that the axes are interchanged -- we want the y -axis to be the vertical one, not the horizontal one.  Another problem that can arise is that it may not be especially easy to solve f(y, alpha) = 0 for the parameter alpha .  We can take care of this using a different plotting command.

> with(plots):
implicitplot(y^6-2*y^4+alpha=0,alpha=-5..5,y=-5..5);

[Plot]

Note that this has managed to pick up some behavior missing in the first plot, where alpha ranged from 0 to -14000.  It has a rather jagged look to it, though.  This is because Maple is not plotting enough points.  That can be fixed using the plot option numpoints.

> implicitplot(y^6-2*y^4+alpha=0,alpha=-5..5,y=-5..5,numpoints=20000);

[Plot]

This, of course, does not give the whole bifurcation diagram -- we need some vertical lines with arrows on them.  However, it is a good first step.  The curve pictured divides the plane into two regions, one including the negative alpha -axis, the other including the positive alpha -axis.  In the first, the righthand side y^6-2*y^4+alpha is negative and the arrows would point down.  In the second, the righthand side is positive and the arrows would point up.  There seem to be three bifurcation points.  We can get Maple to find them.

> solve({y^6-2*y^4+alpha=0,diff(y^6-2*y^4+alpha,y)=0},{alpha,y});

{y = 0, alpha = 0}, {y = 0, alpha = 0}, {y = 0, alpha = 0}, {alpha = 32/27, y = 2*RootOf(3*_Z^2-1)}

The origin appears three times, and the last one is really two different points.

> solve(3*Z^2-1=0);

1/3*sqrt(3), -1/3*sqrt(3)

This gives the two points as (32/27, 2*sqrt(3)/3 ) and (32/27, -2*sqrt(3)/3 ).  Numerically,

> evalf({[32/27,2*sqrt(3)/3],[32/27,-2*sqrt(3)/3]});

{[1.185185185, 1.154700539], [1.185185185, -1.154700539]}