# 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. – Evgeny Sep 11 '15 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}$. – Jason S Sep 11 '15 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. – Oskar Limka Sep 11 '15 at 20:55
• And definitely double variables by taking $\eta:=\dot\theta$. – Oskar Limka Sep 11 '15 at 20:59
• Oskar, that sounds exactly what I should do, but I have no idea how to go about doing it. :-( – Jason S Sep 11 '15 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). – Oskar Limka Sep 11 '15 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? – Evgeny Sep 12 '15 at 5:03