3
$\begingroup$

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.

Improve this question
$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$
    – Michael E2
    Commented Sep 3, 2020 at 3:19
  • $\begingroup$ 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] $\endgroup$
    – Michael E2
    Commented Sep 3, 2020 at 3:37
  • $\begingroup$ 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]]]} . $\endgroup$
    – Akku14
    Commented Sep 3, 2020 at 5:10
  • $\begingroup$ 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]. $\endgroup$
    – Akku14
    Commented Sep 3, 2020 at 5:27
  • $\begingroup$ 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! $\endgroup$
    – JamesVR
    Commented Sep 4, 2020 at 0:22

2 Answers 2

Reset to default
2
$\begingroup$

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

enter image description here

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} *)
Improve this answer
$\endgroup$
1
$\begingroup$

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

plot

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

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.