How does constraint resolution affect the stability/accuracy of numerical integration?

Question

I understand some basic analysis techniques (local truncation error, global error, zero-stable, absolute stable, etc.) of numerical integration. But I find it hard to apply these techniques in practice.

In my case (a hair/cloth simulation for computer games), I usually do a numerical integration step followed by constraints resolving steps in every frame. I wonder how does constraint resolution affect the numerical integration?

Particularly, I am concerned with the following situation: in most cases the exact solution of constraints cannot be reached due to numerical error, insufficiently many iterative steps, etc.

When an error is introduced during constraint resolution, what will happen to the numerical integration in the next frame? Does this error act as external sources (inhomogeneous part) to the differential equations system?

---

EDIT: more background

I am using Verlet integration in my simulation. To resolve constraints I directly modify the positions. In the next temporal integration step the error in positions would contribute "ghost" velocity to the system.

My intention is to analyze the impact of this error and find a way to control the error to make my temporal integration stable and accurate to some desired order.

Answer 1

One way to look at it is to say that you are interested in solving $$ y'(t) = f(y(t)), \qquad g(y(t)) = 0, $$ and you only approximately enforce the constraints, thus solving $$ y'(t) = f(y(t)), \qquad g(y(t)) = O(\epsilon). $$

If the constraints are well-behaved (say differentiable, although that's a bit much maybe), then (by the implicit function theorem) $g(y(t)) = O(\epsilon)$ implies that $y(t)$ is "close" to a true solution $ y(t) = \tilde y(t) + \delta y(t)$ with $g(\tilde y(t)) = 0$ and the error being on the order of $$ |\delta y(t)| = O(|\nabla g^{-1}|\epsilon). $$ (The inverse of possibly-singular Jacobian is, say, in the sense of least squares.)

Thus if both approximate $y$ and the true $\tilde y$ satisfy the same ODE, then $$ (y-\tilde y)' = f(y) - f(\tilde y) \sim (\nabla f)(y-\tilde y). $$ Conversely, solving the approximate system is like solving the perturbed system $$ y' = f(y) + (\nabla f)\delta y $$ with $\delta y$ being an unknowable function on the order of $|\nabla g^{-1}|\epsilon$.

So what matters is not the numerical integrator you use, but the sensitivity of your ODE to small perturbations. In other words, if the integrator solves the right equation correctly, then it is only the properties of the ODE that should matter. What you care about is whether the error $|y-\tilde y|$ stays small over longer time periods.

If the ODE is stable, i.e., solutions starting with similar initial conditions remain close in some sense, the effect of the error will be on the order of $O(|\nabla f||\nabla g^{-1}|\epsilon)$, so hopefully small enough. If the ODE is backward stable, i.e., for any approximate solution you can find a "nearby" initial condition that would have produced the true solution, then in that sense the errors would be small enough.

In both cases, the sensitivity of the ODE to perturbations is the issue, not the numerical method. So I would say it should be sufficient to check for any approximate solution $y$ whether there is indeed a "nearby" true solution $\tilde y$, which can be done by analyzing the constraint function. You can, of course, come up with very badly behaved constraints.

In general, perturbation theory of ODEs is a very well known subject.

Answer 2

Most of the DAE solvers I'm aware of don't take advantage of symplectic structure, but it's entirely possible that one of them does.

Typically, if you use a DAE solver, you are using an implicit method, or maybe a semi-implicit method; the algebraic variables are determined through a nonlinear solve, and I don't see how you get around that. I don't know how you enforce constraints. If you modify the positions via a nonlinear solve, then a semi-implicit DAE solver might be competitive; if modifying positions is cheaper than a nonlinear solve, then the DAE solver is probably more expensive. A fully implicit DAE solver is more likely to be slower than your current approach; however, DAE solvers have solid convergence theory behind them (see Hairer and Wanner's classic text; Brenan and Petzold's book, and the book by Roswitha Marz), which would put the numerical error analysis on more solid footing.

Other approaches to constraint satisfaction also exist. You could try looking at what the molecular dynamics community does to enforce constraints in their mechanical systems, such as via Lagrange multiplies, or specialized algorithms geared towards enforcing things like bond lengths and bond angles of atoms in molecules. It sounds like what you're doing is solving a projected dynamical system, and you want to be careful with that projection operation. I remember one of my advisors discussing how projection can affect stability, such that perturbations off of the constraint manifold due to, well, floating-point arithmetic, can yield unstable solutions because the underlying ODE system (before projection) is unstable in the off-manifold directions. You could also look at differential-variational inequalities. Both variational inequality solvers and DAE solvers are available in PETSc; you can find an implicit DAE solver in SUNDIALS.

One thing you want to watch out for with a mechanical system is the index of your DAE. (There are many definitions of index; I'm just going to assume "index" means "differential index", even though I'm vaguely aware there are pitfalls in using that definition.) Mechanical systems may have multiple DAE formulations, and some of them are high-index (in practice, this means the index is greater than 2). Most DAE solvers tend to solve index-1 or semi-explicit Hessenberg index-2 DAEs; if you have a system not in that form, you'll need a specialized solver (there are a few index-3 solvers specialized for mechanical systems) or you'll need to reformulate your problem via something like a change in coordinates.