Con el comando "pplane" de MATLAB se puede dibujar un retrato de fase bastante bien. Otra forma es esta bonita herramienta en línea: http://comp.uark.edu/~aeb019/pplane.html
En tu caso obtienes esta bonita trama
Puedes ver que tus puntos fijos estables no son asintóticamente estables. Esto es bastante típico de los sistemas hamiltonianos, ya que estas órbitas alrededor de un punto crítico estable se comportan periódicamente como "círculos". Esto se puede ver como el Hamiltoniano te da una función de Lyapunov que te da estabilidad pero NO estabilidad asintótica.
Por ejemplo, tomemos el punto fijo $(1,-\pi/2)$ como se ve en el gráfico. Es estable pero no asintóticamente estable ya que las otras órbitas en la vecindad no se mueven hacia este punto fijo.
Tenemos $$ Df(p,q) = \begin{pmatrix} -\cos(q) & p \sin(q) \\ 1 & \cos(q) \end{pmatrix} $$
y en el equilibrio $(1,-\pi/2)$ obtenemos
$$ Df\left(1,-\frac{\pi}{2}\right) = \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix} $$ con valores propios $\pm i$ . Se trata de puntos fijos no hiperbólicos, ya que su parte real es cero y la estabilidad no puede observarse inmediatamente. Para ello introducimos el concepto de funciones de Lyapunov. Pero antes veamos algunos puntos fijos inestables. Por ejemplo como se ve en el gráfico $(0,0)$ . Usted tiene $$ Df\left(0,0\right) = \begin{pmatrix} -1 & 0 \\ 1 & 1 \end{pmatrix} $$ y por tanto los valores propios $-1$ y $1$ . Dado que un valor propio tiene parte real positiva obtenemos que $(0,0)$ es inestable. Usted ve la estabilidad de los equilibrios hiperbólicos son bastante fáciles de calcular.
Ahora, volvamos a nuestro punto fijo no hiperbólico.
Teorema (funciones de Lyapunov): Sea $\dot x=f(x)$ s.t. $f$ es un $C^1$ campo vectorial und $f(x_e)=0$ dado. Sea $U \subset \mathbb{R}^n$ sea un nbh de $x_e$ . Si una función $V \in C^1(U, \mathbb{R})$ e \begin{equation} \begin{cases} V(x_e)=0 \\ V(x)>0 &\text{for } x \in U\backslash\{x_e\} \\ \left\langle \nabla V(x),f(x)\right\rangle = 0 &\text{for } x \in U\backslash\{x_e\} \end{cases} \label{lyapunov} \Fin. entonces $x_e$ es estable en el sentido de Lyapunov pero no asintóticamente estable. Llámese $V$ Función de Lyapunov.
Por lo tanto, ya se puede ver que los hamiltonianos son una buena elección para las funciones de Lyapunov. Sea $f(p,q)=(\dot p, \dot q)$ y $f(p_e,q_e)=0$ . El Hamiltoniano tiene la propiedad útil $$\left\langle \nabla H(p,q),f(p,q) \right\rangle = \left \langle \binom{H_{p}(p,q)}{H_{q}(p,q)}, \binom{\dot p}{\dot q} \right \rangle=0$$ es decir, se cumple automáticamente la tercera propiedad para una función de Lyapunov. Para una constante $C \in \mathbb{R}$ y $$\nabla(H(p,q)+C)=\nabla H(p,q),$$ y por lo tanto $C$ es arbitraria. Puede elegirse de forma que $H(p_e,q_e)=0$ Formulemos, pues, un teorema que dé la última propiedad que falta para una función de Lyapunov.
Teorema: Sea $f(p,q)=(\dot p,\dot q)=(H_q(p,q),-H_p(p,q))$ y $f(p_e,q_e)=0$ . Si $H$ tiene un mínimo local estricto en $(p_e,q_e)$ entonces $(p_e,q_e)$ es estable en el sentido de Lyapunov pero no asintóticamente estable.
Así pues, basta con comprobar si el equilibrio es un mínimo local estricto de $H$ por lo tanto $H(p,q) > 0$ para $(p,q) \in U \backslash \{(p_e,q_e) \}$ es decir, la segunda propiedad.
Así que tomemos $H(x,y)=\frac{p^2}{2}+p\sin(q)+C$ y elija $C$ tal que $H(1,-\pi/2)=\frac{1}{2}-1+C=0$ es decir $C=\frac{1}{2}$ .
Desde $$\nabla H\left(1,-\frac{\pi}{2}\right)=\binom{p+\sin(q)}{p\cos(q)}\bigg|_{(1,-\frac{\pi}{2})}=\binom{0}{0} \hspace{0.1cm}$$ $$\text{D}^2 H\left(1,-\frac{\pi}{2}\right)=\begin{pmatrix} 1 & \cos(q) \\ \cos(q) & -p\sin(q) \end{pmatrix}\bigg|_{(1,-\frac{\pi}{2})}=\begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}$$
y por tanto su hessiano es positivodefinitivo y $(1,-\pi/2)$ es un mínimo local estricto. Según nuestro teorema anterior $(1,-\pi/2)$ es un equilibrio estable (pero no asintóticamente estable) de su sistema.
Ya que he mencionado que esto es típico de los sistemas hamiltonianos, también quiero mencionar el famoso Duffing-oscilador. Está descrito por el sistema
$$\binom{\dot x}{\dot y}=\binom{ y}{x-x^3} $$
con el Hamiltoniano $H(x,y)=\frac{1}{2} y^2 - \frac{1}{2} x^2 + \frac{1}{4} x^4$ . Obtendrá $(\pm 1,0)$ son puntos fijos estables (pero no asimétricos) y $(0,0)$ es inestable (la inestabilidad se deduce inmediatamente observando los valores propios del Jacobiano).