It's not entirely clear (to me) what you mean by "applying an angular velocity". Angular velocity is an instantaneous property; if you only have a single angular velocity vector, you only have the first derivative of the motion at one point in time, and thus can't say anything about the orientation at other times.
I'll assume that you're assuming the angular velocity to be constant, and that by "applying" it you mean calculating the orientation at later times of a body that has the given initial orientation and undergoes a motion with the given constant angular velocity.
To do this, you need to know a) how to find the unit quaternion corresponding to a rotation vector and b) how to compose two unit quaternions in order get the unit quaternion corresponding to the composition of the corresponding rotations.
For a), let's denote the angular velocity vector by $\vec{\omega}$. Over a time $t$, a motion with this (constant) angular velocity results in a rotation that can be described by the vector $\vec{w}=\vec{\omega}t$, where the direction and magnitude of $\vec{w}$ specify the axis and angle of the rotation, respectively. This rotation corresponds to the unit quaternion
$$q=\left(\cos\frac{|\vec{w}|}{2},\frac{\vec{w}}{|\vec{w}|}\sin\frac{|\vec{w}|}{2}\right)\;.$$
For b), since applying the rotation corresponding to a unit quaternion $q$ to a three-dimensional vector $\vec{v}$ consists in conjugating a quaternion $\nu$ corresponding to $\vec{v}$ by $q$, composing rotations just corresponds to multiplying the corresponding quaternions:
$$q(\vec{v})=q\nu q^{-1}\;,$$
$$q_1(q_2(\vec{v}))=q_1(q_2\nu q_2^{-1})=q_1q_2\nu q_2^{-1} q_1^{-1}=(q_1q_2)(\vec{v})\;.$$
(This is part of what makes quaternions so useful for representing three-dimensional rotations). So now you just multiply $q$ and your initial unit quaternion, say, $q_0$. Since the initial rotation is applied first, the corresponding unit quaternion is on the inside of the nested conjugations, so it has to be the right-hand factor: $q'=qq_0$.
For more on all this, see this Wikipedia article.
After a lot of searching, I have finally found a reference that has what I needed. Although it doesn't appear in a 'Quaternion to Euler' or 'Quaternion to Tait-Bryan' google search, when I gave up and started looking for 'Quaternion to Axis-Angle' with the intention of going through that representation as an intermediate step, I came across the following wonderful document:
Technical Concepts
Orientation, Rotation, Velocity and Acceleration,
and the SRM
Version 2.0, 20 June 2008
Author: Paul Berner, PhD
Contributors: Ralph Toms, PhD, Kevin Trott, Farid Mamaghani, David Shen,
Craig Rollins, Edward Powell, PhD
http://www.sedris.org/wg8home/Documents/WG80485.pdf
It covers a lot of the formalisms, but most importantly, shows derivations and solutions for 3-1-3 and 3-2-1 Euler angle representation. It also seems to cover inter-conversion between pretty much every other rotation representation I'm aware of, and so I would also recommend it as a good general reference.
Oh, and the actual solution for a 3-2-1 ($z-y-x$) Tait-Bryan rotation convention from that reference:
$$
\phi = \operatorname{arctan2}\left(q_2 q_3 + q_0 q_1,\frac{1}{2}-(q_1^2 + q_2^2)\right) \\
\theta = \arcsin(-2(q_1 q_3 - q_0 q_2)) \\
\psi = \operatorname{arctan2}\left(q_1 q_2 + q_0 q_3,\frac{1}{2}-(q_2^2 + q_3^2)\right)
$$
Note that the gimbal-lock situation occurs when $2(q_1 q_3 + q_0 q_2) = \pm1$ (which gives a $\theta$ of $\pm \frac{\pi}{2} $), so it can be clearly identified before you attempt to evaluate $\phi$ and $\psi$.
(Convention for arctan2 is $\operatorname{arctan2}(y, x)$, as hinted on page 3.)
Best Answer
You have a curve $q\colon I \to \mathbb S^3$, where $I$ is a time interval and $\mathbb S^3\subseteq\mathbb R^4$ is the set of all unit quaternions. Here, we represent quaternions as 4-dimensional vectors. Note that if a quaternion should encode an orientation or rotation, it has to be a unit quaternion.
Since $q(t)\in\mathbb S^3\subseteq \mathbb R^4$ lives in the linear space $\mathbb R^4$ you can calculate its time derivative $\dot q(t) = \tfrac d{dt}q(t)$ by \begin{align}\tag{1} \dot q(t) = \tfrac d{dt}q(t) = \lim_{h\to0} \frac{q(t+h) - q(t)}h, \end{align} meaning that $\bigl(q(t+\delta)-q(t)\bigr)/\delta$ is an approximation to $\dot q(t)$ for small $\delta$. We can see that this approximation lives somewhere in $\mathbb R^4$. More specifically, it is an element of $T_{q(t)}\mathbb S^3$, the tangent space of the sphere at the element $q(t)\in\mathbb S^3$, because $q(\tau)\in\mathbb S^3$ for all $\tau\in I$.
Since $\mathbb S^3$ is a Lie group it is favourable to represent the velocity of $q(t)$ by a vector $\Omega(t)\in\mathbb R^3$ that fulfills $$\dot q(t) = \frac12 q(t) * \begin{bmatrix}0\\\Omega(t)\end{bmatrix}.$$ This can be thought of as mapping the tangent space $T_{q(t)}\mathbb S^3$ to the tangent space $T_e\mathbb S^3 = \{[0,x^T]^T\in\mathbb R^4\}$ with the neutral element $e=[1,0,0,0]^T$. The $\Omega(t)$ is actually the angular velocity in the body frame. We can calculate $\Omega(t)$ from $\dot q(t)$ by $$\tag{2} \Omega(t) = \operatorname{Im}\bigl(2 \overline{q(t)}*\dot q(t)\bigr), $$ where the overline represents quaternionic conjugation (which is actually the inversion on $\mathbb S^3$) and $\operatorname{Im}$ extracts the imaginary part, hence just drops the first component (which has to be zero here). Now we can put the limit expression for $\dot q(t)$ in here an get \begin{align*} \Omega(t) &= \operatorname{Im}\bigl(2 \overline{q(t)}*\dot q(t)\bigr)\\ &= \operatorname{Im}\left(2 \overline{q(t)}*\lim_{h\to0} \frac{q(t+h) - q(t)}h\right) \\ &= \operatorname{Im}\left(2 \lim_{h\to0} \frac{\overline{q(t)}*q(t+h) - \overline{q(t)}*q(t)}h\right) \\ &= \operatorname{Im}\left(2 \lim_{h\to0} \frac{\overline{q(t)}*q(t+h) - e}h\right) \\ &= \lim_{h\to0} 2\operatorname{Im}\frac{\overline{q(t)}*q(t+h)}h, \end{align*} since $\operatorname{Im}e = [0,0,0]^T$. This means that $2\operatorname{Im}\overline{q(t)}*q(t+\delta)/\delta$ is an approximation to $\Omega(t)$, the body frame angular velocity, for small $\delta$.
With equations (1) and (2) from above, you can transform the derivative $\dot q(t) = \tfrac d{dt}q(t)$ into the body frame angular velocity $\Omega(t)$ and back.
Note that $2\operatorname{Im}q(t+\delta)*\overline{q(t)}/\delta \approx \omega(t)$ gives an approximation to the angular velocity with respect to the inertial frame. It simply holds $\omega(t)=\Omega(t)^{q(t)}$, where $x^{q(t)}\in\mathbb R^3$ is the application of the rotation $q(t)$ to $x$ and is defined by $$ \begin{bmatrix}0\\ x^{q(t)}\end{bmatrix} = q(t)*\begin{bmatrix}0\\ x\end{bmatrix}*\overline{q(t)}. $$
The simplest way to numerically integrate the attitude is to use a forward Lie-group Euler method. If you somehow calculate $\dot q(t)$ or $\Omega(t)$ you can get the new attitude $q(t+\delta)$ by $$ q(t+\delta) = q(t)*\widetilde{\exp}(\delta\cdot \Omega(t)), $$ where $\widetilde{\exp}\colon\mathbb R^3\to\mathbb S^3$ is essentially a Lie group exponential function which is defined by a power series, but can, for unit quaternions, be written in closed form as $$ \widetilde{\exp}(v) = \cos(\tfrac12 \|v\|) + \frac{v}{\|v\|}\sin(\tfrac12 \|v\|) $$ for $v\in\mathbb R^3$, where $\|\bullet\|$ is the standard Euclidean norm.
If you're interested in more sophisticated Lie group time integration, you can check out one of my papers [Arnold, Hante 2016] or a nice paper (without quaternions though) of my supervisor [Arnold, Cardona, Brüls 2016] or the preliminary version.
Hope I could clear things up a little.