On the impact of explicit or semi-implicit integration methods over the stability of real-time numerical simulations

by   Teodor Cioaca, et al.

Physics-based animation of soft or rigid bodies for real-time applications often suffers from numerical instabilities. We analyse one of the most common sources of unwanted behaviour: the numerical integration strategy. To assess the impact of popular integration methods, we consider a scenario where soft and hard constraints are added to a custom designed deformable linear object. Since the goal for this class of simulation methods is to attain interactive frame-rates, we present the drawbacks of using explicit integration methods over inherently stable, implicit integrators. To help numerical solver designers better understand the impact of an integrator on a certain simulated world, we have conceived a method of benchmarking the efficiency of an integrator with respect to its speed, stability and symplecticity.



There are no comments yet.


page 9


Numerical time integration of lumped parameter systems governed by implicit constitutive relations

Time-integration for lumped parameter systems obeying implicit Bingham-K...

Efficient IMEX Runge-Kutta methods for nonhydrostatic dynamics

We analyze the stability and accuracy (up to third order) of a new famil...

Linearly Implicit Multistep Methods for Time Integration

Time integration methods for solving initial value problems are an impor...

Time-Accurate and highly-Stable Explicit operators for stiff differential equations

Unconditionally stable implicit time-marching methods are powerful in so...

Symbolic Multibody Methods for Real-Time Simulation of Railway Vehicles

In this work, recently developed state-of-the-art symbolic multibody met...

Walking into the complex plane to 'order' better time integrators

Most numerical methods for time integration use real time steps. Complex...

Modern theory of hydraulic fracture modeling with using explicit and implicit schemes

The paper presents novel results, obtained on the basis of the modified ...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Intricate physical and mathematical models for simulating soft or rigid bodies are the centre of intensive research efforts. Computer animation, laparoscopic haptic surgery simulations, graphical special effects or robotic manipulation are just a few of the fields where such models play a key role. The focus of our research is to address the issue of carrying out stable, physics-based simulations at interactive update rates. For this purpose, we have conceived an elementary framework for building deformable objects with soft or stiff constraints. Using this framework, we can test the efficiency and impact of several, well-known numerical integrators. To address our goal, we look find an equilibrium between an integrator’s computational overhead, its precision, symplecticity and, most importantly, its inherent stability.

Following the aforementioned goals, we present a short literature survey on the numerical simulation of soft or semi-rigid bodies in section 2. The mathematical apparatus for creating and simulating constrained objects is explained in section 3, and our own, deformable linear object model case study is presented in section 4. We present a general update strategy that supports any type of explicit or semi-implicit integration method in section 5. Finally, the efficiency of several popular integrators is analysed and discussed in section 6, and we conclude this research in section 7.

2 Related work

Stability as a central attribute to Physics-based object simulation within the Computer Graphics world was addressed by directly consider inherently stable integrators. The most popular method employed is the first order, implicit Euler scheme. Servin et al. [SLM06] treated infinitely stiff springs as kinematic constraints and developed a system capable of simulating elastic behaviour at the cost of solving a large sparse linear system for each iteration. The hybrid method of Schroeder et al. [SKZF11] uses explicit updates for the elastic forces and implicit strategies for other components. This complex idea tries to reduce the effect of suppressing material vibrations due to the use of pure implicit integrators. By alternating between implicit, semi-implicit and explicit methods, Volino et al. [VMT05] demonstrated how different cloth materials can be modeled to overcome the inabilities of one integrator to support a certain material property. Some object models allow using for computing an analytical force approximation for an implicit solver as demonstrated by Mesit et al. [MGC07]. Their method of simulating gas filled soft bodies can thus easily support any integrator, but it relies on the topological structure for all analytical derivations. Finally, structured mass-spring models using pure implicit integrators were popularized by Baraf and Witkin [BW98] in their seminal paper on stable cloth simulation and by Desbrun et al. [DSB99] where inverse dynamics were used to tackle outstretching artefacts. For a discussion of the performance of some popular implicit solvers we direct the user to the paper of Hauth et al. [HE01].

