4
$\begingroup$

As I was following a similar question I came across a beautiful answer here. However, in that code, he used InternalBag[] which heavily depends on ram that I can't afford and I need to make a data file too. So Please Help me edit this code so that I can export the data to my hard disk.

Along with that, this code only stores the extremum values of Z as in WhenEvent there is only one condition i.e. z'[t]==0. I want to store minimum and maximum values in two different columns. How do I do that?

res = Internal`Bag[];(*a place to store results*)
tmax = 100;(*how long to run for each r value*)
{x0, y0, z0} = {10, 10, 10};(*initial ICs*)

Do[sol = With[{a = 10.0, c = 8/3},
    NDSolve[{x'[t] == a (y[t] - x[t]), 
      y'[t] == r x[t] - y[t] - x[t] z[t], z'[t] == x[t] y[t] - c*z[t],
       x[0] == x0, y[0] == y0, z[0] == z0,(*save extrema of z[t]*)
      WhenEvent[z'[t] == 0, Internal`StuffBag[res, {r, z[t]}]]}, {x, 
  y, z}, {t, 0, tmax}, Method -> "StiffnessSwitching",MaxSteps -> \[Infinity]]][[1]];
 (*save end value for next ICs*){x0, y0, z0} = {x[tmax], y[tmax], z[tmax]} /. sol;, {r, 100, 10, -0.1}];

ListPlot[Internal`BagPart[res, All], PlotStyle -> {Red, Opacity[0.1], PointSize[0.001]}, AxesLabel -> {"r", "z"}]

I was suggested to use Reap or Sow, but I am unable to edit the code properly

$\endgroup$
2
  • $\begingroup$ Glad you got an answer to your question. BTW, where did you see that Internal`Bag is obsolete? $\endgroup$
    – Chris K
    Aug 23, 2022 at 12:26
  • $\begingroup$ @ChrisK Sorry, I was wrong there. It's not obsolete but it heavily depends on ram. $\endgroup$
    – user444
    Aug 24, 2022 at 3:40

1 Answer 1

7
$\begingroup$

Here is an example how you would use Sow and Reap. Note, I decreased the r interval, as I did not want to wait a long time.

tmax = 100;(*how long to run for each r value*)
{x0, y0, z0} = {10, 10, 10};(*initial ICs*)

res = Reap[
    Do[sol = 
      With[{a = 10.0, c = 8/3}, 
        NDSolve[{x'[t] == a (y[t] - x[t]), 
          y'[t] == r x[t] - y[t] - x[t] z[t], 
          z'[t] == x[t] y[t] - c*z[t], x[0] == x0, y[0] == y0, 
          z[0] == z0,(*save extrema of z[t]*)
          WhenEvent[z'[t] == 0, Sow[{r, z[t]}]]}, {x, y, z}, {t, 0, 
          tmax}, Method -> "StiffnessSwitching", 
         MaxSteps -> \[Infinity]]][[1]];
     (*save end value for next ICs*){x0, y0, 
       z0} = {x[tmax], y[tmax], z[tmax]} /. sol;, {r, 11, 
      10, -0.1}]][[2]];

ListPlot[res, PlotStyle -> {Red, Opacity[0.1], PointSize[0.001]}, 
 AxesLabel -> {"r", "z"}]
$\endgroup$
2
  • $\begingroup$ Thank you. How do I separate minimum and maximum values in WhenEvent? $\endgroup$
    – user444
    Aug 23, 2022 at 9:15
  • $\begingroup$ In the WhenEvent you only have access to 3 data: t, z[t] and z'[t]. That is not enough to determine if we have a min or max. But "sol" contains the whole function and you may use it to determine if it is min or max. $\endgroup$ Aug 23, 2022 at 10:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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