[[ cos(omega * time), -sin(omega * time), 0] [ sin(omega * time), cos(omega * time), 0] [ 0, 0, 1]]
The final first derivative update matrix ends up looking like this:
[[ cos(omega * time), -sin(omega * time), vx * time] [ sin(omega * time), cos(omega * time), vy * time] [ 0, 0, 1]]
With vx and vy being the linear velocity terms expressed in local space.
Then you just do your position/orientation matrix * first derivative update matrix. However, the problem here is that the first derivative matrix has to be pre-integrated, because for angular velocity the time term is inside the trig functions, and for linear velocity the time term is just a scalar multiply.
Which isn't a huge problem if the time step is fixed. However, I don't see how to add the final term: + 1/2 * alpha * t^2. That is, the second derivative matrix. And I'm not sure how to build such a matrix from forces being applied to produce torque and acceleration.