# Finding the period of a function of an interpolating function

I have two quantities of interest, x and y which are functions of $$\theta$$ and thus implicitly of time. They also depend on two parameters a and b. However, y is merely a dependent equation. They are defined and solved by

x[a_, b_] := a Cos[θ[t]]^2 + b Sin[θ[t]]^2 + θ'[t];
y[a_, b_] := 3 a^2 (Cos[2 θ[t]] + Sin[θ[t]]^2) - b Sin[θ[t]];
sol1 = First @ NDSolve[{x[1, 2] == 0, θ[0] == 0}, θ, {t, 0, 10}];


Clearly, y is a periodic function of time. However, if I try to find the period of this function using

FunctionPeriod[y[1, 2] /. sol1, t]


I get 0 every time. This seems to be a problem due to the nature of the solution as an interpolating function, but I was wondering if there is a good way to get the period.

• An exact solution is {\[Theta] -> Function[{t}, -ArcTan[Sqrt[b] Cos[Sqrt[a] Sqrt[b] t], Sqrt[a] Sin[Sqrt[a] Sqrt[b] t]]]} and then FunctionPeriod will work. – Michael E2 Sep 3 at 3:19
• Another way to get the exact period: The time for θ to go from 0 to 2π is equal to Integrate[1/(a Cos[\[Theta][t]]^2 + b Sin[\[Theta][t]]^2), {\[Theta][t], 0, 2 Pi}, Assumptions -> a > 0 && b > 0] – Michael E2 Sep 3 at 3:37
• An exact solution without the discontinuities of the ArcTan is {\[Theta] -> Function[{t}, -Sqrt[a] Sqrt[b] t + 2 ArcTan[ Cot[Sqrt[a] Sqrt[b] t] (-1 + Sqrt[Sec[Sqrt[a] Sqrt[b] t]^2])] - ArcTan[(Sqrt[a] Tan[Sqrt[a] Sqrt[b] t])/Sqrt[b]]]}  . – Akku14 Sep 3 at 5:10
• Another way to find the period of the numerical solution is define yy[t_] = y[1, 2] /. sol1[[1]]  and minimize the quadratic difference nmin = NMinimize[{(yy[3] - yy[3 + d])^2, 1 < d < 6}, d]  to get {4.84296*10^-14, {d -> 4.44288}}  with 4.4428 as period. Choose suitable starting point from a plot of y[1,2]. – Akku14 Sep 3 at 5:27
• So I realize that this equation has an exact solution - it was really chosen to illustrate a more general case where no exact solution exists, but that particular example would have been far too long to write out. Thank you for the suggestions though! – JamesVR Sep 4 at 0:22

Another approach

Plot the absolute value

absPlot = Plot[Abs@y[1, 2] /. sol1, {t, 0, 10}];


Extract points from the plot

points = absPlot // Cases[#, Line[x_] :> x, All] & // First;


Find peaks of the y values

yPeaks = FindPeaks[points[[All, 2]]];


Points corresponding to the peaks

peakPoints = points[[yPeaks[[All, 1]]]];


Plot to check

Plot[Abs@y[1, 2] /. sol1, {t, 0, 10}, Epilog -> {Red, PointSize[Large], Point[#] & /@ peakPoints}]


Depending on your definition of "period"

peakPoints // Select[Last@# > 2.5 &] // Part[#, All, 1] & // Differences
(* {1.56407, 2.87684, 1.56662, 2.87797} *)

peakPoints // Select[Last@# < 2.5 &] // Part[#, All, 1] & // Differences
(* {4.44382} *)


re Let's look at expressions you working with in more detail than you provide.

x[a_, b_] := a Cos[θ[t]]^2 + b Sin[θ[t]]^2 + θ'[t]
y[a_, b_] := 3 a^2 (Cos[2 θ[t]] + Sin[θ[t]]^2) - b Sin[θ[t]]


I think it's easier to work with θ as pure function, so I use NDSolveValue to solve for it.

θ = NDSolveValue[{x[1, 2] == 0, θ[0] == 0}, θ, {t, 0, 3 π}]


Now we can plot y[1, 2] over the domain of θ:

Plot[y[1, 2], {t, 0, 3 π}]


It does look periodic, doesn't it? But, of course, it's not. It's just a good polynomial approximation to a periodic function over a small closed interval of the reals. However, you can get a decent numerical approximation of the period of the function being approximated by computing the interval between the two negative minima shown int plot. Like so:

Subtract @@ Flatten[FindArgMin[y[1, 2], {t, #}] & /@ {8, 33/10}]

4.44288


Note: Should you want to create a real periodic function from y[1, 2], you will find methods for doing that by following this link