Analyzing Zeeman's heartbeat equations. Proving that there is one and only limit cycle

Question

From Poincare Bendixon Theorem we know that a non empty compact limit set of a continuously differentiable dynamical system in the plane which contains no equilibrium point is a closed orbit.

I'm trying to analyse Zeeman's heartbeat equations, \begin{align}ϵ\,\frac{dx}{dt}&=−(x^3−Tx+b),\;T>0\\ \frac{db}{dt}&=(x−x_0),\end{align} and I have to prove that there exist a closed orbit (limit cycle) and that it's the only one. So I need to use the above mentioned theorem to prove that there exist a limit cycle.

I've already found out that the origin of my system is totally unstable (which is (0,0) for $\epsilon = 0.2$, $T=1$ and $x_0=0$). Also all the trajectories flow away from the origin. Generally the critical point of my system is $$(x,b)=(x_0,Tx_0-x_0^3).$$ Do I consider this point to be a unique equilibrium point? Theorem says that there must be no equilibrium point, so what if there exists totally unstable critical point? Is there any other method for proving that there exist one and only limit cycle?

Answer 1

Yes, by the structure of the equation it is the only stationary point.

If $b+x^3-Tx>\sqrtϵ$, then $\dot x\sim -\frac1{\sqrtϵ}$ moves the flow fast to the left stopping at the falling branch of $b=Tx-x^3$. Similarly for $b+x^3-Tx<-\sqrtϵ$. For a point close (meaning with distance $O(ϵ)$) to the falling segments of the graph of the cubic, and assuming that $x_0$ is between the locations of the maxima and minima of the cubic, the first equation forces the flow towards the cubic while the second one moves $b$ towards the extremum next to the current position, downwards for negative $x$ and upwards for positive $x$. Thus the flow of the system creeps down the negative falling branch until it reaches the minimum where it leaves the cubic and is blown to the positive falling branch. There it creeps upwards until it reaches the maximum, is blown to the left and so on.

To find a contour where the vector field points uniformly inwards seems complicated.

---

The following is incomplete. While the theory of Lienard systems is more general, the described Lyapunov functions work only under the assumption $c(u)\ge 0$.

---

Consider the second order form $ϵ\ddot x+(3x^2-T)\dot x+(x-x_0)=0$ and compare to the form of a Lienard equation $$ u''+c(u)u'+g(u) $$ as in Stability in a Liénard system. This is achieved by setting $u=x-x_0$, $g(u)=\frac uϵ$, $c(u)=\frac1ϵ(3(x_0+u)^2-T)$. Then the Lyapunov function proposed in the answer is, using $p=\dot u$ $$ V(u,p)=\frac12\left(p+C(u)\right)^2+\frac12p^2+2G(u),\\~\\C'=c,~G'=g,~C(0)=0=G(0) $$ with derivative along the solutions of \begin{align} \dot V&=-(p+C(u))\,g(u)-p\,(c(u)p+g(u))+2g(u)\,p \\ &=-C(u)\,g(u)-p^2\,c(u) \\ &=-\frac1{ϵ^2}(u^2+3ux_0+3x_0^2-T)u^2-\frac1ϵ(3(x_0+u)^2-T)p^2 \end{align} Assuming $3x_0^2<T$, this is positive around $(u,p)=(0,0)$ and negative for $(u+\frac32x_0)^2>T$ or $|p|$ large enough. Thus you can identify a region with a hole around $(0,0)$ where the system has no stationary points and the vector field flows inwards on the boundary. Now apply Poincaré-Bendixon.

---

How to get this (Lyapunov?) function, start with the usual, multiply with $u'$ and integrate, you get that $$\frac12u'^2+G(u)=-\int c(u)u'^2\,dt.$$ Use partial integration on the right side, $$\int c(u)u'^2\,dt=C(u)u'-\int C(u)u''=C(u)u'+\int C(u)(c(u)u'+g(u))\,dt=C(u)u'+\frac12 C(u)^2+\int C(u)g(u)\,dt.$$ Now insert and combine to complete squares $$ \frac12(u'+C(u))^2+G(u)=-\int C(u)g(u)\,dt $$ The first right side has an integrand that is a multiple of $u'^2$, the second integrand is a multiple of $u^2$. Use the sum of both to get both coordinates involved. Then use handwaving to make the segment where $c(u)<0$ negligible?

---

Appendix: Code for the plots

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

eps=0.1; T=1.0;
def odesys(y): x,b = y; return [ -(b+x**3-T*x)/eps, x-x0 ];

x_lim = [-2,2]; y_lim = [-1,1] 
X, B = np.meshgrid(np.arange(x_lim[0], x_lim[1]+0.01, 0.1), np.arange(y_lim[0], y_lim[1]+0.01, 0.1))

x_cub = np.linspace(x_lim[0], x_lim[1], 501);
b_cub = T*x_cub-x_cub**3;
t_sol = np.linspace(0, 20, 2501);

for k in range(3):
     x0 = [ -0.5, 0, 0.5][k]*T**0.5;
     a= (T/3.0)**0.5 # extrema at -a, +a; -a < x0 < a
     plt.subplot(1,3,k+1);
     plt.ylim(y_lim);
     plt.plot(x_cub, b_cub,'r', lw=2);
     dotX, dotB = odesys([X,B]);
     plt.streamplot(X, B, dotX, dotB, density=2.5);
     y_sol = odeint(lambda y,t: odesys(y), [(T/3)**0.5, 2*(T/3)*(T/3)**0.5], t_sol, atol=1e-9, rtol=1e-10);
     x_sol, b_sol = y_sol.T; plt.plot(x_sol, b_sol, 'k', lw=1.5);
     plt.grid(); plt.title('$\epsilon=%.2g$, $x_0=%.2f$, $T=%.2f$'%(eps,x0,T));
plt.show();