3
$\begingroup$

I'm trying to plot a bifurcation map for a Lorenz System. I have a function (called estables in the code that follows) which, for a certain value r, returns the list of stationary points for the system with that parameter r.

Then I do this:

r0 = 0;
rf = 10;
dr = 0.5;
ys = Table[estables[ra], {ra, r0, rf, dr}]

ListPlot[ Flatten[ys, {2}], PlotRange -> All , DataRange -> {r0, rf}, 
PlotStyle -> {Blue, PointSize[Small]}, Joined -> False ]

And I get a weird thing where the points are plotted out of order. That is, the points in the first line are not the ones plotted above the 0, as it should be.

What would be the best way to achieve this?

Edit: This are the differential equations:

ecf = σ (p[t] - f[t]);
ecp = -p[t] + d[t] f[t];
ecd = b (r - d[t] - f[t] p[t]);

and this is the function

t0 = 100;
tf = 105;
dt = 0.05;

estables[ra_] :=

Module[{par, solnum, fs}, 
    par = {σ -> 3., b -> 1, r -> ra}; 
    solnum = 
      NDSolve[{Derivative[1][f][t] == ecf, 
               Derivative[1][p][t] == ecp, 
               Derivative[1][d][t] == ecd, 
               f[0] == 0.00001, p[0] == 0., 
               d[0] == 1} /. par, 
               {f, p, d}, {t, t0, tf}, 
               MaxSteps -> 100000
    ];
    fs = Flatten[Table[Evaluate[f[t] /. solnum], {t, t0, tf, dt}]]; 
    Union[FindPeaks[fs], -FindPeaks[-fs]]][[;; , 2]]
$\endgroup$
  • $\begingroup$ Including the definition of estables would certainly help us to help you. $\endgroup$ – J. M. will be back soon May 21 '16 at 13:08
  • $\begingroup$ I didn't include it for brevity, but yes, you're probably right. $\endgroup$ – carllacan May 21 '16 at 13:12
  • $\begingroup$ Where are t0, tf, and dt? $\endgroup$ – J. M. will be back soon May 21 '16 at 13:19
  • 1
    $\begingroup$ I have a suggestion: why not use WhenEvent[] instead to capture the points where f'[t] == 0? $\endgroup$ – J. M. will be back soon May 21 '16 at 13:50
  • $\begingroup$ Mmm, I suspected there would be a more sophisticated way to capture the stable points, that looks good. I'm not sure how to use it, though. Which action should be triggered when f'[t] ==0? $\endgroup$ – carllacan May 21 '16 at 13:56
4
$\begingroup$

Thanks to J.M. I found a more elegant way to obtain the statinoary values that returned far simpler data, which was easy to plot. This is how my code looks now:

ecf = s (p[t] - f[t]);
ecp = -p[t] + d[t] f[t];
ecd = b (r - d[t] - f[t] p[t]);

par = {s -> 3., b -> 1, r -> 100};
t0 = 800;
tf = 850;



estables[ra_] :=

 Module[{par, solnum,  points}, 
  par = {s -> 3., b -> 1, r -> ra}; 
  {solnum, points} = Reap@NDSolve[{
           Derivative[1][f][t] == ecf, 
           Derivative[1][p][t] == ecp, 
           Derivative[1][d][t] == ecd, 
           f[0] == 0.00001, p[0] == 0., 
           d[0] == 1, 
           WhenEvent[{f'[t] == 0 && t > t0}, Sow[f[t]]]} 
           /. par, {f, p, d}, {t, t0, tf}, MaxSteps -> 100000];
  Flatten[points]
];


r0 = 0;
rf = 105;
dr = 0.2;

ys = Table[{ra, e}, {ra, r0, rf, dr}, {e, estables[ra]}];
ListPlot[ Flatten[ys, 1], PlotRange -> All , DataRange -> {r0, rf}, 
 PlotStyle -> {Blue, PointSize[0.002]}, Joined -> False, 
 AxesLabel -> {r, f} ]
$\endgroup$
  • $\begingroup$ Great job on self answer! When I copy and run your code it appears that points 2 and 3 generate an error. The rest of the points appear to be OK. $\endgroup$ – Jack LaVigne May 22 '16 at 13:50
  • $\begingroup$ Actually, I've just found out that my code is wrong. It gives a close enough diagram, but a lot of "spurious" points appear that shouldn't be there. Is that what you mean? I've almost managed to fix it, I will edit in a moment. $\endgroup$ – carllacan May 22 '16 at 20:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.