Given the system
$$\dot{u_1} = \eta_1 u_2 u_3, \dot{u_2} = \eta_2 u_1 u_3, \dot{u_3} = \eta_3 u_2 u_3$$
with constants
$$\alpha = \frac{12}{m(a^2 + c^2)}, \beta = \frac{12}{m(a^2 + b^2)}, \gamma = \frac{12}{m(b^2 + c^2)}$$
$$\eta_1 = (\gamma - \beta), \eta_2 = (\alpha - \gamma), \eta_3 = (\beta - \alpha)$$
without loss of generality, we let $\left|u(0)\right| = 1$ and solved for equilibira, resulting in
$$(\pm1,0,0),(0,\pm1,0),(0,0,\pm1)$$
From here, the stability near equilibria is analyzed by defining $F(u)=\alpha u^2+\beta u^2+\gamma u^2$, where $u=\begin{bmatrix} u_1 & u_2 & u_3\end{bmatrix}$, using a $u(0)$ near each equilibria, and observing the evolution of $F(u)$ on the unit sphere about the chosen equilibrium. For example, choosing a $u(0) \approx (0,0,1)$, we get $F(u) = F(u(0)) = \delta + \gamma$ for some small $\delta>0$. This results in a closed orbit about $(0,0,1)$. Hence, we get $(0,0,1)$ is neutrally stable.
My question is why do we not reduce the problem to a 2-dimensional system by substituting $u_2=\sqrt{1-u_1^2-u_3^2}$ and then analyzing stability by using the eigenvalues of the Jacobian matrix at the equilibria $(u_1,u_3)=(\pm1,0)$ and $(u_1,u_3)=(0,\pm1)$? I used this process to remove a variable from an SIRS problem earlier, and it worked well. However, when I attempted to reduce the dimension here, I ended up with undefined values in the Jacobian, and I'm not sure why.
Also, is it possible to find and use the eigenvalues of a 3-dimensional Jacobian to analyze the stability? If so, why is this process not used either?
The short answer is that one does not reduce the system two a two dimensional one by expressing one variable in terms of the rest, because the qualitative picture of the dynamics is much clearer when things are left in their symmetric form and one analyzes the system on the two sphere. Furthermore, linearization techniques, as the one you suggest, will fail miserably here, especially for four out of the six equilibrium points, because the system is conservative and in fact it is a totally integrable system in every sense of the word - it has a full set of commuting conserved quantities and you can write down the solutions explicitly in terms of elliptic functions or theta function, whatever your preference is.
First of all, the system is called the Euler's top. It describes the rotation of a solid object around it's center of mass in the absence of any other external forces (think of a satellite in space, twisting and turning in zero gravity).
Next, you have two conserved quantities, $J(u) = u_1^2+u_2^2+u_3^2$ and $F(u) = \alpha\, u_1^2 + \beta\, u_2^2 + \gamma\, u_3^2$. The latter, $F(u)$ is the so called Hamiltonian of the system. The conservation of $F(u)$ comes from the conservation of energy. The conservation of $J(u)$ comes from the conservation of angular momentum of the system.
An orbit, starting from a point on any of the level spheres $J(u) = c_1$ stays on that sphere forever. Analogously, any orbit, starting from a point on any of the level surface $F(u) = c_2$ stays on it forever. Thus the curves defined by $J(u) = c_1, \, F(u) = c_2$ are the orbits of your system.
Finally, choose a sphere, say $J(u) = 1$, and study the dynamics on this sphere. The orbits of your system are obtained as the intersection curves of the sphere with $F(u) = c_2$ for various values of $c_2$. The surface $F(u) = c_2$ is an ellipsoid. So an orbit of the system is an intersection curve of the sphere with an ellipsoid. Think of an ellipsoidal balloon that you start inflating homogeneously (i.e. all ellipsoidal balloons are similar, given by the equations $F(u)=c_2$ for various values of $c_2$). Then you star with an ellipsoid inside the sphere, not intersecting it at all (a vary small $c_2$). At some point (for some value of $c_2$) the ellipsoid touches the sphere at two antipodal points for the first time and these two points are the first two equilibrium points of the system (touches means that the sphere and the ellipsoid shear a pair of common tangent planes at these two points). Then, as $c_2$ grows, the intersections become pairs of ovals, encircling the two equilibrium points. These ovals are the orbits of your system, as discussed earlier. The two equilibrium points are stable, but not asymptotically. The function $F(u)$ is then a Lyapunov function and can be used to justify the stability. This process continues until the ellipsoid becomes tangent to two new antipodal equilibrium points (and at the same time intersects the sphere in a pair of curve with merging corners, in other words, the ellipsoid and the sphere share a pair of tangent planes at two antipodal points and at the same time the two surfaces intersect at other points). The intersection curves then become two pairs of heteroclinic orbits connecting the two new equilibira. These equilibria are saddle points. You can successfully apply a linearization procedure, as you suggest, for these two points, concluding non-stability.
As the ellipsoid expands further, the orbits turn again into pairs of ovals, which shrink and shrink as $c_2$ grows, and this happens until the ellipsoid touches the sphere for the last time at a pair of antipodal points (these are the only two intersection points between the sphere and the ellipsoid), which are also stable but not asymptotically, with Lyapunov function $F(u)$.
Edit. The SIRS model you are referring to is probably this one: \begin{align} \dot{S} &= - \alpha \, S\, I\\ \dot{I} &= \alpha \, S\, I - \beta \, I\\ \dot{R} &= \beta \, I \end{align} One can immediately see that the first two equations are not influenced by the variable $R$ which is present only in the third one. Hence, if one solves only the system of the first two equations \begin{align} \dot{S} &= - \alpha \, S\, I\\ \dot{I} &= \alpha \, S\, I - \beta \, I\\ \end{align} the third one is easy to solve, you just have to integrate $I$ one more time to get $R$. Alternatively, you can check that the quantity $F = S + I + R$ is conserved, so your system restricts to the family of planes $ S + I + R = c$ for any real $c$, just like in the previous example, the system restricts to families of spheres (and ellipsoids). The mechanism behind it is similar. Just like in the case of the Euler top, where you had the spheres and the ellipsoids, in this model you have another conserved quantity, i.e. you have the conserved quantity $F = S+I+R$ and there is another one, a bit more complicated, which allows you to completely analyze the qualitative behavior of your original system and in fact to solve it completely (I haven't done the calculations, but I can see more or less see it).