How to calculate and plot instantaneous phase in Mathematica

I would like a plot of the instantaneous phase difference between a frequency-swept drive and the nonlinear oscillator it is driving. x[t] below is the instantaneous displacement of the oscillator and plotx provides a plot.

Thanks, Carey

s =
 NDSolve[{x''[t] + x[t] - 0.167 x[t]^3 == 
    0.005 Cos[t - 0.5*0.0000652*t^2], x[0] == 0, x'[0] == 0}, 
  x, {t, 0, 3000}, MaxSteps -> 35000]

plotx = Plot[Evaluate[x[t] /. s], {t, 0, 3000}, PlotPoints -> 10000, 
  Frame -> {True, True, False, False}, FrameLabel -> {"t", "开发者_如何学运维x"}, 
  FrameStyle -> Directive[FontSize -> 15], PlotLabel -> "(a)", 
  Axes -> False]

---

(Response, take 2)

You can get a reasonable approximation of the phase with

f[tt_?NumericQ] := -(ArcTan @@ ({x[t], x'[t]}/
    Sqrt[x[t]^2 + x'[t]^2]) /. s[[1]]) /. t -> tt

Here are some plots. First we show the driving term and the result together. It indicates they are a bit out of phase.

plotx2 = Plot[
  Evaluate[{x[t], Cos[t - 0.5*0.0000652*t^2]/5} /. s], {t, 0, 100}, 
  Frame -> {True, True, False, False}, FrameLabel -> {"t", "x"}]

Now we show the two phases together. I plot over a slightly different range this time.

phaseangles = 
 Plot[{f[t], Mod[t - 0.5*0.0000652*t^2, 2*Pi, -Pi]}, {t, 100, 120}, 
  Frame -> {True, True, False, False}, FrameLabel -> {"t", "x"}]

Last we show the phase differences.

phasediffs = 
 Plot[{f[t] - Mod[t - 0.5*0.0000652*t^2, 2*Pi, -Pi]}, {t, 100, 120}, 
  Frame -> {True, True, False, False}, FrameLabel -> {"t", "x"}]

Possibly I'm off by something additive (those Mod[] terms get bothersome), but this should give an idea of how one might proceed.

Daniel Lichtblau Wolfram Research

---

I'd look very closely at the method of averaging. In Strogatz's implementation, both the average envelope and phase of a nonlinear oscillator are found. Since you are looking for something a little bit beyond the first order, I'd consider looking at this paper from the Air Force Academy.