Explicit methods are more popular due to their reduced complexity. Finite element simulations using explicit integrators were performed by Fierz et al. [FSH10]. To stabilize the numerical process, the authors proposed modifying the stiffness matrices of ill-conditioned tetrahedral elements. This allowed their simulation to use higher time steps in spite of the explicit updates. Although less computationally demanding and stable than their implicit counterparts, explicit integrators can share a fair amount of stability, energy conserving capabilities and time reversibility. In this respect, Tsai et al. [TKL04] present a comprehensive survey on symplectic integrators used in molecular dynamics.

In a different class of their own, position-based dynamics (see Bender et al. [BMOT13] for a survey) offer a workaround for any force or impulse based system, avoiding the intricacies of using integrators. This family of methods is, however, inaccurate for scenarios where velocities and forces need to be measured, hence we mention it as an alternative for special effects applications.

For a comprehensive list of simulation methods involving deformable objects, we invite the reader to consult the work of Nealen et al. [NMK06] or that of Jimenez [Jim12].

3 Constrained object animation

To better understand the elements involved in most physics-based animation scenarios, we will briefly introduce a simple constraint enforcing system. Supporting this kind of simulation mechanism requires a discrete geometrical sampling of the initial shape of the object. For exemplification purposes, we consider a deformable linear object, embedded in a 3D space, whose cross-section is considerably smaller than its length. Its elastic properties that determine how its discrete structure changes are implemented by adding geometric constraints enforced by potentials (see the work of Teschner et al. [THMG04] for a more detailed and generalized application). These potentials provide a direct measure of how the structure of points differs from its initial configuration.

Generally, if are the vertices of a constrained group, then a constraint function is defined as follows:


This function incorporates additional information relating the current vertex configuration of the group to the initial geometrical image through specific scalar measures of length, area, angle, volume, etc. For example, a length based constraint incorporates the initial or rest lengths of directly connected vertices as a reference for the measure of deviation. Generally, an energy function produces only positive amounts and can be written as:


For a configuration, we now consider the restriction of the energy function at a node . This function is written as:


Since the gradient vector coincides with the direction of maximum potential increase, it is natural to consider a penalty function that points in the opposite direction:


This last equation can be written in an equivalent form:


The simplest example of such a behaviour is the case of a linear spring connecting two points and . The natural constraint function is defined as:


where denotes the spring’s rest length. From this constraint function, a corresponding elastic deformation potential energy can be expressed as:


We define the following operator:


also known in some works as the variational or functional derivative operator.

An elastic spring force at node can be derived using equation 7:

After applying several derivative computation rules (chain rule and the derivative of a product of two functions), we find the general elastic force expression:


where is the linear spring’s stiffness coefficient.

A deformable object defined as a closed surface can be discretized by dividing its interior volume into tetrahedral cells. Apart from linear springs, it is desirable to enforce local constraints aimed at preserving the volumes under deformation. Such forces are easy to introduce by deriving them from a volume preserving constraint function involving tetrahedral cells:


where is the volume of the tetrahedron in its rest configuration. Subsequently, the volume preserving potential is:


Using the same method as for linear springs, we compute the force at the vertex by using the operator:


After conveniently arranging the results from the derivation of equation (12), we can write down the expanded formula of this force:


where is an added-in stiffness coefficient.

Figure 1: Volumetric force on a tetrahedral cell.

Examining figure 1, we can see how the volumetric force acts to prevent volume changes. For example, if point is shifted and the volume increases as a consequence, the direction of the volumetric force is given by the cross product vector . Inconsistent or ”flipped” cell configurations, as well as degenerate tetrahedra can and are likely to be encountered during a simulation. Compressible bodies prevent themselves from being completely flattened by acting as non-linear spring elements. When their state is close to a collapse, the elastic forces should increase asymptotically towards infinity. Since the volumetric force component acts like a vertex to face spring in our system, we add a non-linear spring component, as suggested in [CM97]. Hence, if is the current volume of a tetrahedron, then we can write the expression of the improved volumetric force as:


where is the sign function defined on the set of real values.

