To update the state vector, $y$, we need to know the acceleration, $\ddot{x}_i$, and the angular acceleration, $\dot\omega_i$, for each object, $i$.
These values can be determined through classical mechanics, which will be reviewed below.
Force $F$, measured in Newtons (N), acts translationally on the centre of mass.
Torque $\tau$, measured in Newton-meters (Nm), acts rotationally about the centre of mass:
In 2D, if force $F$ is applied at position $r$ from the centre of mass, it acts for the purposes of translational acceleration as though it were applied at the centre of mass.
But such a force also induces a torque $$\tau = |r|\ F_\perp$$
where $F_\perp$ is the component of $F$ perpendicular to $r$:
In 3D, the torque $\tau$ is a vector about which the rotational acceleration occurs and $$\tau = r \times F$$
The cross product automatically takes only the perpendicular component of $F$ into account, since $$r \times F = r \times (F_\perp + F_\parallel) = r \times F_\perp + r \times F_\parallel = r \times F_\perp$$
since $r$ and $F_\parallel$ are parallel and thus have a cross product of zero.
The classical Newton's laws are $$F = m\ \ddot{x}$$
where $m$ is the mass. And $$\tau = I\ \dot\omega + \omega \times (I \omega)$$
where $I$ is a $3 \times 3$ inertia tensor which describes the distribution of mass in the object.
In what follows, we are (somewhat confusingly, given the notation of $x$ as a vector above) using $x,y,z$ to denote 3D coordinates.
An inertia tensor has the form $$I = \left[\begin{array}{lll} I_{xx} & I_{xy} & I_{xz} \\ I_{yx} & I_{yy} & I_{yz} \\ I_{zx} & I_{zy} & I_{zz} \\ \end{array}\right]$$
with moments of inertia $$I_{xx} = \int_M (y^2 + z^2)\ \rho(x,y,z)\ dx\ dy\ dz$$ $$I_{yy} = \int_M (x^2 + z^2)\ \rho(x,y,z)\ dx\ dy\ dz$$ $$I_{zz} = \int_M (x^2 + y^2)\ \rho(x,y,z)\ dx\ dy\ dz$$
and products of inertia $$I_{xy} = I_{yx} = \int_M xy\ \rho(x,y,z)\ dx\ dy\ dz$$ $$I_{xz} = I_{zx} = \int_M xz\ \rho(x,y,z)\ dx\ dy\ dz$$ $$I_{yz} = I_{zy} = \int_M yz\ \rho(x,y,z)\ dx\ dy\ dz$$
where $M$ is the volume of the object and $\rho(x,y,z)$ is the density at location $(x,y,z)$. The location is taken with respect to the centre of mass.
For a cuboid (i.e. an axis-aligned rectangular solid) of dimensions $a \times b \times c$: $$I_{xx} = {1 \over 12}\ m\ (b^2 + c^2) \ \ \ \ \ \ \ \ \textrm{(with } I_{yy} \textrm{ and } I_{zz} \textrm{ similar)}$$ $$I_{xy} = I_{xz} = I_{yz} = 0$$
For a sphere of radius $r$: $$I_{xx} = I_{yy} = I_{zz} = {2 \over 5}\ m\ r^2$$ $$I_{xy} = I_{xz} = I_{yz} = 0$$
Inertia tensors for other parameterized objects can be found online. Or you can compute them by integrating over the object.
For non-deformable objects, the inertia tensor is constant.
For translations with $F = m\ \ddot{x}$, if $F = 0$ then $\ddot{x} = 0$.
But this is not the case for rotations with $\tau = I\ \dot\omega + \omega \times (I \omega)$.
If the object is not rotating, $\omega = 0$ and $\tau = I\ \dot\omega$. In this case, if $\tau = 0$ then $\dot\omega = 0$.
What if $\omega \neq 0$ (i.e. the object is rotating)? Is it still possible have $\omega \times (I \omega) = 0$ and so have $\tau = I\ \dot\omega$?
Yes! $\omega \times (I \omega) = 0$ when $\omega$ and $I \omega$ are parallel. That occurs when $\omega$ is an eigenvector of $I$, since $I \omega = \lambda \omega$ (for some $\lambda$) in this case.
The eigenvectors of the inertia tensor, $I$, are the object's axes of rotational symmetry. For a cuboid, these directions are the axis directions. For a sphere, these are all directions, since the inertia tensor is a scalar multiple of the identity matrix.
The above shows that, if an object is rotating about an axis that is not an axis of rotational symmetry, the object will undergo angular acceleration even if no external torque is applied.
We can determine $\ddot{x}$ and $\dot\omega$ for an object: $$\textstyle \ddot{x} = {F \over m}$$ $$\dot\omega = I^{-1} (\tau - \omega \times (I \omega))$$
Thus, given the total force and total torque on each object, $i$, we can determine $\ddot{x}_i$ and $\dot\omega_i$.
From there, we can fill in the elements of the derivative, $\dot{y}$, of the state vector, $y$.
And from there, we can integrate to compute the change in state, thus animating our objects.