<p>Here is the same first order Taylor series method as given in minorlogic's answer (using the JOML API):</p>
<pre><code>public static Quaterniond integrateAngularVelocity(
        double avx, double avy, double avz, double dt) {
    double len = Math.sqrt(avx * avx + avy * avy + avz * avz);
    double theta = len * dt * 0.5;
    if (len &gt; 1.0e-12) {
        double w = Math.cos(theta);
        double s = Math.sin(theta) / len;
        return new Quaterniond(avx * s, avy * s, avz * s, w);
    } else {
        return new Quaterniond(0, 0, 0, 1);
    }
}
</code></pre>
<p>The origin and meaning of this method is explained in the background section of <a href="https://www.researchgate.net/publication/257451700_A_novel_Quaternion_integration_approach_for_describing_the_behaviour_of_non-spherical_particles" rel="nofollow noreferrer">A novel Quaternion integration approach for describing the behaviour of non-spherical particles</a> by Facheng Zhao and Berend van Wachem (2013), and in significant detail in sections 4.5 (Time Derivatives) and 4.6 (Time-integration of Rotation Rates) of the excellent <a href="https://arxiv.org/pdf/1711.02508.pdf" rel="nofollow noreferrer">Quaternion kinematics for the error-state Kalman filter</a> by Joan Solà (2017).</p>
<p>This is only a first-order Taylor series expansion of the time integral of angular velocity. A second-order correction is given in Eq 227 of Solà:</p>
<pre><code>double avgAVX = 0.5 * (avPrev.x + avCurr.x);
double avgAVY = 0.5 * (avPrev.y + avCurr.y);
double avgAVZ = 0.5 * (avPrev.z + avCurr.z);
Quaterniond integral = IMUIntegration.integrateAngularVelocity(
        avgAVX, avgAVY, avgAVZ, dt);
Vector3d correction = avPrev.cross(avCurr, new Vector3d()).mul(dt * dt / 24.0);
Quaterniond qCorrection = new Quaterniond(
        correction.x, correction.y, correction.z, 0.0);
Quaterniond deltaOrientation = integral.add(qCorrection).normalize();
</code></pre>
<p>However the quaternion addition in the final step has no meaningful geometric interpretation, and increases the numerical error over the error introduced by integration at discrete time steps. Therefore, it is better to use a multiplication-only method, as introduced as the predictor–corrector
direct multiplication (PCDM) method in Zhao and Wachem.</p>
<p>This method was improved upon in <a href="https://link.springer.com/content/pdf/10.1007/s00707-016-1670-x.pdf" rel="nofollow noreferrer">Improved quaternion-based integration scheme for rigid body motion</a> by L. J. H. Seelen, J. T. Padding, and J. A. M. Kuipers (2016). If I'm reading the math right, and if you simplify by assuming there are no external torques or forces, and if you throw away the linear velocity components to focus only on angular velocity and rotation, and if you work only in local coordinates rather than switching back and forth between local and world coordinates (which the authors do in order to add in the contributions of external torques, velocities and forces), then you end up with a method that simply integrates the angular velocities between the previous and current IMU readings, and then between the current and next IMU readings, sampling at any desired resolution along these two interpolations, and applying the <code>integrateAngularVelocity</code> routine shown above with the smaller <code>dt</code> value induced by the integration sampling interval:</p>
<pre><code>private static final int NUM_INTEGRATION_STEPS = 6; // Should be even

public static List&lt;IMUIntegrated&gt; integrateAngularVelocities(List&lt;IMURaw&gt; raw,
            double dt) {
    List&lt;IMUIntegrated&gt; integrated = new ArrayList&lt;&gt;();
    for (int i = 0; i &lt; raw.size(); i++) {
        // Find average between curr and prev/next angular velocities
        IMURaw prevRawAV = i == 0 ? null : raw.get(i - 1);
        IMURaw currRawAV = raw.get(i);
        IMURaw nextRawAV = i == raw.size() - 1 ? null 
                : raw.get(i + 1);
        Vector3d prevAV = prevRawAV == null ? new Vector3d()
                : prevRawAV.getAngularVelocity();
        Vector3d currAV = currRawAV.getAngularVelocity();
        Vector3d nextAV = nextRawAV == null ? new Vector3d()
                : nextRawAV.getAngularVelocity();
        Vector3d prevAvgAV = prevAV.add(currAV, new Vector3d()).mul(0.5);
        Vector3d nextAvgAV = currAV.add(nextAV, new Vector3d()).mul(0.5);

        // Find integration interval
        double integrationDt = dt / NUM_INTEGRATION_STEPS;

        // Linearly interpolate and integrate angular velocities
        Quaterniond deltaOrientation = new Quaterniond();
        for (int j = 0; j &lt;= NUM_INTEGRATION_STEPS; j++) {
            double frac;
            Vector3d startAV, endAV;
            if (j &lt; NUM_INTEGRATION_STEPS / 2) {
                frac = (j / (double) (NUM_INTEGRATION_STEPS / 2));
                startAV = prevAvgAV;
                endAV = currAV;
            } else {
                frac = ((j - NUM_INTEGRATION_STEPS / 2)
                        / (double) (NUM_INTEGRATION_STEPS / 2));
                startAV = currAV;
                endAV = nextAvgAV;
            }
            // Linearly interpolate angular velocities
            Vector3d interpolatedAV = startAV.mul(1.0 - frac, new Vector3d())
                    .add(endAV.mul(frac, new Vector3d()));
            // Integrate to produce a quaternion
            deltaOrientation = integrateAngularVelocity(
                    interpolatedAV.x, interpolatedAV.y, interpolatedAV.z,
                    integrationDt)
                // Concatenate onto cumulative transformation
                .mul(deltaOrientation);
        }
        integrated.add(new IMUIntegrated(currRawAV.timestamp,
                deltaOrientation.normalize()));
    }
    return integrated;
}
</code></pre>
Efficient quaternion angular velocity
<p>I have an orientation expressed with a quaternion and an angular velocity expressed as either a quaternion or a number (radians per second around the original orientation). I understand how to do this using conversion to axis-angle but that method is rather computationally expensive and is not a realistic option. How would I go about modifying the orientation quaternion given a time interval (in seconds)? I need a solution for both cases (the quaternion and the number). However, converting one case into the other is acceptable and may be preferable depending on the computational complexity of the various algorithms/formulae required for conversions.</p>