Games Series #1 Simulating a Rope

I believe that in order to better use an object in a video game, one needs to have some general knowledge about it in the real world. This helps with design decisions—where, when, how do we use this in the game? What mechanics can it enable? Understanding the real world properties of the object and of the materials that the object may be made of might also influence various implementation decisions—how should this object look? How should it interact with other objects in the game?

Consequently, I believe a bit of background about ropes is in order. To make a rope, we start with some fibers. Those may be natural (hemp, linen, cotton, etc) or syntehtic (nylon, polypropylene, etc). We twist those fibers together into what we call yarns, and those yarns into strands. Finally, we twist those together to make rope.

The origin of rope dates back to prehistoric times, and through the years it played an important role in the development of humankind. Rope can be made by hand, but since about 4000 BC, people started using machines to make it. It is believed that the Egyptians were the first to invent such devices, probably because rope was indispensable to their agriculture and the building of their structures, including their temples and pyramids. They usually made rope out of materials such as water reed, date palms, flax, grass, papyrus, leather, or even camel hair.

One well-known method of making rope employs the concept of a ropewalk—a technique used mostly by Europeans in the 17th and 18th centuries. You can see how it works in Science Channel’s short video from below, from their How It’s Made series.

The Euler method is one of the simplest and most popular numerical integrators. In essence, it is a way of approximating solutions to a first order ordinary differential equation with a given initial value, but it can approximate solutions to ordinary equations of any order, since all of them can be represented as systems of multiple first order equations. The equations of motion are just one example of this kind.

To understand how it works, first we need to understand the idea of approximating with the tangent. For simplicity, say we have a function \( f: \mathbb{R}\rightarrow \mathbb{R} \) that is differentiable at some point \( t \in \mathbb{R}. \) We would like to know the value of f at point \( t + \Delta t, \) but for some reason this information not available to us. To approximate it, we can use the tangent to the graph of f at t. Imagine that instead of moving along the graph of f from t towards \( t + \Delta t\), we move along the tangent. If we don't go very far, we should still be close to the graph of f, i.e. the smaller \( \Delta t, \) the smaller the error in the approximation. Written mathematically, this means that \( f(t) + \Delta t \, f'(t) \) can be used as an approximation of \( f(t + \Delta t). \)

This concept is illustrated in the example below. You can drag the blue tangent along the graph of a function f. Keep an eye on the red error, which you can drag to adjust \( \Delta t. \)

