Stability analysis of non-linear system outside of equilibrium point I have a non-linear system that comes from the transformation matrix for angular rotational rate from Euler -> Body frame:
$$
\dot{r} = A + sin(r)tan(p)B + cos(r)tan(p)C \\
\dot{p} = cos(r) B - sin(r) C
$$
where A, B, and C are constants. The equilibrium point is at:
$$
r_e = arctan(B/C)\\
p_e = arctan(-\frac{A}{B sin(r_e) + C cos(r_e)})
$$
When taking the Jacobian matrix and substitute the equilibrium point, I found that the real part of the eigenvalue is $$B cos(r)tan(p) - C sin(r) tan(p)$$
which always equals zero at the equilibrium point, so the system is marginally stable at the equilibrium point. This lead me to thinking that I need to check the system stability outside of the equilibrium point using Lyapunov theory. However, I've been struggling with finding the right candidate function.
When plotting the phase portrait of the system at some A, B, C values, I found that $r$ and $p$ seem to exhibit cyclical behavior like sine and cosine waves. For example, here's the
phase portrait at A = 0, B = 0, C = 1
However, after closer examination, they are not exactly sine and cosine waves. So something like $r^2 + p^2$ for the candidate function wouldn't work. Here's the plot of r^2 + d^2 when simulated with dt = 0.0001. It turns out that the radius is slower growing, so for this choice of A, B, and C, the system is unstable. For A = 0.1, B = 0.1, C = -1.0, the system is also unstable (plot of r^2 + d^2), but going a slower rate. I'm trying to find a candidate function where we can show that its derivative is positive depending on the choice of A, B, and C, but
I'm not sure how I should go about finding this candidate function, and if Lyapunov theory is not the suitable approach here, what other analysis should I look into?
Thanks
 A: This isn't a full answer, as I wasn't able to find a Lyapunov function, but hopefully this might at least might provide some insights into the dynamics of this system.
By performing a coordinate transformation it might lead to simpler dynamics, which allows for more insights of the dynamics of the system. Since common terms in the original dynamics are $\sin(r)$, $\cos(r)$ and $\tan(p)$, I propose the following coordinate transformation
\begin{align}
x &= B\,\sin(r) + C\,\cos(r), \\
y &= B\,\cos(r) - C\,\sin(r), \\
z &= \tan(p).
\end{align}
Which can be shown to have to following dynamics
\begin{align}
\dot{x} &= \left(B\,\cos(r) - C\,\sin(r)\right) \left(A + \left(B\,\sin(r) + C\,\cos(r)\right) \tan(p)\right) \\
&= y \left(A + x\,z\right), \\
\dot{y} &= -\left(B\,\sin(r) + C\,\cos(r)\right) \left(A + \left(B\,\sin(r) + C\,\cos(r)\right) \tan(p)\right) \\
&= -x \left(A + x\,z\right), \\
\dot{z} &= \sec(p)^2 \left(B\,\cos(r) - C\,\sin(r)\right) \\
&= y \left(1 + z^2\right).
\end{align}
By using $\sin(r)^2 + \cos(r)^2 = 1$, it can also be noted that
$$
x^2 + y^2 = R^2,
$$
with $R^2 := B^2 + C^2$. Using this proposed coordinate transformation it can be noted that the original coordinates $r$ and $p$ can only be recovered uniquely when assuming that $-\pi \leq r < \pi$ and  $-\pi/2 \leq p < \pi/2$. The equilibrium points of this transformed system can be shown to be $(x_{eq}, y_{eq}, z_{eq}) = (R, 0, -A/R)$ and $(x_{eq}, y_{eq}, z_{eq}) = (-R, 0, A/R)$. Evaluating the Jacobian at both equilibria yields the same following result
$$
J_{eq} =
\begin{bmatrix}
y\,z & A+x\,z & x\,y \\
-A-2\,x\,z & 0 & -x^2 \\
0 & 1+z^2 & 2\,y\,z
\end{bmatrix}_{eq} = 
\begin{bmatrix}
0 & 0 & 0 \\
A & 0 & -R^2 \\
0 & 1+A^2/R^2 & 0
\end{bmatrix},
$$
which similar to the original system only has eigenvalues for zero real part. So again nothing can be said about the local stability of the equilibria.
If one prefers to have a two instead of three dimensional system, one can reduce this system while introducing back some trigonometric functions. Since $x^2 + y^2 = R^2$ holds and $y=0$ holds for both equilibria a natural choice would be to define the reduced coordinate as $x = -R\cos(\theta)$ and $y = R\sin(\theta)$ and scaling $z$ by $R$, i.e. $\sigma = R\,z$, results in slightly neater dynamics, resulting in
\begin{align}
\dot{\theta} &= A - \sigma \cos(\theta), \\
\dot{\sigma} &= \left(R^2 + \sigma^2\right) \sin(\theta),
\end{align}
with equilibria $(\theta_{eq},\sigma_{eq}) = (2\,k\,\pi, A)$ and $(\theta_{eq},\sigma_{eq}) = ((2\,k+1)\,\pi, -A)$ for all $k \in \mathbb{Z}$. The corresponding Jacobians again have eigenvalues with zero real part, as shown below
$$
J_{eq} =
\begin{bmatrix}
\sigma \sin(\theta) & -\cos(\theta) \\
\left(R^2 + \sigma^2\right) \cos(\theta) & 2\,\sigma \sin(\theta)
\end{bmatrix}_{eq} = 
\begin{bmatrix}
0 & \mp 1 \\
\pm\left(R^2 + A^2\right) & 0
\end{bmatrix}.
$$
I tried to find a Lyapunov function for each of the proposed transformed systems, but without success. I also tried limit my search to just Lyapunov functions whose time derivative is zero (using as starting point $V(x) = (x-x_{eq})^\top P\,(x-x_{eq})$ with $P\succ0$ and $J_{eq}^\top\,P + P\,J_{eq} = 0$), to at least ensure Lyapunov stability, but also without success.
