Dirichlet Enforcement for Ordinary Differential Equations

Unconstrained Ordinary Differential Equations

Consider a system described by then following linear, second order differential equation:

\[\begin{split}\begin{gather} \textbf{M} \; \ddot{\textbf{x}}(t) + \textbf{C} \; \dot{\textbf{x}}(t) + \textbf{K} \; \textbf{x}(t) = \textbf{0} \\ \textbf{x}(0) = \textbf{x}_0 \\ \dot{\textbf{x}}(0) = \dot{\textbf{x}}_0 \end{gather}\end{split}\]

In order to numerically integrate this differential equation, we can start by expressing the the future state as an extrapolation of the current state (e.g. using backward Euler):

\[\begin{split}\begin{gathered} \begin{cases} \textbf{x}(t + \Delta t) = \textbf{x}(t) + \dot{\textbf{x}}(t) \Delta t + \ddot{\textbf{x}}(t + \Delta t) \Delta t^2 \\ \dot{\textbf{x}}(t + \Delta t) = \dot{\textbf{x}}(t) + \ddot{\textbf{x}}(t + \Delta t) \Delta t \end{cases} \end{gathered}\end{split}\]

Substituting this assumption back in to the differential equation evaluated at time \(t + \Delta t\) produces an expression that relates a future state: \(\{\textbf{x}(t + \Delta t), \dot{\textbf{x}}(t + \Delta t)\}\) to the current state \(\{\textbf{x}(t), \dot{\textbf{x}}(t)\}\):

\[\textbf{M} \; \ddot{\textbf{x}}(t+\Delta t) + \textbf{C} \, \big(\dot{\textbf{x}}(t) + \ddot{\textbf{x}}(t + \Delta t) \Delta t\big) + \textbf{K} \; \big(\textbf{x}(t) + \dot{\textbf{x}}(t) \Delta t + \ddot{\textbf{x}}(t + \Delta t) \Delta t^2\big) = \textbf{0}\]

After rearranging terms to put the known quantities on the righthand side, we are left with

\[\big(\textbf{M} + \textbf{C} \Delta t + \textbf{K} \Delta t^2\big) \; \ddot{\textbf{x}}(t+\Delta t) = - \textbf{C} \; \dot{\textbf{x}}(t) - \textbf{K} \big(\textbf{x}(t) + \dot{\textbf{x}}(t) \Delta t\big)\]

From here, provided that the coefficient matrix is invertible, a system of equations can be solved to determine \(\ddot{\textbf{x}}(t+\Delta t)\) and advance the solution state to time \(t + \Delta t\).

Imposing Single Degree-of-Freedom Constraints

Now, let us consider how the solution process changes if some components of \(\textbf{x}(t)\) are controlled externally. In this case, our differential equation looks like

\[\begin{split}\begin{gather} \textbf{M} \; \ddot{\textbf{x}}(t) + \textbf{C} \; \dot{\textbf{x}}(t) + \textbf{K} \; \textbf{x}(t) = \textbf{0} \\ \color{red}{\textbf{G} \; \textbf{x}(t) = \textbf{g}(t)} \\ \textbf{x}(0) = \textbf{x}_0 \\ \dot{\textbf{x}}(0) = \dot{\textbf{x}}_0 \end{gather}\end{split}\]

where \(\textbf{G}\) is a matrix with 1 on the diagonals corresponding to constrained components, and 0 elsewhere, and \(\textbf{g}(t)\) is the vector of prescribed values for the constrained degrees of freedom. There are several equivalent ways to interpret the additional constraint equations (we assume the constraints are consistent with the initial conditions):

\[\textbf{G} \; \textbf{x}(t) = \textbf{g}(t) \; \iff \; \textbf{G} \; \dot{\textbf{x}}(t) = \dot{\textbf{g}}(t) \; \iff \; \textbf{G} \; \ddot{\textbf{x}}(t) = \ddot{\textbf{g}}(t)\]

