0

I have a system of two kinds of interacting oscillators. I represent these oscillator types with phi and theta. I want to extract, using Reap and WhenEvent, the times and values of the maxima of these oscillations. I know how to do this for phi OR for theta. I don't know how to do this for phi AND for theta.

I can extract maxima times and values, but I don't know how to keep them separate. I.e.--I can extract the set {t_phi_or_theta_max, phi[t_phi_max}, theta[t_theta_max]}, but I do not know how to extract {t_phi_max, phi[t_phi_max]} and {t_theta_max, theta[t_theta_max]}.

Phi and theta cannot be combined into a single indexed variable. Here, I present each as a single variable, but in the full code each are vectors of variables. The code below finds the peaks for phi, but I want to know how also to find peaks for theta.

tmax = 200;

osc1 = theta''[t] == -.468 theta[t] - .0044 theta'[t];
oscInit1 = {theta[0] == .66, theta'[0] == 0};

osc2 = phi''[t] == -Sin[phi[t]] - .0101 phi'[t] ((phi[t]/.33)^2 - 1) - .0007 Cos[phi[t]]  theta''[t];
oscInit2 = {phi[0] == .66, phi'[0] == 0};

vars1 = theta;
vars2 = phi;

(* defineEvent defines the peak finding conditional *)
defineEvents =
WhenEvent[phi'[t] == 0 && phi[t] > 0,
Sow[{t, phi[t]}]]; 

sol = Reap@NDSolve[{
 {osc1, osc2},
 {oscInit1, oscInit2},
 defineEvents
 },
{vars1, vars2},
{t, 0, tmax}];

I can plot the oscillations and the peaks of phi with the following code:

phiOUT = Evaluate[phi[t]] /. sol[[1]];
thetaOUT = Evaluate[theta[t] /. sol[[1]]];
Show[Plot[{phiOUT, thetaOUT}, {t, 0, 50}, 
PlotLegends -> {"phi", "theta"}],
ListPlot[sol[[2]], PlotMarkers -> Red]]

Giving the following plot:

enter image description here

KBL
  • 643
  • 3
  • 10

0 Answers0