# How numerical integration for multibody systems works

I’ve got some great references and quotes here that help me explain how exactly we are able to numerically calculate the dynamics of a multibody system.

## What is a state space? Cornelius Lanczos helps explain

I’m currently reading the superb 1949 book from Cornelius Lanczos: “The Variational Principles of Mechanics”. At page 12 – 13, he provides a great explanation for the so-called configuration space:

“No matter how numerous the particles constituting a given mechanical system may be, or how complicated the relation existing between them, the entire mechanical system is pictured as a single point of a many-dimensional space, called the “configuration space”. For example, the position of a rigid body – with all the infinity of mass points which form it – is symbolized as a single point of a 6-dimensional space. This 6-dimensional space has, to be sure, nothing to do with the physical reality of the rigid body. It is merely correlated to the rigid body in the sense of a one-to-one correspondence. The various positions of the rigid body are “mapped” as various points of the 6-dimensional space, and conversely, the various points of a 6-dimensional space can be physically interpreted as various positions of a rigid body.”

When we include dynamic values – usually velocity – in the coordinates, we move to a so-called phase-space or state-space. There appear to be some subtle differences between them that I can’t really grasp (e.g. whether or not this space is euclidean, whether each coordinate corresponds to one actual degree of freedom or additional boundary conditions are used?). For now, it’s enough to understand that each possible state of our system can be mapped to a certain set of variables.

## How can we calculate our system dynamics? Erin Catto provides a great introduction to numerical solvers

Erin Catto, author of the box2d physics engine, has published an easy-to-follow introduction to numerical solvers (PDF) at the box2d downloads section:

The first 20-something pages give a nice introduction to differential equations and newton’s law for deriving them. On page 25, you’ll find the concept of the phase space. Pages 29-31 show how to turn differential equations into their state-space form – which is needed for most numerical solvers (including the family of ode solvers in Matlab and GNU Octave).

Numerical solutions are introduced from page 31 on. The basic idea is that if you know your system’s current state (you do if you have a set of values for your state space!) and where it’s supposed to be moving next (you do if you have a set of differential equations of your system in state-space form!), then you can incrementally update your system space and move forward (or backward) in time. This is what a numerical solver for differential equations is doing. The explanation includes different solvers, focusing on the relatively simple explicit, implicit and semi-implicit euler algorithms. This is a good start to get an idea on what numerical stability and numerical damping is.

## Where do we go from here? Ideas for self-study and further reading

This blog focuses on multibody dynamics with a science and engineering background. For most applications beyond games, the euler algorithms described above have way too much error – but they are easy to play with and introduce you to all the important concepts of numerical integration. As starting points for self-study, you might want to:

• Implement the different Euler algorithms in Excel (or any other spreadsheet software) for systems with one, two or more degrees of freedom.
• Code them into a programming language of your choice, no matter if you want to use C, Java, Matlab / GNU Octave or anythign else.

Going beyond, you may want to have a look at how the 4th-order Runge-Kutta algorithm works – see for example my implementation for Arduino. In a nutshell, it takes the same approach as explicit euler, but using 4th-order polynomials (well… sort of) instead of constant slopes for extrapolation, making it way more accurate. I haven’t had a look at how other, more complex methods, work, but this Wikipedia article appears to be a good starting point. For solving ODEs and trying different solvers, Matlab is a great testbed and I can especially recommend their documentation here on how to solve differential equations.

## Update: Example in Processing

I’ve just updated my code in processing for interactive double-pendulum simulations to include a parallel simulation of explicit euler, implicit euler and runge-kutta in parallel. You can use the buttons to toggle simulation speed and variants. For me, this seems a great example for comparison of the three possible integration algorithms! Also, I was surprised to see how fast the results diverge (with symplectic euler being way closer to runge-kutta than the horrible explicit euler): (state after ~one pendulum oscillation)

Advertisements