2
$\begingroup$

I have been trying to find the peaks of a function I have plotted using ParametricNDSolve. I have to find these peaks to calculate the amplitude of all the various waves in the observed output. By amplitude here I mean the difference between an adjacent crest and trough. I shall attach the program here to clarify the problem.

Clear["Global`*"]
om = 1;
k = 1;
L = 1/1000;
P = 1.3;

eqns = {a'[t] == 
I*((de/om)*a[t] - (b[t] + Conjugate[b[t]])*a[t] - 
    1/2) - (k/(2*om))*a[t], 
b'[t] == -I*((P*(Abs[a[t]])^2)/2 + b[t]) - (L/(2*om))*b[t], 
b[0] == 0, a[0] == 0};

s = ParametricNDSolve[eqns, {a, b}, {t, 0, 100}, de, 
MaxSteps -> \[Infinity]];

x[de_][t_] = (b[de][t] + Conjugate[b[de][t]]) /. s;

Manipulate[
Plot[x[de][t], {t, 0, 99}, PlotRange -> {-4, 4}], {de, -3, 3, .1, 
Appearance -> "Labeled"}]

Now I have to find the value of the amplitude of x for every value of de that I have plotted. While I have been able to find the peaks for individual values of de and a single peak in for these values, I have been unable to find all the peaks and amplitudes for all the values of de. I used various function including FindPeaks[] to try to find such values but was unable to do so. Normally I would have used for loops but that seems to advised against in Mathematica. I would also like to know whether such a code would be unwieldy and whether I use brute force to find such values.

I would be very grateful to anyone who could help me solve this dilemma. I would like to thank you for it beforehand

$\endgroup$

2 Answers 2

2
$\begingroup$

It might be easiest to detect the extrema within the NDSolve using WhenEvent.

Clear["Global`*"]
om = 1; k = 1; L = 1/1000; P = 1.3;
tmax = 10000;

eqns = {
  a'[t] == I*((de/om)*a[t] - (b[t] + Conjugate[b[t]])*a[t] - 1/2) - (k/(2*om))*a[t], 
  b'[t] == -I*((P*(Abs[a[t]])^2)/2 + b[t]) - (L/(2*om))*b[t], 
  b[0] == b0, a[0] == a0,
  WhenEvent[b'[t] + Conjugate[b'[t]] == 0, Sow[{t, b[t] + Conjugate[b[t]]}]]
};

de = -3;
{a0, b0} = {0, 0};
{s, ex} = Reap[NDSolve[eqns, {a, b}, {t, 0, tmax}, MaxSteps -> \[Infinity]][[1]]];
(* plot saved extrema *)
ListPlot[ex[[1]], AxesLabel -> {"t", "extrema"}]

Mathematica graphics

If you're interested in the long-term behavior, then it looks like tmax>100 is probably a good idea.

Then you can loop over de, saving the difference between the last two extrema as the amplitude (assuming a simple cycle) and using the final values of a and b as initial conditions for the next value of de.

Clear[de];
{a0, b0} = {a[tmax], b[tmax]} /. s;
tmax = 1000;
res = {};
Do[
  {s, ex} = Reap[NDSolve[eqns, {a, b}, {t, 0, tmax}, MaxSteps -> \[Infinity]][[1]]];
  {a0, b0} = {a[tmax], b[tmax]} /. s; (* save final values *)
  AppendTo[res, {de, Abs[ex[[1, -1, 2]] - ex[[1, -2, 2]]]}];
, {de, -3, 3, 0.1}];
(* plot amplitude vs de *)
ListPlot[res, AxesLabel -> {"de", "amp"}]

Mathematica graphics

That jump at de = 0.7 makes me think there might be alternative attractors. Indeed, if we generate the bifurcation diagram starting at de = 3 and moving towards lower de, we get a different result.

de = 3;
tmax = 10000;
{a0, b0} = {0, 0};
{s, ex} =  Reap[NDSolve[eqns, {a, b}, {t, 0, tmax}, MaxSteps -> \[Infinity]][[1]]];

Clear[de];
{a0, b0} = {a[tmax], b[tmax]} /. s;
tmax = 1000;
res = {};
Do[
  {s, ex} = Reap[NDSolve[eqns, {a, b}, {t, 0, tmax}, MaxSteps -> \[Infinity]][[1]]];
  {a0, b0} = {a[tmax], b[tmax]} /. s;
  AppendTo[res, {de, Abs[ex[[1, -1, 2]] - ex[[1, -2, 2]]]}];
, {de, 3, -3, -0.1}];
ListPlot[res, AxesLabel -> {"de", "amp"}, PlotStyle -> Orange]

Mathematica graphics

So a full understanding of the model's dynamics will require a bit of care!

$\endgroup$
2
$\begingroup$

Not a complete answer, hopefully helpful:

Trying to visualize your problem I get

pic=Plot3D[x[de][t], {t, 0, 99}, {de, -3, 3}, PlotRange -> {-4, 4},MaxRecursion -> 4]

enter image description here

a surface with many peaks.

If you split the region t,de NMaximize finds the global(!) peak in the subregion

NMaximize[{x[de][t], 0 < t < 10, -3 < de < 3}, {t, de}]
(*{0.042377, {de -> 2.20983, t -> 5.55602}}*)
$\endgroup$

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.