Suppose now that we don't even know the value of f at t; suppose we don't know the value of f at any other point. Instead, suppose we know that if the value of f equals y for some input x and some output y, then we know what \( f'(x) \) is. Furthermore, suppose we have this information for all \( x, y \in \mathbb{R}. \) We would like to find f by using this information. This is called a first order ordinary differential equation—an equation involving x, f, and \( f'. \) Mathematically, this can be written as \[ F(x, f, f') = 0, \] where f is the unknown function. This kind of equation usually comes with some initial condition for f, which is needed in order to find an exact or approximate solution. This initial condition states that \( f(x_0) = y_0, \) for some initial choice of the point \( P_0 = (x_0, y_0). \)

The Euler method begins at this point \( P_0. \) We know that \( f(x_0) = y_0. \) As described above, if we know \( f(x_0) \) and \( f'(x_0), \) we can approximate \( f(x_0 + \Delta t) \) to be close to \( f(x_0) + \Delta t \, f'(x_0). \) This gives us a new point \[ \begin{align} P_1 &= \left(x_0 + \Delta t, y_0 + \Delta t \, f'(x_0)\right)\\ &= (x_1, y_1). \end{align} \] Starting with \( P_1, \) the point \( P_2 \) is found in a similar manner, by using the same approximation technique. So are \( P_3, \) \( P_4 , \) etc. These discrete points put toghether make an approximate solution to the equation.

A good way to visualise first order ordinary differential equations is by the means of a slope field. You can see one in the example below—the small lines in the background. One such line represents the slope of a potential solution f, should its graph pass through that point. In the case of the Euler method, you could think of them as indicators of where to go, given your current position. For example: if you find yourself at, say, the point (3, 2) and you know that the slope there is 1.5, you might choose to go one unit to the right, in which case you would also have to go one and a half units above, ending up at (4, 3.5). You might as well choose to go only 0.5 units to the right, in which case you would have to go \( 0.5 \times 1.5 = 0.75 \) units above, ending up at (3.5, 2.75). The number of units you choose to go to the right is called step size. In our equations, we called it \( \Delta t. \) If you start at the point given by the initial condition of the equation and take multiple such steps, you end up with what we call a discrete approximate solution.

In this specific example below, the differential equation was chosen to be \( f'(x) = \sin(x) + f(x). \) The exact analytic solution is a function of the form \( f(x) = c e^x - \sin(x) / 2 - \cos(x) / 2, \) where c depends on the initial condition. The blue curve represents one such solution for a chosen initial condition, which may be changed by moving the blue handle. Each red circle represents a point found by the Euler method. You can drag any of them to change the size of the step, \( \Delta t. \)

In the Euler method we use the approximation with the tanget. How can we improve on that? Well, suppose we know the value of some function g at \( 0 \) and that it’s also infinitely differentiable at \( 0. \) We can say that \[ g(\Delta t) = g(0) + \int_{0}^{\Delta t} g'(t) \, dt. \] Note that since t goes from \( 0 \) to \( \Delta t \) inside the integral, we can replace it with \( \Delta t - t \) without changing the end result: \[ g(\Delta t) = g(0) + \int_{0}^{\Delta t} g'(\Delta t - t) \, d t. \] \[ g(\Delta t) = g(0) + \\ + \int_{0}^{\Delta t} g'(\Delta t - t) \, d t. \] Integrate by parts to get \[ g(\Delta t) = g(0) + \Delta t \, g'(0) - \int_{0}^{\Delta t} t \, g''(\Delta t - t) \, d \Delta t. \] \[ g(\Delta t) = g(0) + \Delta t \, g'(0) - \\ \int_{0}^{\Delta t} t \, g''(\Delta t - t) \, d \Delta t. \] If we keep integrating by parts and then take \( f(t + \Delta t) = g(\Delta t) , \) we obtain what we call the Taylor expansion of f at t: \[ f(t + \Delta t) = f(t) + \Delta t \frac{f'(t)}{1!} + \Delta t^2 \frac{f''(t)}{2!} + \Delta t^3 \frac{f^{(3)}(t)}{3!} + \dots \] \[ f(t + \Delta t) = \sum_{n=0}^{\infty} \frac{f^{(n)}(t)}{n!} \Delta t^n. \] Note that if we only keep terms linear in \( \Delta t, \) we obtain the formula used in the Euler method. We can improve the accuracy by adding more and more higher order terms. However, if we want to approximate the position, we can’t go higher than quadratic, since that would require us to know \( s^{(3)}(t) = a'(t). \) We could choose to go only with terms that are quadratic or lower: \[ s(t + \Delta t) \approx s(t) + \Delta t \, v(t) + \frac{1}{2} \Delta t^2 \, a(t), \] \[ s(t + \Delta t) \approx s(t) + \Delta t \, v(t) + \\ + \frac{1}{2} \Delta t^2 \, a(t), \] and the local error would be \( \mathcal{O}(\Delta t^3). \) Unfortunately, we would not be able to use the same approach for \( v \) because, again, it would require us to know \( a'(t). \) The good news is that we can employ a mathematical trick to get rid of \( v \) and all other odd order derivatives entirely. If we expand \( f(t - \Delta t) \): \[ f(t - \Delta t) = f(t) - \Delta t \frac{f'(t)}{1!} + \Delta t^2 \frac{f''(t)}{2!} - \Delta t^3 \frac{f^{(3)}(t)}{3!} + \dots \] \[ f(t - \Delta t) = \sum_{n=0}^{\infty} \frac{f^{(n)}(t)}{n!} (-\Delta t)^n \] and then add both expansions together, we get \[ f(t + \Delta t) + f(t - \Delta t) = 2 \left( f(t) + \frac{\Delta t^2 f''(t)}{2!} + \frac{\Delta t^4 f^{(4)}}{4!} + \dots \right). \] \[ f(t + \Delta t) + f(t - \Delta t) = \\ =2 \sum_{n=0}^{\infty} \frac{f^{(2n)}(t)}{(2n)!}\Delta t^{2n}. \] By subtracting \( f(t - \Delta t) \) from both sides and giving up all the terms that are higher than quadratic, we obtain the final approximation formula \[ f(t + \Delta t) \approx 2f(t) - f(t - \Delta t) + \Delta t^2 \, f''(t). \] \[ f(t + \Delta t) \approx \\ 2f(t) - f(t - \Delta t) + \Delta t^2 \, f''(t). \] This is the formula used in Verlet integration. Similarly to the Euler method, the Verlet method says that we can repeatedly apply it to integrate, in our case, the equations of motion. Note that we only need to do it for the position—the velocity is implicit.

To get a better understanding, let’s take a closer look at it. The first thing to note is that there is no first derivative term, we only use the second derivative. If \( f = s, \) the position of some particle, note how the right-hand side can be rewritten as \[ s(t) + \Delta t \frac{s(t) - s(t - \Delta t)}{\Delta t} + \Delta t^2 \, a(t). \] \[ s(t) + \Delta t \frac{s(t) - s(t - \Delta t)}{\Delta t} + \\ + \Delta t^2 \, a(t). \] But \( s(t) - s(t - \Delta t) \) divided by \( \Delta t \) is just the average velocity of the particle between times \( t - \Delta t \) and \(t\). In other words, to approximate \( s(t + \Delta t), \) with Euler we had to add the current instantaneous velocity times \( \Delta t \) to \( s(t), \) but with Verlet we have to add the average velocity between the previous interval, from \( t - \Delta t \) to t, times \( \delta t \) and the instantaneous acceleration acceleration times \( \Delta t^2 \).

It turns out that making these two simple changes improves the global error from \( \mathcal{O}(\Delta t) \) to \( \mathcal{O}(\Delta t^2). \) Not perfect, but this together with ease of implementation, the implicit velocity term, and low computational cost, makes Verlet an especially good choice for an integrator in a soft body simulation.

In pseudocode, a procedure using the Verlet integration method to update the positions of a set of particles might look as follows:

Thomas Jakobsen devised his method as the head of research at IO Interactive, while working on Hitman: Codename 47, a game which I replayed many times.

The game featured realistic ragdolls that the player could interact with. A first glimpse of this new system was present in the tutorial level, where the player was taught how to use a couple of items on an actual ragdoll, setting the expectations for what was to come.

The paper was published in the format of an HTML page, and was presented at Game Developers Conference 2001. It has since been cited many times, becoming somewhat of a de facto resource when diving into the world of soft bodies simulation.

The paper doesn't argue why the method works, formalizations were published later (more about those in the next section). In this section we'll take a closer look at the Jakobsen method in the hopes of obtaining some insights at an intuitive level.

As it can be seen in the table, after \(n\) passes, the second constraint is always satisfied. Only the first constraint remains unsatisfied, with an error equal to \( \left( \Delta d_1 + 2\Delta d_2 \right) / 2^{2n} \) that goes to \(0\) as \(n\) goes to infinity.

In this particularly small case the algorithm works. Below you can find an animation of how it works when there are 10 particles in a one dimensional system. Segments are represented as red rectangles when their length is far from the target length d. They turn blue as their length gets closer to d. Note that this process, theoretically, doesn't really ever end, but at some point it gets so close to satisfying all constraints completely that you won't be able to notice its progress anymore.

The idea of pass can be observed clearly in this example. You can see how a completely relaxed constraint in complete blue goes from left to right as the algorithm completes the pass. You can see how, with each pass, the system moves closer to satisfying all of its constraints. Error is eliminated from the system at each first and last steps of the pass.

Of course, a trivial example with only one dimension may not do much to help us understand how it works at an intuitive level. Things get more interesting in two-dimensional environments. Let us visualize one more example—this time in 2D.

The classic way of simulating dynamic systems is force based. What this means is that in order to find the positions of some particles in a system, one looks at the forces acting on those particles and infers their positions by means of numerical integration. Because it can offer very high accuracy, this approach is usually used in fields where accuracy is important, such as in computational chemistry or computational physics. In computer games, however, accuracy is not necessary. Believability and simplicity are more desirable features.

In 2006, Matthias Müller, Bruno Heidelberger, Marcus Hennix, and John Ratcliff published a paper titled Position Based Dynamics. They argued that in the case of computer games and other systems that do not require a high degree of accuracy, we may skip the forces and instead work directly on positions. This is exactly what we did with the Jakobsen method when we enforced constraints by directly manipulating the positions of the particles. Jakobsen, along with Müller, actually ended up being considered a founder of the field.

With this approach, we forego some accuracy in exchange for speed and controllability. Both of these are useful in the case of video games, since games need to be both fast and interactive.

At each step we want to simulate the physics of the particles in the system based on their positions and velocities while keeping the constraints satisfied. The first part may be achieved by the use of a numerical integrator, such as Euler or Verlet. The problem is that this simulation substep may break the constraints. Thus, before the final result is returned, we need to add in another substep in which the positions are adjusted such as to keep the constraints satisfied. Usually, in practice, we won't be able to satisfy all of them perfectly, so we resort to an approximation.

Mathematically, to fix the constraints, we need to solve (i.e. find an approximate solution to) the system of non-linear equations \[ f_c(x) \succ_c 0, \forall c \in C. \]

Even though the system we are dealing with is non-linear, let us for a moment take a step back and present a set of methods that deal with system of linear equations. Then we are going to see how we can generalize them for non-linear systems.

Suppose we have a square linear system of equations \( Ax = b. \) Suppose, however, that the matrix A is too complicated to work with directly. What we can do instead is replace it with a new matrix M and move the difference \( A - M \) to the other side of the equation: \[ Mx = (A - M)x + b. \] Now, if we need an exact solution, we are still forced to deal with A, nothing changes. On the other hand, if an approximation suffices, we may solve this new equation iteratively as follows. We begin with an initial guess \( x^{(0)} \) and at each step we take the current \( x^{(k)} \) and find an even better approximation \( x^{(k + 1)} \) by solving \[ Mx^{(k)} = (A - M)x^{(k + 1)} + b. \] This equation may be simpler, but it never gives the exact solution—we need to keep solving it over and over. In practice, we keep solving it until the error is below some desired threshold, or as long as the time budget allows.

On top of that, not every M will work. Intuitively, the closer it is to A, the faster the iteration converges to the exact solution. \( M = A \) converges in just one interation. On the other hand, M also needs to be as simple as possible. There are, in fact, three choices of M that are both simple and with good convergence rates when they do converge:

• \( M_1 = \text{the diagonal part of } A \) (the Jacobi method),
• \( M_2 = \text{a triangular part of } A \) (the Gauss-Seidel method),
• \( M_1 + M_2 \) (successive overrelaxation).

Keep in mind that the convergence of \( x^{(k)} \) when \( k \rightarrow \infty \) is not guaranteed. However, if \( x^{(k)} \) does converge to some \( x^{(\infty)} \) as \( k \rightarrow \infty, \) it is true that \( A x^{(\infty)} = b, \) i.e. \( x^{(\infty)}, \) if it exists, is the correct solution to our initial system of equations.

Let us direct our attention to the choices \( M_1 \) and \( M_2 \) for the matrix M, and how the Jacobi and the Gauss-Seidel methods act on individual components of \( x^{(k)}. \) In the case of Jacobi, we have \[ x_i^{(k + 1)} = \frac{1}{a_{ii}} \left( b_i - \sum_{j \neq i } a_{ij} x_j^{(k)} \right). \] Note here that to find \( x_i^{(k + 1)}, \) Jacobi only uses the components of \( x^{(k)}. \) Thus, the components of \( x^{(k + 1)} \) can be computed independently and in parallel. On the other hand, with Gauss-Seidel we have \[ x_i^{(k + 1)} = \frac{1}{a_{ii}} \left( b_i - \sum_{j = 1 }^{i - 1} a_{ij} x_j^{(k + 1)} - \sum_{j = i + 1 }^{n} a_{ij} x_j^{(k + 1)} \right). \] \[ x_i^{(k + 1)} = \frac{1}{a_{ii}} \left( b_i - \sum_{j = 1 }^{i - 1} a_{ij} x_j^{(k + 1)} - \\ \sum_{j = i + 1 }^{n} a_{ij} x_j^{(k + 1)} \right). \] Let us take a moment and understand how this formula is different from Jacobi’s. They are quite similar, but there is one key difference. With Gauss-Seidel, after a component \( x_i^{(k + 1)} \) is computed, it is then used in the computation of all subsequent components \( x_j^{(k + 1)}, j > i. \) We make use of new approximations as they become available. This is similar to what we did with the Jakobsen method.

As you might have noticed, all of these methods are for linear systems, but we are dealing with a non-linear system. The good news is that we can employ a mathematical trick to overcome this impediment.

As a reminder, the system we are trying to solve is \[ f_c(x) \succ_c 0, \forall c \in C, \] where x represents the concatenation of the position vectors of the particles in the system. Instead of solving the entire system at once, we are going to solve each equation \( f_c(x) \succ_c 0, \) for a specific c, individually. In similar Gauss-Seidel fashion, after we find a solution x given by the condition of some constraint c, we use it as a starting point to find a solution to the equation given by the next constraint. In this way, we keep finding values for x that, as we'll see, converge to a state where all the constraints are perfectly satisfied.

After x is updated to satisfy the condition of some constraint c, we say that the constraints c is relaxed. With this in mind, we say that our algorithm works by successively relaxing all the constraints in the system.

But how do we relax a constraint? Each one yields one scalar equation \( f_c(x) \succ_c 0 \) with up to \( 3n \) unknowns. This is an underdetermined non-linear system. Let us frame the problem as follows: given an initial value for \( x, \) say \( x^{(k)}, \) we want to find a correction \( \Delta x \) such that \( f_c(x^{(k)} + \Delta x) \succ_c 0 . \) Then we say \( x^{(k + 1)} \) is \( x^{(k)} + \Delta x. \) Now we need to find \( \Delta x \) instead, but it’s still a non-linear system. Fortunately we may employ the concept of linearization around a point, in this case \( x^{(k)} \): \[ f_c(x^{(k)} + \Delta x) \approx f_c(x_k) + \nabla f_c(x_k) \cdot \Delta x \succ_c 0, \] \[ f_c(x^{(k)} + \Delta x) \approx f_c(x_k) + \\ + \nabla f_c(x_k) \cdot \Delta x \succ_c 0, \] where \( \nabla f_c(x^{(k)}) \) is the gradient of \( f_c \) at \( x^{(k)}, \) a vector made up of all the partial derivatives with respect to each component of \( x^{(k)}. \) The system is now linear, but still underdetermined. This is mended by restricting \( \Delta x \) to be in the direction of \( \nabla f_c(x^{(k)}). \) The problem is now reduced to finding one scalar \( \lambda \) such that \[ f_c(x^{(k)}) + \lambda \lVert \nabla f_c(x^{(k)}) \rVert^2 \succ_c 0. \] This results in \[ \lambda = -\frac{f_c(x^{(k)})}{ \lVert \nabla f_c(x^{(k)}) \rVert^2 } \] and thus \[ \Delta x = -\frac{f_c(x^{(k)})}{ \lVert \nabla f_c(x^{(k)}) \rVert^2 } \nabla f_c(x^{(k)}). \] To get a better understanding, let us zoom in and inspect what happens in the case of just one particle \( p \in P. \) First, note that the numerator, which is the squared length of the gradient, can be rewritten as the sum of the squared lengths of the subvectors of the gradient: \[ \lVert \nabla f_c(x^{(k)}) \rVert^2 = \sum_{p \in P} \Big\lVert \left[ \nabla f_c(x^{(k)}) \right]_p \Big\rVert^2, \] where \( \left[ \nabla f_c(x^{(k)}) \right]_p \) is the subvector of the gradient that contains only the partial derivatives with respect to components of \( x_p. \) Consequently, the correction for a single particle p is \[ \Delta x_p = -\frac{f_c(x^{(k)})}{ \sum_{r \in P} \lVert \left[ \nabla f_c(x^{(k)}) \right]_r \rVert^2 } \left[ \nabla f_c(x^{(k)}) \right]_p. \] \[ \Delta x_p = -\frac{f_c(x^{(k)})}{ \sum_{r \in P} \lVert \left[ \nabla f_c(x^{(k)}) \right]_r \rVert^2 } \\ \left[ \nabla f_c(x^{(k)}) \right]_p. \] This formula is something we could use in practice. In fact, we did use it when we implemented the Jakobsen method. Let us again consider the case of two particles and a single constraint with the function \[ f(x = [ x_1; x_2 ]) = \lVert x_1 - x_2 \rVert - d \] which is satisfied when \( f(x) = 0. \) We have that \[ \nabla f(x) = \frac{[ x_1 - x_2; x_2 - x_1 ]}{\lVert x_1 - x_2\rVert} \] and thus \[ \begin{align} \big[ \nabla f(x) \big]_1 &= \frac{x_1 - x_2}{\lVert x_1 - x_2\rVert}, \\ \big[ \nabla f(x) \big]_2 &= \frac{x_2 - x_1}{\lVert x_1 - x_2\rVert}. \end{align} \] For \( \lambda \) we have \[ \lambda = - \frac{ \lVert x_1 - x_2 \rVert - d }{ 2 }. \] Let us say that \( \Delta d = \lVert x_1 - x_2 \rVert - d. \) Finally, we get the individual corrections: \[ \begin{align} \Delta x_1 &= \frac{ \Delta d (x_2 - x_1) }{ 2 \lVert x_1 - x_2 \rVert }, \\ \Delta x_2 &= \frac{ \Delta d (x_1 - x_2) }{ 2 \lVert x_1 - x_2 \rVert }, \end{align} \] which are the exact same corrections we used in our implementation of the Jakobsen method.

There are many more features we could add to our simulation. Masses may be incorporated by using the inverse of a \(3n\) by \(3n\) mass matrix \[ M = \begin{bmatrix} m_1 & 0 & 0 & 0 & \dots & 0 \\ 0 & m_1 & 0 & 0 & \dots & 0 \\ 0 & 0 & m_1 & 0 & \dots & 0 \\ 0 & 0 & 0 & m_2 & \dots & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & 0 \\ 0 & 0 & 0 & 0 & \dots & m_n \end{bmatrix} \] (where \( m_p \) is the mass of the particle p) to transform the gradient in the computation of \( \Delta x \) as such \[ \Delta x = \lambda M^{-1} \nabla f_c(x^{(k)}). \] Stiffness may be introduced in various ways. For example we may multiply the correction by a stiffness constant \( k \in [0 \dots 1] \) when applying the correction as such: \[ x^{(k + 1)} \leftarrow x^{(k)} + k \Delta x. \] This is, unfortunately, dependent on the number of iterations we do, but can be fixed by doing \[ x^{(k + 1)} \leftarrow x^{(k)} + \big(1 - \left(1 - k\right)^s\big) \Delta x, \] instead, where \( s \) is the number of iterations.

We saw how to implement a basic rope simulation. We also briefly examined the mathematical framework that sits behind the implementation we used. Let me end this by highlighting that there are many more features one could add to this implementation. Of course, there are also other ways, including completely different mathematical frameworks, to achieve the same effect. For example, instead of position based dynamics, one might use a spring-mass model—this would be more like the “traditional” way of doing a physics simulation, one where we don't skip the force and velocity layers when enforcing constraints.

A word of warning, though—soft bodies are ferocious creatures. They can quickly devour the CPU and / or GPU if used excessively. Using them in complex three-dimensional worlds where many other things need to be going on at the time needs to be done with care. They may drive down performance significantly, especially on older and less powerful machines.

On the other hand, they can provide that extra something to your game, that extra kick that can make the difference between a dull and an immersive environment. They can make the world you're building seem more lively by bringing a touch of softness in a sea of rigid objects.