Another scenario that can prevent the simulation from recovering its initial rest shape in the absence of perturbing forces is the inversion phenomenon. The authors of [ST08] present an improved mechanism, capable of handling the inversion of a tetrahedral structure for finite element simulations. This process relies on finding a rotation that best aligns the deformed cell with its undeformed counterpart and then deriving penalty forces. While we could have used a similar approach, the above modification works for cases where negative volumes are reported. The elastic penalty forces described by equation (14) are capable of acting against the inversion. Additionally, other constraint-based forces can be derived (e.g. area preserving forces for the triangular faces of a tetrahedron or angle preserving forces).

4 Deformable linear object model

Using constrained mass point configuration, we can now describe the steps required to build a tetrahedral cell-based deformable object along the geometric image of a support curve:

  1. A curve discretization: where is a sample 3D point.

  2. A set of frames: define the orientation vectors , and . At each point , a coordinate frame is attached.

  3. Volumetric cells: for each pair of neighbouring vertices, , three tetrahedra are constructed: , , and (as depicted in figure 2 ).

    Figure 2: Tetrahedral cell division of a DLO segment.

The tetrahedral cells allow retrieving consistent information about local torsion and curvature changes from one link to another adjacent segment. A direct reference for how the object twists around the line is given by the relative orientation of to . Structural resistance is added by substituting the edges of the tetrahedra with linear springs. These constraints tend to act like curvature springs since they oppose the bending of the object around the line.

To account for plausible twisting behaviour, we introduce a quaternion based constraint system that acts like a torsional spring at each node of the object. Considering there are three connected nodes, and , in this order, the torsional spring tends to reposition the point such that the resulting configuration is closer to the initial, rest configuration. Such a process requires finding a suitable axis and computing the relative orientation of the vector with respect to the and vectors. We achieve this behaviour by computing a torsion compensating quaternion:


where , . The angle represents the angle between the and vectors with respect to the axis , and the superscript designates values for the initial, undeformed state of the object. The torsion quaternion, can then be used to rotate the vector to minimize the torsion offset. Instead of directly rotating this vector, a force will be applied to the node, thus mimicking the effects of an angular spring. Figure 3 depicts how the torsion quaternions are derived with respect to local geometry.

Figure 3: Torsion compensation using quaternions.

5 Simulation update logic

In general, the steps required for a complete update of the object’s state can be assembled as follows:

  • Numerical integration: compute the current acceleration from the current state , and add any force contributions from the collision solver stage. Using an explicit integration method, the new state is computed.

  • Approximate forces: . These values are to be used in the collision response stage.

  • Correct positions: applying a position based dynamics constraint projection, the linear springs are additionally relaxed. Such a process is equivalent to an iterative Gauss-Seidel solver and it requires translating the object vertices by small displacements. This step helps satisfying stiff constraints and avoids outstretching artifacts.

  • Collision resolution: pairwise link collisions are detected and response forces and impulses are derived. The positions and velocities are corrected, storing response forces in accumulators, . These force residues are fed back to the numerical integration scheme and are included in the acceleration component the next iteration will use.

These steps must then be performed in this order, leading to a plausible update effect at interactive frame-rates.

6 Results

We have used our model to simulate a simple laparoscopic suturing task. Given a physical substrate, a thread was driven through a tissue-like structure (see figure 4). The wire model we used is ideal for testing the behaviour of an integrator where both stiff and soft constraints drive the simulation. A simple methodology was employed to compute a score sheet for each integration scheme. We tracked the speed and stability of several methods while varying the time-step. The results were quantized by considering the explicit Euler method as a reference. The order of accuracy is not particularly important as it is orthogonal to our stability goals. Due to their energy-conserving properties, symplectic integrators are favoured when competing with other methods that achieve similar scores.

(a) Driving a wire through a severed tissue layer
(b) Tightening the suturing wire
Figure 4: Suturing simulation
Method Max Time Symplectic Order Score
Explicit Euler 0.0098s 8ms NO 1.225
Symplectic Euler 0.0294s 8ms YES 3.675
Midpoint 0.0153s 8ms NO 1.9125
Half-Step 0.0168s 10ms NO 1.68
Verlet 0.023s 10ms YES 2.3
Forest-Ruth 0.021s 12ms YES 1.75
Symplectic Midpoint 0.028s 8ms YES 3.5
Runge-Kutta 4 0.027s 14ms NO 1.928
Moified Half-Step 0.0216s 10ms NO 2.16
Table 1: Integrator benchmark results

