In a certain sense, @Geoff Oxberry is correct in saying that stability and preservation of quadratic invariants are not directly related. For instance, there exist explicit methods that will preserve energy for your problem (and they are certainly not $A$-stable).
However, in another sense there is a relation between the two.
Well-posedness
First, note that your problem as described is not well posed, and the invariant you give does not follow from the problem description. I believe you should have instead
$$\ddot{u} + u = 0$$
with initial conditions
$$u(0) = u_0$$
and
$$\dot{u}(0) = v_0$$
where $u_0^2 + v_0^2 = 1$. Then the problem is well-posed and the invariant condition you have written follows.
Conservation of energy
Now, if you rewrite this second-order equation as a first-order system
$$\begin{pmatrix} \dot{u} \\ \dot{v} \end{pmatrix} =
\begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix}
\begin{pmatrix} u \\ v\end{pmatrix},$$
you find that the eigenvalues are $\pm i$. If we write the system above as $\dot{w} = A w$, the exact solution is given by
$$w(t) = \exp(At)w(0)$$
and the reason you have energy conservation is that $|e^{\lambda t}|=1$ for both eigenvalues of $A$.
Discrete conservation
After discretizing with a one-step method, you will have instead an equation of the form
$$w_{n+1} = R(\Delta t A)w_n,$$
where $\Delta t$ is the numerical time step and $R(z)$ is the stability function of the method (generally a polynomial or rational function). Since $A$ is normal, to have discrete conservation you need that $$|R(\pm i \Delta t)|=1.$$.
The connection
One way to interpret the condition just stated is to say that the eigenvalues of the ODE system, scaled by the numerical timestep, must lie on the boundary of the absolute stability region of the numerical integrator. That's because the stability region is defined precisely as the set of complex numbers $z$ for which $|R(z)| \le 1$.
Thus, with just the right time step, even the 4th-order Runge-Kutta method can be used to preserve energy for this problem, since its stability region boundary intersects the imaginary axis. Of course nobody recommends this, because it only works for one special value of $\Delta t$.
The trapezoidal rule conserves energy for any time step, because its region of absolute stability is exactly the left half-plane; i.e., the entire imaginary axis is on its stability region boundary, so the condition above is fulfilled for any $\Delta t$.
Meanwhile, it's impossible to conserve energy with forward or backward Euler, since those methods have a stability region whose boundary only touches the imaginary axis at the origin.
