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 4k, 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 
;
).