At the end of the day, we are solving for a single acceleration vector, \(\ddot{\textbf{x}}(t+\Delta t)\), and as a result, we cannot hope to simultaneously satisfy more than one interpretation of these constraints. The different Dirichlet Enforcement Methods in serac relate to the following different interpretations of the constraint equations.

"Direct Control"

Direct control proceeds to solve for accelerations by effectively replacing the constrained rows of the differential equation by the constraint equations,

\[(\textbf{1} - \textbf{G})(\textbf{M} \; \ddot{\textbf{x}}(t+\Delta t) + \textbf{C} \; \dot{\textbf{x}}(t+\Delta t) + \textbf{K} \; \textbf{x}(t+\Delta t)) + \textbf{G} \; \textbf{x}(t+\Delta t)\;=\;\textbf{g}(t+\Delta t)\]

Now, like before, we substitute in expressions for extrapolated state and solve for \(\ddot{\textbf{x}}(t+\Delta t)\). In this interpretation of the constraint, the resulting accelerations ensure that \(\textbf{G} \; \textbf{x}(t+\Delta t) = \textbf{g}(t+\Delta t)\) is satisfied.

"Rate Control"

Rate control uses a similar idea, but with the rate-form interpretation of the constraint:

\[(\textbf{1} - \textbf{G})(\textbf{M} \; \ddot{\textbf{x}}(t+\Delta t) + \textbf{C} \; \dot{\textbf{x}}(t+\Delta t) + \textbf{K} \; \textbf{x}(t+\Delta t)) + \textbf{G} \; \dot{\textbf{x}}(t+\Delta t)\;=\;\dot{\textbf{g}}(t+\Delta t)\]

Here, the end-step accelerations are determined such that the derivative of the constraint equations is exactly satisfied, \(\textbf{G} \; \dot{\textbf{x}}(t+\Delta t) = \dot{\textbf{g}}(t+\Delta t)\). This tends to have better stability than the "Direct Control" option, but is prone to solution "drift" over long periods of time.

"Full Control"

This last option starts by additively decomposing the solution vector into two parts: constrained and unconstrained:

\[\begin{split}\begin{cases} \textbf{x}(t) = \textbf{x}_c(t) + \textbf{x}_u(t) = \textbf{G} \; \textbf{x}(t) + (\textbf{1} - \textbf{G}) \; \textbf{x}(t) \\ \dot{\textbf{x}}(t) = \dot{\textbf{x}}_c(t) + \dot{\textbf{x}}_u(t) = \textbf{G} \; \dot{\textbf{x}}(t) + (\textbf{1} - \textbf{G}) \; \dot{\textbf{x}}(t)\\ \ddot{\textbf{x}}(t) = \ddot{\textbf{x}}_c(t) + \ddot{\textbf{x}}_u(t) = \textbf{G} \; \ddot{\textbf{x}}(t) + (\textbf{1} - \textbf{G}) \; \ddot{\textbf{x}}(t) \end{cases}\end{split}\]

From here, the constrained terms \(\{\textbf{x}_c(t), \dot{\textbf{x}}_c(t), \ddot{\textbf{x}}_c(t)\}\) are replaced by their prescribed values, \(\{\textbf{g}(t), \dot{\textbf{g}}(t), \ddot{\textbf{g}}(t)\}\) and only the unconstrained acceleration components, \(\ddot{\textbf{x}}_u(t+\Delta t)\), are solved for. For situations where \(\textbf{g}(t)\) has a discontinuous derivative, this approach may be preferable to Rate Control.

Note: time derivatives of \(\textbf{g}(t)\) are currently computed by central finite difference stencils

\[\begin{split}\begin{gather} \dot{\textbf{g}}(t) \approx \frac{\textbf{g}(t + \epsilon) - \textbf{g}(t - \epsilon)}{2 \; \epsilon} \\ \ddot{\textbf{g}}(t) \approx \frac{\textbf{g}(t + \epsilon) - 2 \, \textbf{g}(t) +\textbf{g}(t - \epsilon)}{\epsilon^2} \end{gather}\end{split}\]