3D Rigidbody Simulation

Started by
15 comments, last by jalex 4 years, 8 months ago

I'm happy to help. It is always recommended to use a very minimalist simulation configuration to test things out. ;)

Advertisement
 
 
 
0
 Advanced issues found
 
2
On 9/29/2018 at 2:20 PM, Irlan Robson said:

The time derivative of an orientation given the angular velocity of the rotation frame represented by the orientation is

q_dot = 0.5 * q_omega * q

Then to integrate the current orientation q over a time step dt you use the formula

q = q + dt * q_dot

This integration scheme is very common for quaternions, and quite wrong as the resulting quaternion `q_next`  after a finite time step `dt` isn't a rotation anymore (with error O(dt))

Try the following: Convert \( \vec{\omega} \) into a rotation, but converting into an angle \( \theta = dt * \| \vec{\omega} \| \) and axis \( \hat{z} = \vec{\omega} / \| \vec{\omega} \| \). Then define a new quaternion `q_rot` using the axis/angle definition and perform the rotation step $$ q_{\rm next} \Rightarrow q_{\rm rot} q $$

In term of the vector+scalar components of ther quaternion the above is:

$$ \begin{pmatrix}\vec q_{\rm next}\\ q_{\next} \end{pmatrix} \Rightarrow\begin{pmatrix}\left(\tfrac{1}{\omega}\sin\frac{dt\,\omega}{2}\right)\vec\omega\\
\cos\frac{dt\,\omega}{2}
\end{pmatrix}\,\begin{pmatrix}\vec q\\ q \end{pmatrix} $$

On 9/29/2018 at 2:20 PM, Irlan Robson said:
 
 
 
0
 Advanced issues found
 
2
On 9/29/2018 at 2:20 PM, Irlan Robson said:

The time derivative of an orientation given the angular velocity of the rotation frame represented by the orientation is

q_dot = 0.5 * q_omega * q

Then to integrate the current orientation q over a time step dt you use the formula

q = q + dt * q_dot

This integration scheme is very common for quaternions, and quite wrong as the resulting quaternion `q_next`  after a finite time step `dt` isn't a rotation anymore (with error O(dt))

Try the following: Convert \( \vec{\omega} \) into a rotation, but converting into an angle \( \theta = dt * \| \vec{\omega} \| \) and axis \( \hat{z} = \vec{\omega} / \| \vec{\omega} \| \). Then define a new quaternion `q_rot` using the axis/angle definition and perform the rotation step $$ q_{\rm next} \Rightarrow q_{\rm rot} q $$

In term of the vector+scalar components of ther quaternion the above is:

$$ \begin{pmatrix}\vec q_{\rm next}\\ q_{\rm next} \end{pmatrix} \Rightarrow\begin{pmatrix}\left(\tfrac{1}{\omega}\sin\frac{\theta}{2}\right)\vec{z}\\
\cos\frac{\theta}{2}
\end{pmatrix}\,\begin{pmatrix}\vec q\\ q \end{pmatrix} $$

The method I showed is just Euler's method applied to quaternions and a derivation can be found here:

http://web.cs.iastate.edu/~cs577/handouts/quaternion.pdf

I've been using this method for years without a single problem. Of course I normalize the quaternion after integration. This is fast.

What you showed seems more like an alternative to Euler's method. And you said it is better than Euler's method. So do you have a proof for this? A paper maybe?

1 minute ago, Irlan Robson said:

The method I showed is just Euler's method applied to quaternions and a derivation can be found here:

http://web.cs.iastate.edu/~cs577/handouts/quaternion.pdf

I've been using this method for years without a single problem. Of course I normalize the quaternion after integration. This is fast.

What you showed seems more like an alternative to Euler's method. And you said it is better than Euler's method. So do you have a proof for this? A paper maybe?

I understand what you are saying, but if you do the math you will see that the angle drifts away with order O(ω*dt) which is the same error order as an Euler step. So if that is ok with the simulation and any long term numerical artifacts, that is ok. But for some cases a higher order integrator is needed (like Runge Kutta) and then this method fails.

In the attached paper there is an example in the end. Starting with a quaternion from the rotation about the \( x \)-axis, apply some rotational velocity about another axis, like \( z \). Do that with many small Euler steps, and you will get the same answer as the method I propose with 1 step.

quaternion_step.pdf

What Jalex suggests is actually what most physics engine use. This is often referred to as 'exponential map'. The funny thing is that we had the same discussion like 15 years ago on the ODE or Bullet forums and now it comes up again as something new :)

The most referenced paper on this (which actually covers a real practical implementation and also deals with numerical issues for very small angular velocities) is this paper by Grassia from 98:

https://www.cs.cmu.edu/~spiff/moedit99/expmap.pdf

14 hours ago, Dirk Gregorius said:

What Jalex suggests is actually what most physics engine use. This is often referred to as 'exponential map'. The funny thing is that we had the same discussion like 15 years ago on the ODE or Bullet forums and now it comes up again as something new :)

The most referenced paper on this (which actually covers a real practical implementation and also deals with numerical issues for very small angular velocities) is this paper by Grassia from 98:

https://www.cs.cmu.edu/~spiff/moedit99/expmap.pdf

Great reference! Thank you.

This topic is closed to new replies.

Advertisement