Taking the equation for simple harmonic motion:
$$
x'' = -\omega^2x
$$
We can introduce the derivative operator:
$$
D = \frac{d}{dt}
$$
This gives us:
$$
x'' =
\frac{d^2\vec{x}}{dt^2} =
\frac{d}{dt} \Big( \frac{d}{dt} \big( \vec{x} \big) \Big) =
D \Big( D \big( \vec{x} \big) \Big) =
D^2 \vec{x}
$$
$$
\implies
D^2 \vec{x} = -\omega^2\vec{x}
$$
This looks like a matrix, so let's see what happens in 2D:
$$
D =
\begin{pmatrix}
a & b \\
c & d
\end{pmatrix}\in \Re
$$
$$
\implies
D^2 =
\begin{pmatrix}
a^2 + bc & b(a + d) \\
c(a + d) & d^2 + bc
\end{pmatrix}
=
\begin{pmatrix}
-\omega^2 & 0 \\
0 & -\omega^2
\end{pmatrix}
$$
Solving \(b(a+d) = 0\) gives us either \(b = 0\) or \(d = -a\).
If we assume that \(b = 0\) then this implies \(a^2 + bc = a^2 = -\omega^2\) which cannot be solved in the reals, so
it must be that \(d = -a\):
$$
D =
\begin{pmatrix}
a & b \\
c & -a
\end{pmatrix}
$$
We still have \(a^2 + bc = -\omega^2\). This only has 2 degrees of freedom (you can only pick 2 of \(a,b,c\) and the
third will be fixed) so lets try rewriting it:
$$
a^2 + bc = -\omega^2 \implies bc = -(\omega^2 + a^2)
$$
Let \(a^2 = \omega^2 (p^2 - 1), |p| \ge 1, \text{sgn}(a) = \text{sgn}(p) \):
$$
bc = -(\omega^2 + a^2) = -p^2\omega^2
$$
$$
\implies b = pq\omega, c = -p\omega/q, q \ne 0
$$
$$
\implies
D_{p,q} =
\omega
\begin{pmatrix}
\text{sgn}(p) \sqrt{p^2 - 1} & pq \\
-p/q & -\text{sgn}(p) \sqrt{p^2 - 1}
\end{pmatrix}
$$
We can see that \(p = q = 1\) gives us the rotation matrix that we expect to see in the perfectly circular case
since we then have that the velocity is 90 degress out of phase with the position, which is itself 90 degress out of
phase with the acceleration:
$$
D_{1,1} =
\begin{pmatrix}
0 & \omega \\
-\omega & 0
\end{pmatrix}
$$
This can be simplified using the identity \( \sec^2(\phi) = \tan^2(\phi) + 1 \):
$$
p = k\sec(\phi), k = \text{sgn}(p)
$$
$$
\implies
D_{p,q} =
\omega
\begin{pmatrix}
k\tan(\phi) & q\sec(\phi) \\
-\sec(\phi)/q & -k\tan(\phi)
\end{pmatrix}
$$
And here we have the general solution for the instantaneous velocity and acceleration of a point undergoing simple
harmonic motion in 2D. Note that \(k\) can be folded into \( \tan(\phi) \) since \( \tan(-\phi) = -\tan(\phi) \) but
\( \sec(-\phi) = \sec(\phi) \).
The demo bit
Since we have a closed form solution for the general case, here's an interactive demo. \(k,\phi,q,\omega\) are as
above, and \(\theta,r\) define the offset of the orbiting body (the black line). \(T\) defines how many seconds to
simulate into the future (the blue line). The red and green lines show the instantaneous velocity and acceleration
respectively. You can click to place the point at an arbitrary position.