Analyzing table 1, the symplectic Euler method (discussed by Cromer [Cro81]) is the most stable for the total iteration time it needs (the score is obtained as the ratio between the maximum time step and the time needed to execute one simulation iteration). We modified the Midpoint method to achieve symplecticity, obtaining the second highest score. Given the fact that the Midpoint family of methods is accurate to the second order for the position terms (while being a first order method for the velocity terms), we recommend it for applications where accuracy is of some importance. The Verlet method ( [Ver67] ), popular for molecular dynamics simulations, also benefits from its relatively high stability, symplecticity and second order accuracy. A third order method, the Forest-Ruth integration scheme [FR90]

, represents the best alternative for applications where accuracy is a key element. The Runge-Kutta fourth order method, although supporting relatively high time steps, is the most time consuming and probably not a good choice for real time applications.

Although symplectic integrators excel in scenarios where their energy preserving features are central (e.g. where only conservative forces are involved), we have found this family of integrators to outperform their explicit integrator counterparts. As a last remark, we have also modified the Half Step method ([HNW93]

) to support a semi-implicit update for the middle estimation

. This modification significantly increases the method’s stability, as seen in table 1.

As an additional benchmark study, we have used the pendulum equation, . The explicit Euler is clearly the most unstable, introducing ghost energies as seen in the phase space diagram comparison with respect to the midpoint method in figure 4(a). On the other hand, in figure 4(b), the symplectic Euler and Midpoint methods have a much higher stability range, with the latter yielding slightly lower energy variations. As second order explicit integrators, the modified half step method is also more stable than the original version (as depicted in figure 5(a)). Nevertheless, both methods tend to add ghost energies, but are much more stable than the first order explicit Euler. As the order of accuracy increases, larger time steps can be used (as it is the case with Runge-Kutta methods), but the performance impact does not justify such a trade-off. Finally, for models where no damping or non-conservative forces are involved, the Forest-Ruth and Verlet methods are the natural choices (comparative phase space plots in figure 5(b)). Nevertheless, in case accuracy can be sacrificed, we still recommend using the first order symplectic Euler method as it has the highest stability versus complexity score.

(a) Explicit Euler (blue) and Midpoint (yellow),
(b) Symplectic Midpoint (blue) and Symplectic Euler (yellow),
Figure 5: First order integrators (phase space diagrams)
(a) Half Step (blue) and Modified Half Step (yellow),
(b) Symplectic Euler (red), Verlet (yellow), Forest-Ruth (blue),
Figure 6: Higher order integrators (phase space diagrams)

7 Conclusion

In this work we have discussed the implications of employing an explicit integration method for updating a soft or semi-rigid body simulation. For real time applications where accuracy is not a goal, we recommend using a symplectic integrator as it has the best performance and stability score. Even when non-conservative forces are involved, this class of integrators is able to cope with stiff constraints. However, for applications where accuracy cannot be sacrificed, both the Runge-Kutta or the Forest-Ruth integrators can be used. The latter choice is accurate up to the third order and conserves energy, while the Runge-Kutta progressively loses small amounts of energy, counting for a slight increase in stability.

On a final note, implicit integrators, mentioned in the related work section 2, are not the usual choice for real time applications, thus motivating our investigation towards an explicit integrator alternative.

Appendix A Explicit Integration Methods

In the following appendix, we present the minor modifications of the explicit integration methods that were tested in our simulation application. These numerical methods make use of an acceleration function, , a fixed time-step, , and compute new positions and velocities, given the previous state .

