<s>[INST] Euler angles and gimbal lock Can someone show mathematically how gimbal lock happens when doing matrix rotation with Euler angles for yaw, pitch, roll?  I'm having a hard time understanding what is going on even after reading several articles on Google.
Is the only work-around to use quaternions?
 [/INST] Many answers have been given that are mathematically very sound and correct. Here is another way, using only high school mathematics:
Suppose, we use the ZXY Euler rotations: first about the Z-axis (yaw), then about the X-axis (roll) and finally about the y-axis(pitch).
$R_z(\psi) = \left(\begin{array}{ccc}
\cos \left(\psi \right) & -\sin \left(\psi \right) & 0\\
\sin \left(\psi \right) & \cos \left(\psi \right) & 0\\
0 & 0 & 1
\end{array}\right)$
$R_x(\phi) = \left(\begin{array}{ccc}
1 & 0 & 0\\
0 & \cos \left(\phi \right) & -\sin \left(\phi \right)\\
0 & \sin \left(\phi \right) & \cos \left(\phi \right)
\end{array}\right)$
$R_y(\theta) = \left(\begin{array}{ccc}
\cos \left(\theta \right) & 0 & \sin \left(\theta \right)\\
0 & 1 & 0\\
-\sin \left(\theta \right) & 0 & \cos \left(\theta \right)
\end{array}\right)$
The rotation matrix in ZXY configuration would be:
$R = R_z(\psi)R_x(\phi)R_y(\theta)$
$R = \left(\begin{array}{ccc}
\cos \left(\psi \right)\,\cos \left(\theta \right)-\sin \left(\phi \right)\,\sin \left(\psi \right)\,\sin \left(\theta \right) & -\cos \left(\phi \right)\,\sin \left(\psi \right) & \cos \left(\psi \right)\,\sin \left(\theta \right)+\cos \left(\theta \right)\,\sin \left(\phi \right)\,\sin \left(\psi \right)\\
\cos \left(\theta \right)\,\sin \left(\psi \right)+\cos \left(\psi \right)\,\sin \left(\phi \right)\,\sin \left(\theta \right) & \cos \left(\phi \right)\,\cos \left(\psi \right) & \sin \left(\psi \right)\,\sin \left(\theta \right)-\cos \left(\psi \right)\,\cos \left(\theta \right)\,\sin \left(\phi \right)\\
-\cos \left(\phi \right)\,\sin \left(\theta \right) & \sin \left(\phi \right) & \cos \left(\phi \right)\,\cos \left(\theta \right)
\end{array}\right)$
Now, let's replace $\phi =\frac{\pi }{2}$
You will get R as:
$R = \left(\begin{array}{ccc}
\cos \left(\psi \right)\,\cos \left(\theta \right)-\sin \left(\psi \right)\,\sin \left(\theta \right) & 0 & \cos \left(\psi \right)\,\sin \left(\theta \right)+\cos \left(\theta \right)\,\sin \left(\psi \right)\\
\cos \left(\psi \right)\,\sin \left(\theta \right)+\cos \left(\theta \right)\,\sin \left(\psi \right) & 0 & \sin \left(\psi \right)\,\sin \left(\theta \right)-\cos \left(\psi \right)\,\cos \left(\theta \right)\\
0 & 1 & 0
\end{array}\right)$
Using trigonometric identities, we get:
$R = \left(\begin{array}{ccc}
\cos \left(\psi +\theta \right) & 0 & \sin \left(\psi +\theta \right)\\
\sin \left(\psi +\theta \right) & 0 & -\cos \left(\psi +\theta \right)\\
0 & 1 & 0
\end{array}\right)$
There is infinite $\psi$ and $\theta$ that create the same rotation matrix.
For $\phi =-\frac{\pi }{2}$:
$R = \left(\begin{array}{ccc}
\cos \left(\psi \right)\,\cos \left(\theta \right)+\sin \left(\psi \right)\,\sin \left(\theta \right) & 0 & \cos \left(\psi \right)\,\sin \left(\theta \right)-\cos \left(\theta \right)\,\sin \left(\psi \right)\\
\cos \left(\theta \right)\,\sin \left(\psi \right)-\cos \left(\psi \right)\,\sin \left(\theta \right) & 0 & \cos \left(\psi \right)\,\cos \left(\theta \right)+\sin \left(\psi \right)\,\sin \left(\theta \right)\\
0 & -1 & 0
\end{array}\right)$
$R = \left(\begin{array}{ccc}
\cos \left(\psi -\theta \right) & 0 & -\sin \left(\psi -\theta \right)\\
\sin \left(\psi -\theta \right) & 0 & \cos \left(\psi -\theta \right)\\
0 & -1 & 0
\end{array}\right)$
You can use Symbolic math tools to create these simplifications. I have used MATLAB:
syms psi

R_z = [ cos(psi), -sin(psi), 0; 
        sin(psi),  cos(psi), 0;
        0       ,  0       , 1
        ];
syms phi 

R_x = [ 1, 0,           0; 
        0, cos(phi), -sin(phi);
        0, sin(phi),  cos(phi)
        ];

syms theta 

R_y = [  cos(theta), 0, sin(theta); 
         0         , 1, 0;
        -sin(theta), 0, cos(theta)
        ];
R = R_z*R_x*R_y;
R_new_1_expr = subs(R, phi, pi/2)
R_new_1_simplified = simplify(R_new_1_expr);
display(R_new_1_simplified)


$\left(\begin{array}{ccc}
\cos \left(\psi +\theta \right) & 0 & \sin \left(\psi +\theta \right)\\
\sin \left(\psi +\theta \right) & 0 & -\cos \left(\psi +\theta \right)\\
0 & 1 & 0
\end{array}\right)$

</s>