# Phase response curve

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

• What is [CapitalOmega]?
– Moo
Commented Mar 19 at 14:37
• \[CapitalOmega] is the forcing frequency. Commented Mar 19 at 15:06
• Right, but don't we need a value for a numerical solution along with two initial conditions?
– Moo
Commented Mar 19 at 15:12
• Sorry it's my mistake. The initial conditions are x[0]=1 and x'[0]=0. I want to plot it for the steady state condition and \[CapitalOmega] is varying from -5 to 5. Commented Mar 19 at 15:24

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]]"}]


• Thanks @Ulrich Neumann. It worked for my problem. Commented Mar 20 at 6:42
• You're welcome! Commented Mar 20 at 9:51
\$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"}]


• Thanks @Bob Hanlon, but I want to plot steady state phase vs. frequency curve Commented Mar 19 at 15:58