In order to integrate the D-NURBS ordinary differential equations of motion (17) in an interactive modeling environment, it is important to provide the modeler or designer with visual feedback about the evolving state of the dynamic model. Rather than using costly time integration methods that take the largest possible time steps, it is more crucial to provide a smooth animation by maintaining the continuity of the dynamics from one step to the next. Hence, less costly yet stable time integration methods that take modest time steps are desirable.

The matrices , , and (and , , and ) are symmetric, sparse, and banded. Several
algorithms are available for the numerical integration of the D-NURBS
ordinary differential equations of motion. The suitability of
implicit or explicit integration algorithms is dependent on the
bandwidth of the matrices, as determined by the dimensionality of the
parametric space and the order of the NURBS basis functions. The
matrices for a D-NURBS curve have a single band which has a
half-bandwidth of 4*k*, where *k* is the order of the NURBS basis.
For D-NURBS surfaces, the matrices become block banded, with each
block containing *n* bands similar to those of dynamic curves, where
*n* depends on the order of the NURBS basis in the opposite parametric
direction.

We integrate the differential equations (17)
through time by discretizing the derivative of over
time-steps . The state of the D-NURBS at time
is integrated using prior states at time *t* and .
Depending on the choice of physical parameters,
(17) may be a stiff system. We use an implicit
time integration method in order to maintain the stability of the
integration scheme. The implicit method employs discrete derivatives
of using backward differences

Making use of the fact that , we obtain the time integration formula

where the superscripts denote evaluation of the quantities at the indicated times, and where the remaining quantities are evaluated at time . For example, we can extrapolate the mass matrix using the formula

and likewise for the other matrices and vectors in (29). The simpler, constant extrapolations , etc., ([15] Section 8.6) also work satisfactorily.

In the interest of efficiency, we do not factorize the matrix expression on the left hand side of (29) in order to solve for . Instead, we employ the conjugate gradient method to obtain an iterative solution [27, 32]. To achieve interactive simulation rates, we limit the number of conjugate gradient iterations per time step to 10. We have observed that 2 iterations typically suffice to converge to a residual of less than . More than 2 iterations tend to be necessary when the physical parameters (mass, damping, tension, stiffness, applied forces) are changed dramatically during interactive sculpting.

Note that when physical parameter values are chosen such that the
equations (17) are not stiff, it is much cheaper
to employ an explicit time integration method using forward
differences. Appendix D discusses the forward
difference approach. Note that the explicit method requires values
for the matrices only at time *t*, hence (30)
is not needed.

For the D-NURBS curve, we simply replace with in (29) and everything proceeds as in the case of surfaces.

In the case of D-NURBS with linear constraints, we discretize the derivatives of (rather than ). Analogous to (29), the discrete version of (22) is

Since there are fewer degrees of freedom in than in , faster numerical implementation of constrained D-NURBS is possible, provided the constraint matrix is sparse. Note that since the conjugate gradient algorithm requires only gradient vectors, we need not compute , and explicitly. The only extra cost is the computation of and the multiplication of with several vectors in (31).

For nonlinear constraints, at each time step we can apply the
conjugate gradient algorithm to solve (27) for the
Lagrange multipliers and the constrained generalized
accelerations (given known and ).
We then integrate and from *t* to to obtain the constrained generalized velocities and
coordinates (e.g., using the simple Euler method ;
).