Geodesics on SO(3) I have two 3D rotations about the origin, represented as
$3 \times 3$ orthogonal matrices $M_1$ and $M_2$
(specified by numerical entries),
and I would like to interpolate (and compute)
a continuous sequence of rotations
between $M_1$ and $M_2$.
Ideally this interpolation would follow a geodesic on SO($3$).
Assuming this is well-known, 
(Q1): I would appreciate a pointer to geodesics on SO($3$),
especially suited to my computational task.
(Q2) [added]. I wonder if the geodesics identified by 
Benjamin,
Paul Siegel,
and Francois Ziegler,
which effectively answer Q1,
can be viewed as equivalent to following a great-circle arc on $S^3$ between
unit quaternions representing the two rotations?
Update (25Jan2019): A paper was just published on this
topic: 

Novelia, Alyssa, and Oliver M. O’Reilly. "On geodesics of the rotation group $SO(3)$." Regular and Chaotic Dynamics 20, no. 6 (2015): 729-738. Springer link.

Update (26Jan2019): The paper above raises the interesting question
whether saccadic motions of the eye follow $SO(3)$ geodesics, 
an issue possibly dating back
to Helmholtz's investigations in the late $19^{\mathrm{th}}$ century:

von Helmholtz, H., "Über die normalen Bewegungen des menschlichen Auges," Archiv für Ophthalmologie, 1863, Vol. 9, No. 2, pp. 153–214; DOI: 10.1007/BF02720895.

 A: If you are assuming the use of the bi-invariant metric, then the geodesics are right/left translations of the one parameter sub-groups $O(t)=O(0) \exp(tG)$ where $G \in \mathfrak{so}(3)$ (i.e. it is anti-symmetric and traceless) and $O(0)\in SO(3)$.
The curve you want will have $O(0)=M_1$ and thus be of the form $O(t) = M_1 \exp(tG)$ with $G=\log(M_1^{-1}M_2)$ (the matrix log). You need to choose a specific matrix log as generically there will be many, you many wish to read about principal matrix logs if you are going to do any numerical computations. This is what will likely come out of Matlab etc. The curve is parameterised by $t \in [0,1]$, it is smooth, is a geodesic of the aforementioned metric and has the end points you hoped for.
A: The geodesics (of the unique connection invariant under left and right translations and inversion) are the translates of one-parameter subgroups: see Helgason, exerc. 6(iii), p. 148. So given your $M_1$ and $M_2$, find (by surjectiveness of $\exp$) a $Z\in\mathfrak{so}(3)$ such that $\exp(Z)=M_1^{-1}M_2$, and put $g(t) = M_1\exp(tZ)$.
