Phase response curve

Question

I just want to plot the phase response curve (phase vs. frequency curve) for the following nonlinear differential equation:

x''[t]*x[t]+3/2*x'[t]^2==x[t]^(-3*k)-1-4*\[Beta]*x'[t]/x[t]+f*Sin[\[CapitalOmega]*t]

Where k=7/5, f=0.2, \[Beta]=0.03 I have no idea that how can I get it numerically. Hoping for help. Thanks

Answer 1

The phase information follows from the optimal trigonometric approximation

(fundamentals of FourierSeries)

Clear["Global`*"]
k = 7/5; f = 1/5; β = 3/100;
tmax = 20;

sys = {
x''[t]*x[t] + 3/2*x'[t]^2 ==x[t]^(-3*k) - 1 - 4*β*x'[t]/x[t] + f*Sin[Ω*t],
x[0] == 1, x'[0] == 0} // Simplify;

X= ParametricNDSolveValue[sys, x, {t, 0, tmax}, {Ω}]

Thanks @BobHanlon for these code lines!

phase[\[CapitalOmega]_?NumericQ] := 
Block[{acs, 
M = NIntegrate[Outer[Times, {1, Cos[\[CapitalOmega] t],Sin[\[CapitalOmega] t]}, {1, Cos[\[CapitalOmega] t],Sin[\[CapitalOmega] t]}], {t, 0, tmax},Method -> "LocalAdaptive"],
rS = NIntegrate[X[\[CapitalOmega]][t] {1, Cos[\[CapitalOmega] t],Sin[\[CapitalOmega] t]}, {t, 0, tmax},Method -> "LocalAdaptive"]}, 
acs = LinearSolve[M, rS ];
ArcTan[acs[[2]], acs[[3]]]]

plot result

zw = Table[{\[CapitalOmega],phase[\[CapitalOmega]]}, {\[CapitalOmega], Subdivide[0 , 5, 25]}]

ListPlot[zw, Joined -> True,AxesLabel -> {\[CapitalOmega], "phase[\[CapitalOmega]]"}]

Answer 2

$Version

(* "14.0.0 for Mac OS X ARM (64-bit) (December 13, 2023)" *)

Clear["Global`*"]

k = 7/5; f = 1/5; β = 3/100;

tmax = 20;

sys = {
    x''[t]*x[t] + 3/2*x'[t]^2 == 
     x[t]^(-3*k) - 1 - 4*β*x'[t]/x[t] + f*Sin[Ω*t],
    x[0] == 1, x'[0] == 0} // Simplify;

sol = ParametricNDSolve[sys, x, {t, 0, tmax}, {Ω}]

Manipulate[
 Plot[x[Ω][t] /. sol, {t, 0, tmax},
  AxesLabel -> (Style[#, 14] & /@ {t, x})],
 {{Ω, 1}, -5, 5, 0.1, Appearance -> "Labeled"}]