# Determine differential equation is periodic and find its period

Assuming I have a second-order ordinary differential equation

$$[a- b\cos^2\theta] {\theta}'' + [b\cos\theta\sin\theta] {\theta '}^2 + g \sin\theta = 0$$

Where $$b>>a>0$$, and $$g<0$$ are constants. How can I determine if its solution is periodic and calculate its period?

I tried to calculate the expression, but it was too complicated and I couldn't get an analytical expression for the solution.

eqn = (a - b Cos[θ[t]]^2) θ''[t] +
b Cos[θ[t]] Sin[θ[t]] θ'[t]^2 +
g Sin[θ[t]] == 0;
sol = DSolve[{eqn, θ[0] == c1, θ'[0] == c2}, θ[t],
t];


The above code caused the application to freeze. Do I have any good way to substitute the condition of $$b>>a>0, g<0$$ into Dsolve in order to obtain an analytical solution?

I tried a set of initial values and found that this differential equation should be periodic.

g = -16;
b = 10;
a = 0;
ft = 10;
initialvalue = {\[Theta][0] == Pi/2 - Pi/1000, \[Theta]'[0] == 0};
eq = (a - b Cos[\[Theta][t]]^2) (\[Theta]^\[Prime]\[Prime])[t] +
b Cos[\[Theta][t]] Sin[\[Theta][t]] Derivative[1][\[Theta]][t]^2 +
g Sin[\[Theta][t]] == 0;
sol = NDSolve[{eq, initialvalue}, {\[Theta], \[Theta]'}, {t, 0, ft}];
theta = Plot[{\[Theta][t]*180/Pi} /. sol, {t, 0, ft},
PlotRange -> Automatic, PlotLabel -> "\[Theta](t)"];


My question is, do I have any good way to prove that this second-order differential equation is periodic? The above code can only reflect that the solution is periodic for a certain initial value, as seen from a numerical perspective.

Update

Perhaps we need to add some conditions to $$\theta$$ in order to prove it, such as $$\theta\in (-\pi/2,\pi/2)$$ or $$(-\pi/4,\pi/4)$$. Using Poincaré mapping for verification may be a good choice, but I am not sure how to substitute the assumption of coefficients $$b>>a>0$$, and $$g<0$$ into the function package of mma.

• Sorry, there are always problems with displaying symbols when using SE.editor in the second mma code.
– lumw
Commented Jun 16, 2023 at 16:03
• No chance. This equation is second order, quadratic first order, nonlinear with periodic coefficients and no first integarl in sight. Looks like the output of a free associative symbol producing automaton. Consider theta'' + sin theta =0, the pendulum equation. Has periodic ocillatory solutions, rotational solutions with positive oscillatory velocity and the 4 arctan soltution top to top in infinite time. Commented Jun 16, 2023 at 21:07
• @RolandF thank you. Perhaps it is necessary to add some restrictions to the initial values of the equation. I am considering it. From what I understand now, some tools like Poincaré mapping may be worth trying.
– lumw
Commented Jun 16, 2023 at 22:29
• why your in 1st block of code eqn does not look like your top math equation -- there is no even parameter a there ? Changing order of terms and breaking lines in random places does not help readability. Commented Jun 17, 2023 at 3:40
• @VitaliyKaurov There are some errors here, thank you for pointing them out. I have already updated them.
– lumw
Commented Jun 17, 2023 at 5:19

My current approach is to draw the vector field diagram through StreamPlot (First, make an intuitive judgment), and then use Poincare Bendixson theorem and Dulac's theorem/Criterion to judge. This is my code:

Firstly, convert it into a first-order equation. Let $$x(t) = \theta(t), y(t) = \theta(t)'$$
so we have $$x'(t)= y(t)$$ $$y(t)' = - (g \sin(x(t)) + b \cos(x(t)) \sin(x(t)) y(t)^2) / (a - b \cos(x(t))^2)$$

Then, Take some initial values ​​in these areas. The complete code is as follows.

g = -16;
b = 10;
a = 0;
ft = 10;
eqNew = (a - b Cos[\[Theta][t]]^2) (\[Theta]^\[Prime]\[Prime])[t] +
b Cos[\[Theta][t]] Sin[\[Theta][t]] Derivative[1][\[Theta]][t]^2 +
g Sin[\[Theta][t]] == 0;
streamPlot =
StreamPlot[{y, -(g Sin[x] + b Cos[x] Sin[x] y^2)/(a -
b Cos[x]^2)}, {x, -Pi/2, Pi/2}, {y, -Pi/2, Pi/2}]
colors = {Red, Green, Blue};
initialConditions = {{Pi/3, -Pi/2.4}, {Pi/2.4, 0}, {-Pi/3, Pi/2}};
solutions =
Table[NDSolve[{eqNew, \[Theta][0] ==
initialConditions[[i, 1]], \[Theta]'[0] ==
initialConditions[[i, 2]]}, {\[Theta], \[Theta]'}, {t, 0,
20}], {i, Length[initialConditions]}];
trajectoryPlots =
Table[ParametricPlot[
Evaluate[{\[Theta][t], \[Theta]'[t]} /. solutions[[i]]], {t, 0,
20}, PlotStyle -> colors[[i]]], {i, Length[initialConditions]}];
Show[streamPlot, trajectoryPlots]


Finally, attempting to apply Poincare Bendixson theorem and Dulac's theorem/criterion.

Update

Here is a more general analysis, by multiplying the original equation by $$2\theta'$$ and calculating the first integral: $$(a-b^2\cos^2θ)θ'^2+4g\sin^2(θ/2)=C$$

The level sets are closed curves around $$(\theta,\theta')=(0,0)$$. more details can see math.SE, the symbols used here are different in order to match the code of Mathematica.

This is just some thoughts to get you started on theoretical solution and to point out some technical tricks. I did not had time to go further yet. My version:

In[]:= $Version Out[]= 13.2.0 for Mac OS X x86 (64-bit) (November 18, 2022)  This version for this equation eqn = (a - b Cos[θ[t]]^2) θ''[t] + b Cos[θ[t]] Sin[θ[t]] θ'[t]^2 + g Sin[θ[t]] == 0;  gives the following sol=DSolveValue[eqn,θ,t,Assumptions->{b>a>0,g<0}] sol // TraditionalForm  So it does at least something with an inverse functions. Maybe this is enough to figure out periodicity or you might use perturbation theory for small $$a$$ etc. Note two details: • You do not need to put general non-numeric constants DSolve - it makes them automatically • Assumptions option is the thing you asked for and probably what makes it to compute Another thing, if you assume in zero approximation $$a=0$$ then you can derive 1st integral, which is a useful quantity in such systems. Then maybe you can use perturbation theory again for small $$a$$. Here is a sketch (could be some generalizations and errors; $$n=b$$): My memory of all this is quite rusty. But here is a (maybe wrong) thought. 1st integral is a conserved quantity $$c$$. If you plot phase-portrait momentum $$p=s'$$ vs. $$s$$ for various $$c$$ you can see the trajectories are bounded and periodic. ContourPlot[p^2-Sqrt[1-s^2],{s,-1,1},{p,-1,1}]  This might translate to $$θ$$ and perhaps small enough perturbation term $$~a$$ will not break periodicity. Another thought - potentially (not sure) you might find a general transformation of coordinates that will allow to keep term $$a$$ and find exact 1st integral. Also you might be able to find Hamiltonian and try to see if trajectories of phase portrait there are periodic. • Thank you, this is a new idea compared to Bendixson's theorem! Great! And I didn't expect DSolveValue to support Assumptions, it's really great. – lumw Commented Jun 17, 2023 at 9:59 • The calculation is correct. By the way, how can we explain the periodicity of the solution for$\frac{(s')^2}{2}+C= k\sqrt{1-s^2}$? as$s=\sin\theta$, The expression$s[t]$is also very complicated. – lumw Commented Jun 17, 2023 at 14:39 • @lumw - I sketched some thoughts at the end. Commented Jun 17, 2023 at 19:45 • You're right, the first integral does tell the story. The first integral can also be calculated by keeping$a$, which can be achieved by multiplying$2\theta(t)'\$ in the original second order equation. But your answer provides a more general way of thinking, that is, when you can't get the first integral. Thank you!
– lumw
Commented Jun 17, 2023 at 23:40