21  Projection with implicit numerical scheme

Author

Affiliations

Enrico Bertolazzi

University of Trento, Department of Industrial Engineering

21.1 Index-1 DAE

Consider a semi-explicit DAE (with index > 1):

\left\{\begin{aligned} \bm{z}' &= \bm{f}(\bm{z}, \bm{w}, t) \\[0.5em] \bm{0} &= \bm{g}(\bm{z}, \bm{w}, t) \end{aligned}\right.

This system can be transformed into an index-1 DAE by index reduction, along with invariants:

\textrm{Index-1 DAE:}\left\{\begin{aligned} \bm{x}' &= \bm{f}(\bm{x}, \bm{y}, t) \\[0.5em] \bm{0} &= \bm{a}(\bm{x}, \bm{y}, t) \end{aligned}\right. \qquad \underbrace{\bm{b}(\bm{x}, \bm{y}, t) = \bm{0}}_{\textrm{invariants}}

The index-1 system is then solved as an implicitly defined ODE:

\bm{x}' = \bm{F}(\bm{x}, t)

where \bm{F}(\bm{x}, t) = \bm{f}(\bm{x}, \bm{y}(\bm{x}, t), t) and \bm{y}(\bm{x}, t) is the solution to the equation \bm{b}(\bm{x}, \bm{y}, t) = \bm{0}.

Advancing the ODE: Explicit Methods

Explicit methods, such as the explicit Euler method, do not require complex handling. The update step is given by:

\bm{x}_{k+1} = \bm{x}_k + \Delta t \bm{F}(\bm{x}_k, t_k)

The evaluation of \bm{F}(\bm{x}_k, t_k) proceeds as follows:

1. Solve \bm{h}({\color{red}\bm{y}}) = \bm{a}(\bm{x}_k, {\color{red}\bm{y}}, t_k) = \bm{0} to compute \bm{y}_k.

2. Compute \bm{F}(\bm{x}_k, t_k) = \bm{f}(\bm{x}_k, \bm{y}_k, t_k).

Advancing the ODE: Implicit Methods

For implicit methods, such as implicit Euler, the update step is:

\bm{x}_{k+1} = \bm{x}_k + \Delta t \bm{F}(\bm{x}_{k+1}, t_{k+1})

This requires solving the implicit equation:

\bm{G}({\color{red}\bm{z}}) = {\color{red}\bm{z}} - \Delta t \bm{F}({\color{red}\bm{z}}, t_{k+1}) - \bm{x}_k = \bm{0}

The solution of \bm{G}({\color{red}\bm{z}}) = 0 requires an iterative method such as Newton-Raphson:

{\color{red}\bm{z}}^{(\ell+1)} = {\color{red}\bm{z}}^{(\ell)} - \left(\frac{\partial \bm{G}({\color{red}\bm{z}}^{(\ell)})}{\partial {\color{red}\bm{z}}}\right)^{-1} \bm{G}({\color{red}\bm{z}}^{(\ell)})

This method requires the Jacobian of \bm{F}(\bm{x}, t):

\frac{\partial \bm{G}({\color{red}\bm{z}})}{\partial {\color{red}\bm{z}}} = \bm{I} - \Delta t \frac{\partial \bm{F}({\color{red}\bm{z}}, t_{k+1})}{\partial {\color{red}\bm{z}}}

RHS Jacobian

Given the system

\bm{F}(\bm{x}, t) = \bm{f}(\bm{x}, \bm{y}(\bm{x}, t), t), \qquad \bm{a}(\bm{x}, \bm{y}(\bm{x}, t), t) = \bm{0},

we can derive the following expressions for the Jacobian:

\left\{ \begin{aligned} \frac{\partial \bm{F}(\bm{x}, t)}{\partial \bm{x}} &= \frac{\partial \bm{f}(\bm{x}, \bm{y}, t)}{\partial \bm{x}} + \frac{\partial \bm{f}(\bm{x}, \bm{y}, t)}{\partial \bm{y}} \frac{\partial \bm{y}(\bm{x}, t)}{\partial \bm{x}} \\[1em] \bm{0} &= \frac{\partial \bm{a}(\bm{x}, \bm{y}, t)}{\partial \bm{x}} + \frac{\partial \bm{a}(\bm{x}, \bm{y}, t)}{\partial \bm{y}} \frac{\partial \bm{y}(\bm{x}, t)}{\partial \bm{x}} \end{aligned} \right.

From this, we can express the derivative of \bm{y}(\bm{x}, t) with respect to \bm{x} as:

\frac{\partial \bm{y}(\bm{x}, t)}{\partial \bm{x}} = -\bm{a}_{\bm{y}}(\bm{x}, \bm{y}, t)^{-1} \bm{a}_{\bm{x}}(\bm{x}, \bm{y}, t),

and the Jacobian of \bm{F}(\bm{x}, t) becomes:

\boxed{ \frac{\partial \bm{F}(\bm{x}, t)}{\partial \bm{x}} = \bm{f}_{\bm{x}}(\bm{x}, \bm{y}, t) - \bm{f}_{\bm{y}}(\bm{x}, \bm{y}, t) \bm{a}_{\bm{y}}(\bm{x}, \bm{y}, t)^{-1} \bm{a}_{\bm{x}}(\bm{x}, \bm{y}, t) }

Practical Computation of the Jacobian

In practice, the Jacobian should be computed as follows:

1. Start with the system:

   \bm{F}(\bm{x}, t) = \bm{f}(\bm{x}, \bm{y}(\bm{x}, t), t), \qquad \bm{a}(\bm{x}, \bm{y}(\bm{x}, t), t) = \bm{0}.

2. Solve for \bm{Z} by solving the linear system:

   \bm{a}_{\bm{y}}(\bm{x}, \bm{y}, t) \bm{Z} = - \bm{a}_{\bm{x}}(\bm{x}, \bm{y}, t)

3. Compute the Jacobian:

   \frac{\partial \bm{F}(\bm{x}, t)}{\partial \bm{x}} = \bm{f}_{\bm{x}}(\bm{x}, \bm{y}, t) + \bm{f}_{\bm{y}}(\bm{x}, \bm{y}, t) \bm{Z}

21.2 Index-3 DAE Special Case

We begin with the following system of differential-algebraic equations (DAEs):

\left\{ \begin{aligned} \bm{p}' &= \bm{v} \\[0.5em] \bm{M}(\bm{p}, \bm{v}, t) \bm{v}' + \bm{\Phi}_{\bm{p}}(\bm{p}, t)^T \bm{\lambda} &= \bm{f}(\bm{p}, \bm{v}, t) \\[0.5em] \bm{\Phi}(\bm{p}, t) &= \bm{0} \end{aligned} \right.

NoteTransformation to ODE with Implicit RHS and Invariants

We aim to transform this system into a set of ordinary differential equations (ODEs) with implicit right-hand sides (RHS) and invariants:

\left\{ \begin{aligned} \bm{p}' &= \bm{v} \\[0.5em] \bm{v}' &= \bm{v}_{\mathrm{dot}}(\bm{p}, \bm{v}, t) \end{aligned} \right. \qquad \left\{ \begin{aligned} \color{blue}\bm{0} &= \color{blue} \bm{h}(\bm{p}, \bm{v}, t) = \frac{\mathrm{d}}{\mathrm{d}t} \bm{\Phi}(\bm{p}, t) \\[0.5em] \color{red}\bm{0} &= \color{red} \bm{\Phi}(\bm{p}, t) \end{aligned} \right.

Here, \bm{v}_{\mathrm{dot}} is derived from the system by solving the following linear system:

\begin{pmatrix} \bm{M}(\bm{p}, \bm{v}, t) & \bm{\Phi}_{\bm{p}}(\bm{p}, t)^T \\[1em] \bm{\Phi}_{\bm{p}}(\bm{p}, t) & \bm{0} \end{pmatrix} \begin{pmatrix} \bm{v}_{\mathrm{dot}} \\[1em] \bm{\lambda} \end{pmatrix} = \begin{pmatrix} \bm{f}(\bm{p}, \bm{v}, t) \\[1em] -\bm{a}(\bm{p}, \bm{v}, t) \end{pmatrix}

Additionally, we define:

\bm{a}(\bm{p}, \bm{v}, t) = \bm{h}_{\bm{p}}(\bm{p}, \bm{v}, t) \bm{v} + \bm{h}_t(\bm{p}, \bm{v}, t)

Defining Functions for RHS and Invariants

Let us introduce the following functions to represent the system:

• Function F:

\bm{F}(\bm{p}, \bm{v}, {\color{red}\bm{v}_{\mathrm{dot}}}, {\color{red}\bm{\lambda}}, t) = \bm{M}(\bm{p}, \bm{v}, t) {\color{red}\bm{v}_{\mathrm{dot}}} + \bm{\Phi}_{\bm{p}}(\bm{p}, t)^T {\color{red}\bm{\lambda}} - \bm{f}(\bm{p}, \bm{v}, t)

• Function G:

\bm{G}(\bm{p}, \bm{v}, {\color{red}\bm{v}_{\mathrm{dot}}}, t) = \bm{\Phi}_{\bm{p}}(\bm{p}, t) {\color{red}\bm{v}_{\mathrm{dot}}} + \bm{a}(\bm{p}, \bm{v}, t)

Thus, the system can be expressed as:

\begin{aligned} \bm{F}(\bm{p}, \bm{v}, {\color{red}\bm{v}_{\mathrm{dot}}(\bm{p}, \bm{v}, t)}, {\color{red}\bm{\lambda}(\bm{p}, \bm{v}, t)}, t) &= \bm{0} \\[1em] \bm{G}(\bm{p}, \bm{v}, {\color{red}\bm{v}_{\mathrm{dot}}(\bm{p}, \bm{v}, t)}, t) &= \bm{0} \end{aligned}

Jacobian Derivatives of Functions F and G

Taking derivatives of \bm{F} and \bm{G} with respect to \bm{p}, we obtain the following system of equations:

\begin{aligned} \bm{F}_{\bm{p}}(\bm{p}, \bm{v}, {\color{red}\bm{v}_{\mathrm{dot}}}, {\color{red}\bm{\lambda}}, t) + \bm{F}_{\bm{v}_{\mathrm{dot}}}(\bm{p}, \bm{v}, {\color{red}\bm{v}_{\mathrm{dot}}}, {\color{red}\bm{\lambda}}, t) {\color{blue}\frac{\partial \bm{v}_{\mathrm{dot}}}{\partial \bm{p}}} + \bm{F}_{\bm{\lambda}}(\bm{p}, \bm{v}, \bm{v}_{\mathrm{dot}}, \bm{\lambda}, t) {\color{blue}\frac{\partial \bm{\lambda}}{\partial \bm{p}}} &= \bm{0} \\ \bm{G}_{\bm{p}}(\bm{p}, \bm{v}, {\color{red}\bm{v}_{\mathrm{dot}}}, t) + \bm{G}_{\bm{v}_{\mathrm{dot}}}(\bm{p}, \bm{v}, {\color{red}\bm{v}_{\mathrm{dot}}}, t) {\color{blue}\frac{\partial \bm{v}_{\mathrm{dot}}}{\partial \bm{p}}} &= \bm{0} \end{aligned}

Expressions for the Jacobians of Functions F and G

From the definition of the functions, we have the following expressions for the Jacobian derivatives:

\begin{aligned} \bm{F}_{\bm{v}_{\mathrm{dot}}}(\bm{p}, \bm{v}, \bm{v}_{\mathrm{dot}}, \bm{\lambda}, t) &= \bm{M}(\bm{p}, \bm{v}, t) \\[1em] \bm{F}_{\bm{\lambda}}(\bm{p}, \bm{v}, \bm{v}_{\mathrm{dot}}, \bm{\lambda}, t) &= \bm{\Phi}_{\bm{p}}(\bm{p}, t)^T \\[1em] \bm{G}_{\bm{v}_{\mathrm{dot}}}(\bm{p}, \bm{v}, \bm{v}_{\mathrm{dot}}, t) &= \bm{\Phi}_{\bm{p}}(\bm{p}, t) \end{aligned}

Final System of Equations for the Jacobians

The system of Jacobian equations can now be written as:

\begin{pmatrix} \bm{M}(\bm{p}, \bm{v}, t) & \bm{\Phi}_{\bm{p}}(\bm{p}, t)^T \\[1em] \bm{\Phi}_{\bm{p}}(\bm{p}, t) & \bm{0} \end{pmatrix} \frac{\partial}{\partial \bm{p}} \begin{pmatrix} \bm{v}_{\mathrm{dot}} \\[1em] \bm{\lambda} \end{pmatrix} = \begin{pmatrix} -\bm{F}_{\bm{p}}(\bm{p}, \bm{v}, \bm{v}_{\mathrm{dot}}, \bm{\lambda}, t) \\[1em] -\bm{G}_{\bm{p}}(\bm{p}, \bm{v}, \bm{v}_{\mathrm{dot}}, t) \end{pmatrix}

Similarly, for derivatives with respect to \bm{v}:

\begin{pmatrix} \bm{M}(\bm{p}, \bm{v}, t) & \bm{\Phi}_{\bm{p}}(\bm{p}, t)^T \\[1em] \bm{\Phi}_{\bm{p}}(\bm{p}, t) & \bm{0} \end{pmatrix} \frac{\partial}{\partial \bm{v}} \begin{pmatrix} \bm{v}_{\mathrm{dot}} \\[1em] \bm{\lambda} \end{pmatrix} = \begin{pmatrix} -\bm{F}_{\bm{v}}(\bm{p}, \bm{v}, \bm{v}_{\mathrm{dot}}, \bm{\lambda}, t) \\[1em] -\bm{G}_{\bm{v}}(\bm{p}, \bm{v}, \bm{v}_{\mathrm{dot}}, t) \end{pmatrix}

Final Transformed Index-3 DAE

Thus, the transformed system of index-3 DAEs is:

\left\{ \begin{aligned} \bm{p}' &= \bm{v} \\[0.5em] \bm{v}' &= \bm{v}_{\mathrm{dot}}(\bm{p}, \bm{v}, t) \end{aligned} \right. \qquad \left\{ \begin{aligned} \color{blue} \bm{0} &= \color{blue} \bm{h}(\bm{p}, \bm{v}, t) = \frac{\mathrm{d}}{\mathrm{d}t} \bm{\Phi}(\bm{p}, t) \\[0.5em] \color{red} \bm{0} &= \color{red} \bm{\Phi}(\bm{p}, t) \end{aligned} \right.

The RHS and its Jacobians are obtained from:

\begin{aligned} \begin{pmatrix} \bm{M}(\bm{p}, \bm{v}, t) & \bm{\Phi}_{\bm{p}}(\bm{p}, t)^T \\[1em] \bm{\Phi}_{\bm{p}}(\bm{p}, t) & \bm{0} \end{pmatrix} \begin{pmatrix} \bm{v}_{\mathrm{dot}} \\[1em] \bm{\lambda} \end{pmatrix} & = \begin{pmatrix} \bm{f}(\bm{p}, \bm{v}, t) \\[1em] -\bm{a}(\bm{p}, \bm{v}, t) \end{pmatrix} \\[3em] \begin{pmatrix} \bm{M}(\bm{p}, \bm{v}, t) & \bm{\Phi}_{\bm{p}}(\bm{p}, t)^T \\[1em] \bm{\Phi}_{\bm{p}}(\bm{p}, t) & \bm{0} \end{pmatrix} \frac{\partial}{\partial \bm{p}} \begin{pmatrix} \bm{v}_{\mathrm{dot}} \\[1em] \bm{\lambda} \end{pmatrix} & = \begin{pmatrix} -\bm{F}_{\bm{p}}(\bm{p}, \bm{v}, \bm{v}_{\mathrm{dot}}, \bm{\lambda}, t) \\[1em] -\bm{G}_{\bm{p}}(\bm{p}, \bm{v}, \bm{v}_{\mathrm{dot}}, t) \end{pmatrix} \end{aligned}