Technique for predicting attractor capture in nonlinear differential equations? (quasi-pendulum equation)

Question

I'm working on simulating this equation (application is motor control, not that it matters):

$$\frac{d^2\theta}{dt^2}+b\frac{d\theta}{dt}=a \sin (x-\theta)$$

where $x = vt$ for $t > 0$, and I'm finding that for given initial conditions $\frac{d\theta}{dt}|_{t=0}$ and $\theta|_{t=0}$, there seems to be a critical value $v_{crit}$ such that:

• if $v < v_{crit}$, $x-\theta$ oscillates but settles down to its equilibrium value $\phi = \sin^{-1} \frac{bv}{a}$ ("capture")
• if $v > v_{crit}$, $x-\theta$ tends to a linearly increasing difference + a small oscillating term.

Predicting the value of $v_{crit}$ is important in my application, and I would like to understand what is going on.

My training in nonlinear differential equations is rather limited and rusty, and I never took advanced classes... I think there might be some insight using energy techniques (Hamiltonians? Lyapunov stability?) since there is similarity to a driven pendulum equation ($\frac{d^2\theta}{dt^2}+\frac{g}{l}\sin \theta = u(t)$), but I can't figure out what.

How can I figure out this critical value?

Can anyone point me towards some reference material (or even the right terms to look up) so I could learn a technique to solve my problem?

It seems like if I can show that $|x-\theta|$ reverses direction before it hits $\pi$, then capture is guaranteed.

---

okay, splitting into 1st-order systems:

$\begin{eqnarray} \dot{\omega} &=& -a \sin(\theta - x) - b\omega \cr \dot{\theta} &=& \omega \end{eqnarray}$

Change of variable $u = \theta - x$ so $\dot{u} = \omega - v$ and $\ddot{u} = \ddot{\theta} = \dot{\omega}$:

$\begin{eqnarray} \dot{\omega} &=& -a \sin u - b\omega \cr \dot{u} &=& \omega - v \end{eqnarray}$

---

If I try to write a Lyapunov equation $E = c\omega^2 + d \cos u$ I get

$$\begin{eqnarray} \dot{E} &=& 2c\dot{\omega}\omega - d \dot{u} \sin u \cr &=&2c\omega(-a\sin u - b\omega) - d (\omega - v) \sin u \cr &=&\omega \sin u (-2ac - d) -2bc\omega^2 + dv \sin u \end{eqnarray} $$

I can make the first term go away if I choose d=-2ac; the second term is negative if $c>0$ but I can't get rid of the third term.

---

Attempt #2: $E = (\omega-v)^2 - 2a \cos u$ I get

$$\begin{eqnarray} \dot{E} &=& 2\dot{\omega}(\omega-v) +2a \dot{u} \sin u \cr &=&2(\omega-v)(-a\sin u - b\omega) +2a (\omega - v) \sin u \cr &=&(\omega-v) \sin u (-2a + 2a) -2b\omega(\omega-v) \cr &=& -2b\omega(\omega-v) \cr &=& -2b\left(\left(\omega-\frac{v}{2}\right)^2 - \frac{v^2}{4}\right) \cr \end{eqnarray} $$

but that's not necessarily negative. Urk.

---

Some numerical sample points: I'm using a = 2.3086177e5, b = 1.78179103 (max v with equilibrium at $a/b \approx 129567$, but $v_{crit}$ tends to be much smaller in practice), and in my simulations I'm seeing:

$\omega|_{t=0} = 0, \theta|_{t=0} = 0 : v_{crit} \approx 958.929$

$\omega|_{t=0} = 0, \theta|_{t=0} = 0.5 : v_{crit} \approx 929.408$

$\omega|_{t=0} = 0, \theta|_{t=0} = 1.0 : v_{crit} \approx 842.336$

$\omega|_{t=0} = 0, \theta|_{t=0} = \pi/2 : v_{crit} \approx 680.156$

For the last case, here's a phase plot (it's a very underdamped system so the turns of the plotted curve come very close together; the cusp on the left is where things slow down for a moment):

and a timeseries plot:

Answer 1

For the variable change $w=vt -\theta$ we have $\dot{w}=v-\dot{\theta}$ and $\ddot{w}=-\ddot{\theta}$. Thus your initial ODE becomes in the new set of variables $(w,\dot{w})$ $$\ddot{w}+b\dot{w}+a\sin w=bv$$ This is a pendulum equation with linear dissipation and constant torque and it appears in the analysis of charged-density-waves (CDW) (see L.-G. Li, Y.-F. Ruan, The analysis on the single particle model of CDW, Physics Letters A, vol. 372, issue 42, pp. 6443–6447,2008). In this paper, they examine the ODE $$\ddot{\phi}+\Gamma \dot{\phi}+\sin\phi=\beta$$ and prove the existence of a critical value $\beta_0$ such that a stable periodic solution occurs whenever $\beta\geq \beta_0$.

No specific value for $\beta_0$ is obtained as their proof is based on the qualitative properties of the ODEs. These results also appear in Lian-Gang Li, arXiv:0807.3288v2 . Problems of this type have been initially considered in (M. Urabe, J. Sci. Hiroshima Univ. A, 18 (1954), p. 379) but I could not find a copy of this paper.

Answer 2

Assuming what you're saying is correct, a wild guess is $v_{\text{crit}}=a/b$ (just because anything above that would have no arcsin, and your first case cannot work)?

Answer 3

Hmm. Empirically I'm noticing something:

Some numerical sample points: I'm using a = 2.3086177e5, b = 1.78179103 (max v with equilibrium at $a/b \approx 129567$, but $v_{crit}$ tends to be much smaller in practice), and in my simulations I'm seeing:

$\omega|_{t=0} = 0, \theta|_{t=0} = 0 : v_{crit} \approx 958.929$

$\omega|_{t=0} = 0, \theta|_{t=0} = 0.5 : v_{crit} \approx 929.408$

$\omega|_{t=0} = 0, \theta|_{t=0} = 1.0 : v_{crit} \approx 842.336$

$\omega|_{t=0} = 0, \theta|_{t=0} = \pi/2 : v_{crit} \approx 680.156$

These numbers are close to $\sqrt{2a(1+\cos \theta|_{t=0})}$

(derived from $v_{crit}^2 = 2a(1+\cos \theta|_{t=0})$, similar to my candidate Lyapunov equation)