Symplectic Midpoint Modified Half Step
Table 2: Slightly modified integrators for improved stability and accuracy
(a) Creating complex knots: tight knots tend correspond to very high elastic potentials. The numerical stability is crucial for plausible knot behaviour.
(b) Highly elastic cable unknotting itself under the action of strong constraint forces
Figure 7: Hose object knot tying


  • [BMOT13] Jan Bender, Matthias Müller, Miguel A. Otaduy, and Matthias Teschner. Position-based methods for the simulation of solid objects in computer graphics. In EUROGRAPHICS 2013 State of the Art Reports. Eurographics Association, 2013.
  • [BW98] David Baraff and Andrew Witkin. Large steps in cloth simulation. In Proceedings of the 25th annual conference on Computer graphics and interactive techniques, SIGGRAPH ’98, pages 43–54, New York, NY, USA, 1998. ACM.
  • [CM97] Lee Cooper Cooper and Steve Maddock. Preventing collapse within mass-spring-damper models of deformable objects. In In The Fifth Internation Conference in Central Europe on Computer Graphics and Visualization, pages 70–78, 1997.
  • [Cro81] Alan Cromer. Stable solutions using the euler approximation. American Journal of Physics, 49(5):455–459, 1981.
  • [DSB99] Mathieu Desbrun, Peter Schröder, and Alan Barr. Interactive animation of structured deformable objects. In Proceedings of the 1999 conference on Graphics interface ’99, pages 1–8, San Francisco, CA, USA, 1999. Morgan Kaufmann Publishers Inc.
  • [FR90] E. Forest and R. D. Ruth. Fourth-order symplectic integration. Phys. D, 43(1):105–117, April 1990.
  • [FSH10] Basil Fierz, Jonas Spillmann, and Matthias Harders. Stable explicit integration of deformable objects by filtering high modal frequencies. Journal of WSCG, 18(1-3):81–88, February 2010.
  • [HE01] Michael Hauth and Olaf Etzmuß. A high performance solver for the animation of deformable objects using advanced numerical methods. Comput. Graph. Forum, 20(3):319–328, 2001.
  • [HNW93] E. Hairer, S. P. Nørsett, and G. Wanner.

    Solving ordinary differential equations I (2nd revised. ed.): nonstiff problems

    Springer-Verlag New York, Inc., New York, NY, USA, 1993.
  • [Jim12] P. Jiménez. Survey on model-based manipulation planning of deformable objects. Robot. Comput.-Integr. Manuf., 28(2):154–163, April 2012.
  • [MGC07] Jaruwan Mesit, Ratan K. Guha, and Shafaq Chaudhry. 3d soft body simulation using mass-spring system with internal pressure force and simplified implicit integration. JCP, 2(8):34–43, 2007.
  • [NMK06] Andrew Nealen, Matthias Müller, Richard Keiser, Eddy Boxerman, and Mark Carlson. Physically based deformable models in computer graphics. Computer Graphics Forum, 25(4):809–836, 2006.
  • [SKZF11] Craig A. Schroeder, Nipun Kwatra, Wen Zheng, and Ronald Fedkiw. Asynchronous evolution for fully-implicit and semi-implicit time integration. Comput. Graph. Forum, 30(7):1983–1992, 2011.
  • [SLM06] Martin Servin, C. Lacoursiére, and N. Melin. Interactive simulation of elastic deformable materials. In SIGRAD 2006 Conference Proceedings, pages 22–32, 2006.
  • [ST08] Ruediger Schmedding and Matthias Teschner. Inversion handling for stable deformable modeling. Vis. Comput., 24(7):625–633, July 2008.
  • [THMG04] Matthias Teschner, Bruno Heidelberger, Matthias Muller, and Markus Gross. A versatile and robust model for geometrically complex deformable solids. In Proceedings of the Computer Graphics International, CGI ’04, pages 312–319, Washington, DC, USA, 2004. IEEE Computer Society.
  • [TKL04] Shan-Ho Tsai, M. Krech, and D.P. Landau. Symplectic integration methods in molecular and spin dynamics simulations. Brazilian Journal of Physics, 34:384 – 391, 06 2004.
  • [Ver67] Loup Verlet. Computer ”experiments” on classical fluids. i. thermodynamical properties of lennard-jones molecules. Phys. Rev., 159:98–103, Jul 1967.
  • [VMT05] P. Volino and N. Magnenat-Thalmann. Accurate garment prototyping and simulation. Computer-Aided Design & Applications, 2(5):645–654, 2005.