How to attain the period of this nonlinear differential equation/system?

Question

Lately, I've been trying to find the period of an angle included in the following differential equations, but only could with the basic model:

Basic or original: $$\mathrm{For}\ (\Phi (0), \Omega (0))=(\Phi_{o},0),\ \frac{d^2\Phi}{dt^2}= \frac{g}{\ell_{o}}\sin{\Phi}-\frac{g}{\ell_{o}}\zeta\ \mathrm{sgn\ \Phi}\ ;$$
Modified: $$\mathrm{For\ the\ same\ initial\ conditions},\ \frac{d^2\Phi}{dt^2}= \frac{g}{\ell_{o}}\frac{\sin{\Phi}}{f(\Phi)}-\frac{g}{\ell_{o}}\zeta \frac{\mathrm{sgn\ \Phi}}{f(\Phi)}\ -2\dot{\Phi}^2 \frac{f'(\Phi)}{f(\Phi)}.$$

Where $g$ is gravity, $\ell_{o}$ is the length of the inverted pendulum, $\zeta$ a group of other constants, $\operatorname{sgn}\left(\cdot\right)$ is the signum function, $\dot{\Phi}=\frac{d\Phi}{dt}$, $f(\Phi)=\sqrt[3]{1-\eta\cos{\Phi}}$ ($\eta$ is another constant) and $f'(\Phi)=\frac{df(\Phi)}{d\Phi}$.

And so, the method I used to get the period was basically this:

Let $F(\Phi)= \frac{g}{\ell_{o}}\sin{\Phi}-\frac{g}{\ell_{o}}\zeta\ \mathrm{sgn\ \Phi}$ , then the diff. eq. reduces to $\frac{d^2\Phi}{dt^2}=F(\Phi).$ And now I just proceed. \begin{align} \int \frac{d^2\Phi}{dt^2}d\Phi &= \int F(\Phi)\ d\Phi\\ \frac{1}{2}\dot{\Phi}^2 &= \int F(\Phi)\ d\Phi\ +C\\ \dot{\Phi} &= \frac{d\Phi}{dt} = \sqrt{2\int F(\Phi)\ d\Phi +C}\\ \frac{T}{4}=\int_{t_{o}}^{t_{1}}dt &= \int_{0}^{\Phi_{o}}\frac{d\Phi}{\sqrt{2\int F(\Phi)\ d\Phi +C}}\\ T &=2\sqrt{2} \int_{0}^{\Phi_{o}}\frac{d\Phi}{\sqrt{\int F(\Phi)\ d\Phi +C}}. \end{align}

This worked for the basic model; but didn't for the modified one. The issue was the integral of $F(\Phi)$ since in the modified version it included all terms divided by $f(\Phi)$ and also the $\dot{\Phi}^2 \frac{f'(\Phi)}{f(\Phi)}$ one too. Can someone tell me any easier way to attain the period of this modified system? Or what approximation could I use to make it easier to deal with?

Answer 1

This likely isn't going to be the answer you want, as it obtains a numeric solution instead of analytic one. But this is still an approach I would take in finding it or searching for a parameter space. The general idea would be to solve this as a pair of first-order ODEs: \begin{align} \dot{\phi} &= \chi \\ \dot{\chi} &= F\left(t,\,\phi,\,\chi\right) \end{align} From which you can use an ODE solving method you like (e.g., RK4). In Python, you could use scipy.integrate.odeint to solve the ODE for you (assuming you have a good time span at the onset). I left a couple parameters ($\eta$, $\zeta$ and $\ell_0$) as ? as their values aren't really indicated in the post.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# define constants..
ETA = ?
ZETA = ?
G = 9.81
L0 = ?
GL0 = G / L0

def df(phi: float) -> float:
  """ derivative of f w.r.t. phi """
  num = np.sin(phi)
  if abs(num) < 1e-10:
    return 0
  return ETA * num / (3.0 * (1.0 - ETA * np.cos(phi))**(2.0 / 3.0))


def f(phi: float) -> float:
  """ function f() """
  return (1.0 - ETA * np.cos(phi))**(1.0 / 3.0)


def evolve(z: np.ndarray, t: float) -> np.ndarray:
  """ evolution function """
  chi, phi = z
  dchi = GL0 * np.sin(phi) - GL0 * ZETA * np.sign(phi) - 2.0 * df(phi)
  dphi = chi
  return np.array([dchi / f(phi), dphi])


def ode(phi_init: np.ndarray):
  """ solve the ODE problem given initial conditions """
  t_span = np.arange(0.0, 6.0 * np.pi, 0.05)
  solution, output = odeint(evolve, phi_init, t_span, full_output = True)
  # plot solution
  phi, phidot = solution.T
  plt.plot(t_span, phi, label=r"$\Phi$")
  plt.plot(t_span, phidot, label=r"$\Omega$")
  plt.legend()
  plt.show()

And then calling the ode function with whatever initial conditions you want. Using $\ell_0=12$, $\eta=0.15$ and $\zeta=0.25$ with $\Phi_0=\pi/4,\,\Omega=0$, I get the following,

Presumably you can replicate the image in the original post using whatever parameters the source indicates. Toying around with some of the parameters, it looks like $\zeta$ is a likely larger than 1.