3

Possible Duplicate:
How to find all the local minima/maxima in a range

I have the solution of the following non-linear system:

sol1 = NDSolve[
  {x'[t] == -(y[t] + z[t]),
   y'[t] == x[t] + 0.2 y[t],
   z'[t] == 0.2 + x[t] z[t] - 5.7 z[t],
   x[0] == 1, y[0] == 1, z[0] == 1
  },
  {x, y, z},
  {t, 0, 100}
  ]

How can I find the $k^{th}$ local maximum of $z(t)$, i.e. $z(k)$, and then plot $z(k+1)$ vs. $z(k)$? There is an example in the "Mapping local maxima" section in Rössler attractor's wiki page. I am working with Wolfram Mathematica 8.0.

Community
  • 1
user5267
  • 321
  • 1
  • 3
  • 5

1 Answers1

4

The maxima will occur at points where the derivative is zero and, except in special cases, they will alternate with minima. You can easily detect where the derivative zero using event detection. In V9, you do this like so.

Clear[x, y, z, sol, pts];
{{sol}, {pts}} = Reap[
   NDSolve[
    {x'[t] == -(y[t] + z[t]),
     y'[t] == x[t] + 0.2 y[t],
     z'[t] == 0.2 + x[t] z[t] - 5.7 z[t],
     x[0] == 1, y[0] == 1, z[0] == 1,
     WhenEvent[z'[t] == 0, Sow[t]]
     },
    {x, y, z}, {t, 0, 100}
    ]];
z = z /. sol;
maxPts = Last /@ Partition[pts, 2];
Plot[z[t], {t, 0, 100}, PlotRange -> All,
 Epilog -> Point[{#, z[#]} & /@ maxPts]]

enter image description here

Note that the extremes are found reliably during the solution of the differential equation and there's no need to numerically solve equations involving interpolating functions afterward.

Now that you've got the maxima in a list, you can do anything you want with them.

Mark McClure
  • 32,469
  • 3
  • 103
  • 161
  • I like the WhenEvent approach, no need to worry about oddities of numerical optimization algorithms :) – ssch Jan 04 '13 at 01:49
  • Thank you for your answer. Do you know how can I do this in the Wolfram Mathematica 8.0? And how can I name the k-th local maximum of the z(t), z(k) and then plot it versus z(k+1)? There is an example in the "Mapping local maxima" of the following link: http://en.wikipedia.org/wiki/R%C3%B6ssler_attractor Thank you again for your help. – user5267 Jan 04 '13 at 16:46
  • @user5267 look at the question Jens linked, the same idea is done in a way that works in v8 there – ssch Jan 04 '13 at 22:56
  • @user5267 and regarding the plotting, the result you'll get will be a list of the values ${z_k}_{k=1}^n$ (after taking away the local minima in the way Mark did, or some other way). Example to plot in that manner pts = Range[10]; ListPlot[{Most@pts, Rest@pts}\[Transpose]] – ssch Jan 04 '13 at 23:00