Lyapunov Function for Competition Equations I am trying to show that the function 
$$Q(x,y)=be(x-\bar{x})^2+2ce(x-\bar{x})(y-\bar{y}) +cf(y-\bar{y})^2$$ is a Lyapunov function for the competition equations:
$$\dot{x}=x(a-bx-cy)$$
$$\dot{y}=y(d-ex-fy).$$
I can find:
$$(\bar{x},\bar{y})=\Big(\frac{af-cd}{ce-bf},\frac{bd-ae}{bf-ce}\Big). $$
I know that for a function, $V(x,y)$ to be Lyapunov, it must be $$V(x,y) \geq 0, \; \forall x,y \in \mathbb{R},$$ how can this be seen? and that $$\dot{V} \leq 0 \; \forall x,y \in \mathbb{R}.$$
Am I missing a trick? I am just getting a lot of messy algebra. I cannot seem to get $Q$ into any kind of addition of squared quantities, nor can i figure out a way to make the derivative negative for all point in the plane.
$$Q(x,y)=be(x-\bar{x})^2+2ce(x-\bar{x})(y-\bar{y}) +cf(y-\bar{y})^2$$
$$\dot{Q}(x,y)=2\Big[\big(be(x-\bar{x})+ce(y-\bar{y})\big)\dot{x}+\big(ce(x-\bar{x}) +cf(y-\bar{y})\big)\dot{y}\Big]$$
 A: The dynamics can also be written as
$$
\dot{z} = \text{diag}(z) (v - M\,z), \tag{1}
$$
with
$$
z = 
\begin{bmatrix}
x \\ y
\end{bmatrix}, \quad
v = 
\begin{bmatrix}
a \\ d
\end{bmatrix}, \quad
M = 
\begin{bmatrix}
b & c \\ e & f
\end{bmatrix}.
$$
From this it is also easier to see what the equilibria ($\dot{z}=0$) are, namely $\bar{z}_1=0$ or $\bar{z}_2 = M^{-1}v$. The second equilibrium actually gives
\begin{align}
\bar{x}_2 &= \frac{a\,f - c\,d}{b\,f - c\,e}, \tag{2a} \\
\bar{y}_2 &= \frac{b\,d - a\,e}{b\,f - c\,e}, \tag{2b}
\end{align}
which differs in a minus sign in front of $\bar{x}$ compared with your equilibrium. The candidate Lyapunov function can also be written into a "more compact" form
\begin{align}
Q(x,y) &= 
\begin{bmatrix}
x - \bar{x} \\ y - \bar{y}
\end{bmatrix}^\top
\begin{bmatrix}
be & ce \\
ce & cf
\end{bmatrix}
\begin{bmatrix}
x - \bar{x} \\ y - \bar{y}
\end{bmatrix} \tag{3a} \\
&= (z - \bar{z})^\top P\,(z - \bar{z}) \tag{3b}
\end{align}
with
$$
P = \begin{bmatrix}
b\,e & c\,e \\
c\,e & c\,f
\end{bmatrix} \tag{4}.
$$
By defining error coordinates as $\varepsilon = z - \bar{z}_2$, then its dynamics can be written as
\begin{align}
\dot{\varepsilon} &= \text{diag}(z) (v - M\,z) \tag{5a} \\
&= \text{diag}(z) (v - M\,(\varepsilon + M^{-1}v)) \tag{5b} \\
&= -\text{diag}(z)\,M\,\varepsilon \tag{5c}
\end{align}
This allows one to simplify the candidate Lyapunov function even further to
$$
Q(x,y) = \varepsilon^\top P\,\varepsilon \tag{6}.
$$
The constraint that this should satisfy $Q(x,y) \geq 0\ \forall\,x,y\in\mathbb{R}$ is equivalent to that all the eigenvalues of $P$ need to be non-negative. However, it can be noted that $Q(x,y)=0$ would also satisfy this requirement. So instead I will use that $Q(x,y) > 0\ \forall\,\varepsilon\neq0$. This is then equivalent to that all the eigenvalues of $P$ need to be positive. The expression for the eigenvalues are a bit cumbersome. Since $P\in\mathbb{R}^{2\times 2}$ we have $\det(P) = \lambda_1\,\lambda_2$ and $\text{tr}(P) = \lambda_1 + \lambda_2$, with $\lambda_1$ and $\lambda_2$ the two eigenvalues of $P$. So instead of looking at the eigenvalues directly we can also use that $\det(P) > 0$ and $\text{tr}(P) > 0$
\begin{align}
\det(P) &= c\,e\,(b\,f - c\,e) > 0, \tag{7a} \\ 
\text{tr}(P) &= b\,e + c\,f > 0. \tag{7b}
\end{align}
By using $(5c)$ and $(6)$ the time derivative of the candidate Lyapunov function can be expressed as
\begin{align}
\dot{Q}(x,y) &= \dot{\varepsilon}^\top P\,\varepsilon + \varepsilon^\top P\,\dot{\varepsilon} \tag{8a} \\
&= -\left(\text{diag}(z)\,M\,\varepsilon\right)^\top P\,\varepsilon - \varepsilon^\top P\,\text{diag}(z)\,M\,\varepsilon \tag{8b} \\
&= -\varepsilon^\top \left(M^\top \text{diag}(z)\,P + P\,\text{diag}(z)\,M\right) \varepsilon \tag{8c}
\end{align}
Checking whether this satisfies the constraint that $\dot{Q}(x,y) \leq 0\ \forall\,x,y\in\mathbb{R}$ becomes a bit more cumbersome. Namely one could look at whether the eigenvalues of $M^\top \text{diag}(z)\,P + P\,\text{diag}(z)\,M$ are always non-negative, but this would be a function of $z$ and thus also $\varepsilon$ so even if an eigenvalue is negative for a given $x$ and $y$ it might still be that $\dot{Q}(x,y) \leq 0$. Only if both eigenvalues are negative then $\dot{Q}(x,y) \leq 0$ would be false. So stability is unsure if only the determinant or the trace (so not both at the same time) is positive. If they both are non-negative it is stable, but when evaluating the determinant and trace it can be seen that that is not the case
\begin{align}
\det(M^\top \text{diag}(z)\,P + P\,\text{diag}(z)\,M) &= 4\,c\,e\,x\,y\,(b\,f - c\,e)^2 > 0, \tag{9a} \\ 
\text{tr}(M^\top \text{diag}(z)\,P + P\,\text{diag}(z)\,M) &= 
2 \left(e\,x\,(b^2 + c^2) + c\,y\,(e^2 + f^2)\right) > 0, \tag{9b}
\end{align}
since for example depending on the sign of $c\,e$ a different pair of opposing quadrants for $x$ $y$ will yield that the determinant is negative. But maybe since you are dealing with a competitive Lotka–Volterra equation you are only interested in $x,y\in\mathbb{R}_{\geq0}$. But I am going to leave it here, so hopefully this answer has given you enough insights and tricks to help to proceed.
