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.
The space of unit quaternions is isomorphic to $SO(3)$ (the space of possible rotations in 3D space), plus a double cover (which is not important in this context). It's easier to picture a sphere for visualization purposes, though. Now imagine you want to change the quaternion vector (which is a vector from the origin to the surface of the (hyper-)sphere). Since you want to keep the length of the vector the same, you have to move it orthogonally to its direction (for a sphere, this means up/down and left/right, but not in/out of the sphere). Thus, you need a vector whose dot product with the $orientation$ quaternion is zero.
But pure imaginary quaternions are orthogonal to pure real ones, and the real axis represents the identity rotation (no rotation at all). Additionally, if $q_1$ and $q_2$ are quaternions, $u$ is a unit quaternion, and $q_1\cdot q_2=0$, then $(q_1u)\cdot(q_2u)=0$ as well (this is because under the interpretation of quaternions as rotations, they are isometries, so they preserve distance and angles), so $orientation\cdot (q_i\ orientation)=0$, if $q_i$ is any pure-imaginary quaternion.
It's a bit more work that the specific choice $q_i=\frac{\Delta t}2(\omega_x{\bf i}+\omega_y{\bf j}+\omega_z{\bf k})$ gives you the amount and axis of rotation you want in your mapping from quaternions to actual rotations, but accepting that, the sum $orientation+(q_i\ orientation)$ gives you a new vector, which is the same length as the old one (to first order), but pointing in a different direction, which is what you wanted. Since there is a second order term that will tend to increase the distance to the the center (which should be intuitively obvious; imagine moving tangent to a sphere from a point on the surface—you stay near the surface for a while, but as you continue straight, you will drift away from the surface), you presumably want a step after those shown to renormalize the vector back onto the sphere.
The answer I just gave is more a programmer's answer than a mathematician's, since your question seemed quite concrete and grounded in some numerical computation. That said, the math angle here has to do with the isomorphism from the space of unit quaternions to the group $\operatorname{Spin}(3)$, which is the simply connected double cover of $SO(3)$ (as I mentioned. Moreover, the Lie algebra $\mathfrak{so}(3)$ generated by (either of) these groups is exactly the space of skew-symmetric $3\times3$ matrices $[\omega]_\times$ for all possible 3-vectors $\omega$, which is exactly what you are doing with your own angular velocity vector. (Note that $[\omega]_\times$ is the 3D Hodge star, defined here.)
Best Answer
This question might go better over at the GameDev SE site; at first glance, though, this formula:
is not representing what's in the comment above it. The formula you're looking for is $\omega(t) = 2q'(t)\bar{q}(t)$ (where as typical I've written $q'(t) = \frac{dq(t)}{dt}$), but what's written in that code is $\omega(t) = 2q'(t)\bar{q'}(t)$; in other words, you're not multiplying by the conjugate of your original orientation quaternion (as you should be) but by the conjugate of (the approximation of) the derivative. The line you want should be something like
but note that this is all predicated on another assumption - that gyroQuater and boxQuater represent the orientation of the same object at two nearby points in time. From your description it sounds like this might not be quite the case, and if it's not so then you may have to be more explicit about just what behavior you're after.