When (and why) does a linear vector field enjoy polynomial conserved quantities? The familiar equations $\sin^2 + \cos^2 = 1$ and $\cosh^2 - \sinh^2 = 1$ can be explained by noticing that the corresponding 2-dimensional linear ODE have particularly simple conserved quantities; if $(x, y)$ are time-dependent variables solving $x' = y$ and $y' = -x$, the polynomial $x^2 + y^2$ is constant over time, while if $x' = y$ and $y' = x$, then $x^2 - y^2$ is constant. Checking that such identities hold is simple, but where do they really come from? A tentative explanation might go as follows. 
Suppose that $x$ and $y$ are implicitly related along a curve in $\mathbb{R}^2$ tangent to our vector field $V = a \,d/dx + b \,d/dy$. Now, consider the "orthogonal" 1-form $w = b \,dx - a \,dy$. The tangency condition is equivalent to asserting that the pullback of $w$ to the curve vanishes. Thus, when $w$ is an exact form on $\mathbb{R}^2$, an integral will give a conserved quantity for our vector field.
When $V$ is a linear vector field given by a matrix $A$, the 1-form $w$ is, in the natural bases, expressed as
$$W = \begin{bmatrix} 0 & 1 \\ -1 & 0 \end{bmatrix} A$$
Our 1-form is closed exactly when this matrix is symmetric---that is, the trace of $A$ is zero---in which case an integral is given by the homogeneous quadratic polynomial $1/2 x^T W x$.
This ad-hoc procedure is apparently able to explain the appearance of conic sections as the solutions to linear ODE, as well as the fact that the integral curves to some (indeed, most) linear ODE are not closely related to the zero sets of polynomials. However, this story does not seem entirely resolved; for example, we have not considered whether quadratic polynomial invariants can be discovered in other ways. My question is: in general, how does one detect the whole collection of polynomial invariants associated with an n-dimensional linear vector field?
 A: The polynomial $P(X)$ is invariant under the autonomous linear differential equation
$X' = A X$
iff 
$ \nabla P(X) A X = 0$ for all $X$.
Of course the terms of each total degree must be zero, so we may as well consider the case where $P(X)$ is homogeneous of a particular degree $d$.
The case $d=0$ is trivial ($P$ is constant).
Next comes $d=1$, where $\nabla P(X) = b$ is a constant row vector, and $b A = 0$ means $b$ is orthogonal to the range of $A$. Thus if $A$ has less than full rank, the motion is confined to affine subspaces parallel to the range of $A$.
Next is $d=2$, where $\nabla P(X) = X^T B$ for some symmetric matrix $B$, and we need $X^T BA X = 0$, which is equivalent to $BA + (BA)^T = B A + A^T B = 0$ (i.e. $BA$ is an antisymmetric matrix).
I don't have a general classification of the case $d=3$, but here is a randomish example with $n=3$.  $$P(x_1, x_2, x_3) = \left( x_{{3}}+x_{{1}} \right)  \left( 2\,{x_{{1}}}^{2}+x_{{1}}x_{{2}
}+x_{{2}}x_{{3}}-{x_{{3}}}^{2} \right) 
$$
is invariant for $$A = \left[ \begin {array}{ccc} a_{{1,1}}&0&2\,a_{{1,1}}+a_{{3,3}}
\\ -8\,a_{{1,1}}-4\,a_{{3,3}}&-4\,a_{{1,1}}-4\,a_{{3
,3}}&2\,a_{{1,1}}+4\,a_{{3,3}}\\ a_{{1,1}}+2\,a_{{3,
3}}&0&a_{{3,3}}\end {array} \right]
$$
