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

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: • Oscar Limka suggestion is reasonable. The only equilibrium of this system is $(-\arcsin \frac{bv}{a} , v )$ and for $\vert v \vert > \frac{a}{b}$ it doesn't exist. It seems to be that Hopf bifurcation happens here, but this needs verifying (at least we should check eigenvalues for this equilibrium). Sorry, can't suggest more ideas now, need to sleep. Sep 11, 2015 at 22:39

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.

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)?

• You are correct that if $|v| > a/b$ then equilibrium is not possible. But I know that $v_{crit}$ is less than this, since the value is dependent on the initial values of $\theta$ and $\frac{d\theta}{dt}$. Sep 11, 2015 at 20:48
• Yes, I would split the study of this problem in two parts: (1) irrespective of initial values, discuss the existence of a stable equilibrium with respect to choices of parameters $a$, $b$ and $\nu$; (2) in those cases that a stable equilibrium exist, find its basin of attraction. Sep 11, 2015 at 20:55
• And definitely double variables by taking $\eta:=\dot\theta$. Sep 11, 2015 at 20:59
• Oskar, that sounds exactly what I should do, but I have no idea how to go about doing it. :-( Sep 11, 2015 at 21:08
• I'm no expert in dynamical systems, but looking at your equation again, I may have spoken too early, and it could be a bit more challenging. You have a genuinely non-autonomous system (because time enters the forcing term) in a nonlinear way. A text that comes to mind is Kloeden & Rasmussen (AMS, 2011). Sep 11, 2015 at 21:22

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)

• So, what is the "candidate Lyapunov equation" in this case? Sep 12, 2015 at 5:03