This class herits from the PHSCore class. It implements the symbolic expression of several numerical methods for the simulation of PHSCore. Below are the additional class attributes and class methods that differ from the parent PHSCore.

Recall of the Port-Hamiltonian structure

$$ \left( \begin{array}{c} \frac{\mathrm{d}\mathbf{x}_{\mathrm{}}}{\mathrm{d} t} \\ \mathbf{w}_{\mathrm{}} \\ \mathbf{y}_{\mathrm{}} \\ \end{array} \right) = \underbrace{\left( \begin{array}{ccc} \mathbf{M}_{\mathrm{xx}} & \mathbf{M}_{\mathrm{xw}} & \mathbf{M}_{\mathrm{xy}} \\ \mathbf{M}_{\mathrm{wx}} & \mathbf{M}_{\mathrm{ww}} & \mathbf{M}_{\mathrm{wy}} \\ \mathbf{M}_{\mathrm{yx}} & \mathbf{M}_{\mathrm{yw}} & \mathbf{M}_{\mathrm{yy}} \\ \end{array} \right)}_{\mathbf M} \cdot \left( \begin{array}{c} \nabla\mathrm{H}\\ \mathbf{z}_{\mathrm{}}\\ \mathbf{u}_{\mathrm{}}\\ \end{array} \right). \qquad (1) $$

The quantites considered as arguments are

• the states $\mathbf{x}(t)$,
• the dissipation variables $\mathbf{w}(t)$,
• the inputs $\mathbf{u}(t)$,
• the control parmameters $\mathbf{p}(t)$
• the observed quantites $\mathbf{o}(t)$

The significances of the control parameters and observed quantites are described latter (consider $\mathbf{p}=\mathbf{o}=\emptyset$ for the moment).

The system functions are

• the Hamiltonian's gradient $\nabla \mathrm H \equiv \mathrm \nabla \mathrm H(\mathbf{x})$,
• the structure matrix $\mathbf M \equiv \mathbf M(\mathbf{x}, \mathbf{w}, \mathbf{u}, \mathbf{p}, \mathbf{o})$,
• the dissipation function $\mathbf z \equiv \mathbf z(\mathbf{x}, \mathbf{w}, \mathbf{u}, \mathbf{p}, \mathbf{o})$.

Discrete state derivative and update

In the sequel, we denote $\mathbf x^k \equiv \mathbf x(t^k)$ on the regular temporal grid $t^k = k\,\mathrm{T_s}$ for the sample period $\mathrm{T_s}$ with associated sample rate $\mathrm{F_s}=\mathrm{T}_\mathrm{s}^{-1}$. Then, a one-stage progressive finite difference approximation of the state derivative is $$\frac{\mathrm{d}\mathbf{x}_{\mathrm{}}(t^k)}{\mathrm{d} t} \simeq \frac{\delta\mathbf{x}^k_{\mathrm{}}}{\mathrm{T_s}}, \qquad (2)$$ where $$\delta\mathbf{x}^k_{\mathrm{}} \triangleq \mathbf{x}_{\mathrm{}}^{k+1} - \mathbf{x}_{\mathrm{}}^{k}, \qquad (3)$$ Now, the update is simply given by $$\mathbf{x}_{\mathrm{}}^{k+1} = \mathbf{x}_{\mathrm{}}^{k}+\delta\mathbf{x}^k_{\mathrm{}}. \qquad (4)$$

Discrete arguments

Quantites considered as arguments are now

• the states $\mathbf{x}^k$,
• the states increments $\delta\mathbf{x}^k$,
• the dissipation variables $\mathbf{w}^k$,
• the inputs $\mathbf{u}^k$,
• the control parmameters $\mathbf{p}^k$,
• the observed quantites $\mathbf{o}^k$,

which are collected in a single vector $$ \mathbf{args^k}_{\mathrm{}} \triangleq \left( \begin{array}{c} \mathbf{x}^k\\ \delta\mathbf{x}^k\\ \mathbf{w}^k\\ \mathbf{u}^k\\ \mathbf{p}^k\\ \mathbf{o}^k \end{array} \right).\qquad \begin{array}{l} \mathrm{States}\\ \mathrm{States~increments}\\ \mathrm{Dissipation~variables}\\ \mathrm{Inputs}\\ \mathrm{Control~parameters}\\ \mathrm{Observed~quantities} \end{array} $$ The significances of the control parameters and observed quantites are described latter (consider $\mathbf{p}^k=\mathbf{o}^k=\emptyset$ for the moment).

Numerical methods

A numerical method consists in the prescription of a numerical equivalent to system (1), that is, in the prescription of numerical equivalents for

• the Hamiltonian's gradient $\nabla \mathrm H \leftarrow \overline{\nabla}\mathrm{H}(\mathbf{args}_{\mathrm{}})$,
• the structure matrix $\mathbf M \leftarrow \overline{\mathbf M}(\mathbf{args}_{\mathrm{}})$,
• the dissipation function $\mathbf z \leftarrow \overline{\mathbf z}(\mathbf{args}_{\mathrm{}})$.

This yields the following discrete system

$$ \left( \begin{array}{c} \frac{\delta\mathbf{x}^k_{\mathrm{}}}{\mathrm{T_s}} \\ \mathbf{w}^k_{\mathrm{}} \\ \mathbf{y}^k_{\mathrm{}} \end{array} \right) = \overline{\mathbf{M}}(\mathbf{args}^k_{\mathrm{}}) \cdot \left( \begin{array}{c} \overline{\nabla}\mathrm{H}(\mathbf{args}^k_{\mathrm{}})\\ \overline{\mathbf z}(\mathbf{args}^k_{\mathrm{}})\\ \mathbf{u}^k_{\mathrm{}} \end{array} \right). \qquad (5) $$

Now, provided

• a state initialization {$\mathbf{x}^k\}_{k = 0}$,
• an input sequence {$\mathbf{u}^k\}_{0\leq k \leq K}$,
• a control sequence {$\mathbf{p}^k\}_{0\leq k \leq K}$,
• an observed sequence {$\mathbf{o}^k\}_{0\leq k \leq K}$,

the simulation is performed by solving (5) iteratively for the state increment $\delta\mathbf{x}^k$ and the dissipation variable $\mathbf{w}^k$ (which are the only unknown in $\mathbf{args}^k$), and then updating according to (4).

The Theta method

For a standard dynamical system, this method reads $$ \dot{\mathbf{x}}=\mathbf{f}(\mathbf{x}) \longrightarrow \frac{\delta\mathbf{x}^{k}}{\mathrm{T_s}} = \mathbf{f}(\mathbf{x}^k + \theta \, \delta\mathbf{x}^k).$$ It includes the explicit (respectively implicit) Euler scheme for $\theta = 0$ (respectively $\theta = 1$), and the midpoint rule for $\theta = \mbox{½}$.

Applying this scheme to the system (1) yields the following numerical equivalents:

• $\overline{\nabla}\mathrm{H}(\mathbf{args}^k_{\mathrm{}}) \triangleq {\nabla}\mathrm{H}(\mathbf x^k+ \theta \, \delta\mathbf{x}^k)$,
• $\overline{\mathbf M}(\mathbf{args}^k_{\mathrm{}}) \triangleq \mathbf M(\mathbf x^k+ \theta \, \delta\mathbf{x}^k, \mathbf{w}^k, \mathbf{u}^k, \mathbf{p}^k, \mathbf{o}^k)$,
• $\overline{\mathbf z}(\mathbf{args}^k_{\mathrm{}}) \triangleq \mathbf z(\mathbf x^k+ \theta \, \delta\mathbf{x}^k, \mathbf{w}^k, \mathbf{u}^k, \mathbf{p}^k, \mathbf{o}^k)$.

The Trapezoidal method

For a standard dynamical system, this method reads $$ \dot{\mathbf{x}} = \mathbf{f}(\mathbf{x}) \longrightarrow \frac{\delta\mathbf{x}^{k}}{\mathrm{T_s}} = \frac{1}{2} \Big(\mathbf{f}(\mathbf{x}^k) + \mathbf{f}(\mathbf{x}^k+ \delta\mathbf x^k)\Big).$$

Applying this scheme to the system (1) yields the following numerical equivalents:

• $\overline{\nabla}\mathrm{H}(\mathbf{args}^k_{\mathrm{}}) \triangleq \frac{1}{2} \Big({\nabla}\mathrm{H}(\mathbf x^k) + {\nabla}\mathrm{H}(\mathbf x^k+ \delta\mathbf{x}^k)\Big)$,
• $\overline{\mathbf M}(\mathbf{args}^k_{\mathrm{}}) \triangleq \frac{1}{2} \Big(\mathbf M(\mathbf{x}^k, \mathbf{w}^k, \mathbf{u}^k, \mathbf{p}^k, \mathbf{o}^k) + \mathbf M(\mathbf{x}^k+ \delta\mathbf x^k, \mathbf{w}^k, \mathbf{u}^k, \mathbf{p}^k, \mathbf{o}^k)\Big)$,
• $\overline{\mathbf z}(\mathbf{args}^k_{\mathrm{}}) \triangleq \frac{1}{2} \Big(\mathbf z(\mathbf{x}^k, \mathbf{w}^k, \mathbf{u}^k, \mathbf{p}^k, \mathbf{o}^k) + \mathbf z(\mathbf{x}^k+ \delta\mathbf x^k, \mathbf{w}^k, \mathbf{u}^k, \mathbf{p}^k, \mathbf{o}^k)\Big)$.

The Discrete Gradient method

For a standard Hamiltonian system with a separate Hamiltonian $\mathrm H(\mathbf x) = \sum_{i=1}^{\mathrm{dim}(\mathbf x)} \mathrm H_i(x_i)$, this method reads $$ \dot{\mathbf{x}} = \mathbf J(\mathbf x)\cdot \nabla\mathrm{H}(\mathbf{x}) \longrightarrow \frac{\delta\mathbf{x}^{k}}{\mathrm{T_s}} = \mathbf J(\mathbf x^k+\theta\, \delta\mathbf x^k)\cdot \overline \nabla \mathrm H (\mathbf x^k, \delta \mathbf x^k),$$ with $$ \left[\overline \nabla \mathrm H (\mathbf x, \delta \mathbf x)\right]_i = \left\{ \begin{array}{l} \frac{\mathrm H_i(x_i+\delta x_i)- \mathrm H_i (x_i)}{\delta x_i} \text{ if } \delta x_i \geq 0, \\ \mathrm H_i^\prime(x_i) \mbox{ else}, \end{array} \right. \quad \mbox{for } {1\leq i \leq \mathrm{dim}(x)}. \qquad (6) $$

Applying this scheme to the system (1) yields the following numerical equivalents:

• $\overline{\nabla}\mathrm{H}(\mathbf{args}^k_{\mathrm{}})$ is given by ( 6),
• $\overline{\mathbf M}(\mathbf{args}^k_{\mathrm{}}) \triangleq \mathbf M(\mathbf x^k+ \theta \, \delta\mathbf{x}^k, \mathbf{w}^k, \mathbf{u}^k, \mathbf{p}^k, \mathbf{o}^k)$ (same as Theta method),
• $\overline{\mathbf z}(\mathbf{args}^k_{\mathrm{}}) \triangleq \mathbf z(\mathbf x^k+ \theta \, \delta\mathbf{x}^k, \mathbf{w}^k, \mathbf{u}^k, \mathbf{p}^k, \mathbf{o}^k)$ (same as Theta method).

Resolution

Arguments

The values for

• the states $\mathbf{x}^k$,
• the inputs $\mathbf{u}^k$,
• the control parmameters $\mathbf{p}^k$,
• the observed quantites $\mathbf{o}^k$,

are known at the begining of each time step. These quantities are constant over a time-step, and are collected in a singe vector: $$ \mathbf{c}_{\mathrm{}} \triangleq \left( \begin{array}{c} \mathbf{x} \\ \mathbf{u} \\ \mathbf{p} \\ \mathbf{o} \\ \end{array} \right). $$ The only unknown quantities in (5) and $\mathbf{args}$ are

• the states increments $\delta\mathbf{x}^k$,
• the dissipation variables $\mathbf{w}^k$.

These unknown quantities are collected in a singe vector: $$ \mathbf{v}_{\mathrm{}} \triangleq \left( \begin{array}{c} \delta\mathbf{x} \\ \mathbf{w}_{\mathrm{}} \end{array} \right),\qquad (7) $$ and $\mathbf{args} \equiv (\mathbf{v}, \mathbf{c})$.

Update

Similarly, the constitutive relations are collected in a single vector valued function: $$ \mathbf{f}_{\mathrm{}}(\mathbf{v}, \mathbf{c}) \triangleq \left( \begin{array}{c} \overline \nabla\mathrm{H}(\mathbf{v}, \mathbf{c})\\ \overline{\mathbf{z}}_{\mathrm{}}(\mathbf{v}, \mathbf{c}) \end{array} \right), \qquad (8) $$ and we define the following matrices (without expressing the dependance on $\mathbf{args} \equiv (\mathbf{v}, \mathbf{c})$): $$ \begin{array}{llc} \overline{\mathbf{M}}_{\mathrm{vv}}= \left( \begin{array}{ccc} \overline{\mathbf{M}}_{\mathrm{xx}} & \overline{\mathbf{M}}_{\mathrm{xw}} \\ \overline{\mathbf{M}}_{\mathrm{wx}} & \overline{\mathbf{M}}_{\mathrm{ww}} \end{array} \right), & \overline{\mathbf{M}}_{\mathrm{vy}} = \left( \begin{array}{ccc} \overline{\mathbf{M}}_{\mathrm{xy}}\\ \overline{\mathbf{M}}_{\mathrm{wy}} \end{array} \right), & (9) \\ \overline{\mathbf{M}}_{\mathrm{yv}} = \left( \begin{array}{ccc} \overline{\mathbf{M}}_{\mathrm{yx}} & \overline{\mathbf{M}}_{\mathrm{yw}} \end{array} \right), & \overline{\mathbf{I}}= \left( \begin{array}{ccc} \mathrm{F_s}\,\mathbf{I_d}(n_\mathbf{x}) & \mathbf 0 \\ \mathbf 0 & \mathbf{I_d}(n_\mathbf{w}) \end{array} \right), & (10) \end{array} $$ where $n_\mathbf{x}=\mathrm{dim}(\mathbf x)$ and $n_\mathbf{w}=\mathrm{dim}(\mathbf w)$. Then, system (2) is rewritten as

1. Solve $\mathbf{F}(\mathbf v^k, \mathbf c^k)=\mathbf 0 $ for $\mathbf v^k$, where $$\mathbf{F}(\mathbf v^k, \mathbf c^k) \triangleq \overline{\mathbf{I}}\cdot\mathbf v^k - \mathbf M_{\mathrm{vv}}\cdot\mathbf{f}_{\mathrm{}}(\mathbf v^k, \mathbf c^k) - \mathbf M_{\mathrm{vy}}\cdot\mathbf u^k, \qquad (11)$$
2. Set $$\mathbf y^k = \mathbf M_{\mathrm{yv}}\cdot\mathbf{f}_{\mathrm{}}(\mathbf{v}^k, \mathbf c^k) + \mathbf M_{\mathrm{yy}}\cdot\mathbf u^k, \qquad (12)$$
3. Set $$\mathbf x^{k+1} = \mathbf x^k + \delta\mathbf x^k. \qquad (13)$$

Explicit and implicit updates

Depending on the original system (1) and on the chosen numerical method, some components of the implicit function $\mathbf{F}(\mathbf v)$ in (11) may be explicited up to a matrix inversion. In order to detect these explicit components, we analyse the structure of the jacobian matrix of $\mathcal J_{\mathbf{F},\mathbf v}(\mathbf v)$ defined by $$\Big[\mathcal J_{\mathbf{F},\mathbf v}(\mathbf v, \mathbf c)\Big]_{i,j} = \frac{\partial \mathbf F_i (\mathbf v, \mathbf c)}{\partial \mathbf v_j}. $$ The result is the following spliting: $$\mathbf x = \left( \begin{array}{c} \mathbf x_\mathrm{E} \\ \mathbf x_\mathrm{I} \end{array} \right), \quad \mathbf w = \left( \begin{array}{c} \mathbf w_\mathrm{E} \\ \mathbf w_\mathrm{I} \end{array} \right), \qquad (14) $$ and $$\mathbf v_\mathrm{E} = \left( \begin{array}{c} \delta \mathbf x_\mathrm{E} \\ \mathbf w_\mathrm{E} \end{array} \right), \quad \mathbf v_\mathrm{I} = \left( \begin{array}{c} \delta \mathbf x_\mathrm{I} \\ \mathbf w_\mathrm{I} \end{array} \right), \qquad (15) $$ such that the columns of $\mathcal J_{\mathbf{F}, \mathbf v_\mathrm{E}}(\mathbf v, \mathbf c)$ associated with $\mathbf v_\mathrm{E}$ does no depend on $\mathbf v\equiv(\mathbf v_\mathrm{E}, \mathbf v_\mathrm{I})$.

According to (7-11), we define the following expressions (without expressing the dependance on $\mathbf{args} \equiv (\mathbf{v}, \mathbf{c})$): $$ \mathbf{f}_\mathrm{E} \triangleq \left( \begin{array}{c} \overline \nabla\mathrm{H}_\mathrm{E}\\ \overline{\mathbf{z}}_\mathrm{E} \end{array} \right), \quad \mathbf{f}_\mathrm{I} \triangleq \left( \begin{array}{c} \overline \nabla\mathrm{H}_\mathrm{I}\\ \overline{\mathbf{z}}_\mathrm{I} \end{array} \right)\qquad (16) $$ $$ \begin{array}{llc} \overline{\mathbf{M}}_{\mathrm{vEvE}} = \left( \begin{array}{ccc} \overline{\mathbf{M}}_{\mathrm{xExE}} & \overline{\mathbf{M}}_{\mathrm{xEwE}} \\ \overline{\mathbf{M}}_{\mathrm{wExE}} & \overline{\mathbf{M}}_{\mathrm{wEwE}} \end{array} \right), & \overline{\mathbf{M}}_{\mathrm{vIvI}} = \left( \begin{array}{ccc} \overline{\mathbf{M}}_{\mathrm{xIxI}} & \overline{\mathbf{M}}_{\mathrm{xIwI}} \\ \overline{\mathbf{M}}_{\mathrm{wIxI}} & \overline{\mathbf{M}}_{\mathrm{wIwI}} \end{array} \right), & (17) \\ \overline{\mathbf{M}}_{\mathrm{vEvI}} = \left( \begin{array}{ccc} \overline{\mathbf{M}}_{\mathrm{xExI}} & \overline{\mathbf{M}}_{\mathrm{xEwI}} \\ \overline{\mathbf{M}}_{\mathrm{wExI}} & \overline{\mathbf{M}}_{\mathrm{wEwI}} \end{array} \right), & \overline{\mathbf{M}}_{\mathrm{vIvE}} = \left( \begin{array}{ccc} \overline{\mathbf{M}}_{\mathrm{xIxE}} & \overline{\mathbf{M}}_{\mathrm{xIwE}} \\ \overline{\mathbf{M}}_{\mathrm{wIxE}}& \overline{\mathbf{M}}_{\mathrm{wIwE}} \end{array} \right), & (18) \\ \overline{\mathbf{M}}_{\mathrm{vEy}}= \left( \begin{array}{ccc} \overline{\mathbf{M}}_{\mathrm{xEy}}\\ \overline{\mathbf{M}}_{\mathrm{wEy}} \end{array} \right), & \overline{\mathbf{M}}_{\mathrm{vIy}}= \left( \begin{array}{ccc} \overline{\mathbf{M}}_{\mathrm{xIy}}\\ \overline{\mathbf{M}}_{\mathrm{wIy}} \end{array} \right), & (19) \\ \overline{\mathbf{I}}_{\mathrm{E}}= \left( \begin{array}{ccc} \mathrm{F_s}\,\mathbf{I_d}(n_{\mathbf{x}_{\mathrm E}}) & \mathbf 0 \\ \mathbf 0 & \mathbf{I_d}(n_{\mathbf{w}_{\mathrm E}}) \end{array} \right), & \overline{\mathbf{I}}_{\mathrm{I}}= \left( \begin{array}{ccc} \mathrm{F_s}\,\mathbf{I_d}(n_{\mathbf{x}_{\mathrm I}}) & \mathbf 0 \\ \mathbf 0 & \mathbf{I_d}(n_{\mathbf{w}_{\mathrm I}}) \end{array} \right), & (20) \end{array} $$ where $n_{\mathbf{x}_{\mathrm E}}=\mathrm{dim}(\mathbf x_{\mathrm E})$, $n_{\mathbf{x}_{\mathrm I}}=\mathrm{dim}(\mathbf x_{\mathrm I})$, $n_{\mathbf{w}_{\mathrm I}}=\mathrm{dim}(\mathbf w_{\mathrm I})$ and $n_{\mathbf{w}_{\mathrm E}}=\mathrm{dim}(\mathbf w_{\mathrm E})$, and $$\mathbf{F}_{\mathrm{E}}(\mathbf v, \mathbf{c}) \triangleq \overline{\mathbf{I}}_{\mathrm{E}}\cdot\mathbf v_{\mathrm{E}} - \mathbf M_{\mathrm{vEvE}}\cdot\mathbf{f}_{\mathrm{E}}(\mathbf v, \mathbf{c})- \mathbf M_{\mathrm{vEvI}}\cdot\mathbf{f}_{\mathrm{I}}(\mathbf v, \mathbf{c}) - \mathbf M_{\mathrm{vEy}}\cdot\mathbf u, \qquad (21)$$ $$\mathbf{F}_{\mathrm{I}}(\mathbf v, \mathbf{c}) \triangleq \overline{\mathbf{I}}_{\mathrm{I}}\cdot\mathbf v_{\mathrm{I}} - \mathbf M_{\mathrm{vIvE}}\cdot\mathbf{f}_{\mathrm{E}}(\mathbf v, \mathbf{c})- \mathbf M_{\mathrm{vIvI}}\cdot\mathbf{f}_{\mathrm{I}}(\mathbf v, \mathbf{c}) - \mathbf M_{\mathrm{vEy}}\cdot\mathbf u. \qquad (22)$$

Now, define $\mathbf{G}(\mathbf v, \mathbf c)$ as the local deviation of $\mathbf{F}(\mathbf v, \mathbf c)$ from its linear behavior w.r.t $\mathbf v_\mathrm{E}$: $$\mathbf{G}(\mathbf v, \mathbf c) \triangleq \mathbf{F}(\mathbf v, \mathbf c) - \mathcal J_{\mathbf{F},\mathbf v_{\mathrm{E}}}(\mathbf v, \mathbf c) \cdot\mathbf v_{\mathrm{E}}, \qquad (23)$$ and correspondingly $$\mathbf{G}_{\mathrm{E}}(\mathbf v, \mathbf c) \triangleq \mathbf{F}_{\mathrm{E}}(\mathbf v, \mathbf c) - \mathcal J_{\mathbf{F}_{\mathrm{E}}, \mathbf v_{\mathrm{E}}}(\mathbf v, \mathbf c)\cdot\mathbf v_{\mathrm{E}}, \qquad (24)$$ $$\mathbf{G}_{\mathrm{I}}(\mathbf v, \mathbf c) \triangleq \mathbf{F}_{\mathrm{I}}(\mathbf v, \mathbf c) - \mathcal J_{\mathbf{F}_{\mathrm{I}}, \mathbf v_{\mathrm{E}}}(\mathbf v, \mathbf c) \cdot\mathbf v_{\mathrm{E}}. \qquad (25)$$

Notice that by definition $\mathcal J_{\mathbf{F}_{\mathrm{E}}, \mathbf v_{\mathrm{E}}}(\mathbf v, \mathbf c)\equiv \mathcal J_{\mathbf{F}_{\mathrm{E}}, \mathbf v_{\mathrm{E}}}(\mathbf c)$ is constant over a time-step. This yields $\mathbf{G}_{\mathrm{E}}(\mathbf v, \mathbf c) \equiv \mathbf{G}_{\mathrm{E}}(\mathbf v_{\mathrm{I}}, \mathbf c) $ for the root $\mathbf{F}_{\mathrm{E}}(\mathbf v^*, \mathbf c)=\mathbf 0$: $$\mathbf v_{\mathrm{E}}^* = -\mathcal J_{\mathbf{F}_{\mathrm{E}}, \mathbf v_{\mathrm{E}}}^{-1}(\mathbf c) \cdot \mathbf{G}_{\mathrm{E}}(\mathbf v_{\mathrm{I}}^*, \mathbf c).\qquad (26)$$

Now, rewrite $\mathbf{F}_{\mathrm{I}}(\mathbf v, \mathbf c)$ as $$\begin{array}{rcl} \mathbf{F}_{\mathrm{I}}(\mathbf v, \mathbf c) &=& \mathbf{G}_{\mathrm{I}}(\mathbf v_{\mathrm{I}}, \mathbf c) + \mathcal J_{\mathbf{F}_{\mathrm{I}}, \mathbf v_{\mathrm{E}}}(\mathbf c)\cdot\mathbf v_{\mathrm{E}} \\ &=&\mathbf{G}_{\mathrm{I}}(\mathbf v_{\mathrm{I}}, \mathbf c) - \mathcal J_{\mathbf{F}_{\mathrm{I}}, \mathbf v_{\mathrm{E}}}(\mathbf c)\cdot\mathcal J_{\mathbf{F}_{\mathrm{E}}, \mathbf v_{\mathrm{E}}}^{-1}(\mathbf c) \cdot \mathbf{G}_{\mathrm{E}}(\mathbf v_{\mathrm{I}}, \mathbf c)\equiv \mathbf{F}_{\mathrm{I}}(\mathbf v_{\mathrm{I}}, \mathbf c), \end{array}$$ with $$ \mathcal J_{\mathbf{F}_{\mathrm{I}}, \mathbf v_{\mathrm{I}}}(\mathbf v_{\mathrm{I}}) = \mathcal J_{\mathbf{G}_{\mathrm{I}}, \mathbf v_{\mathrm{I}}}(\mathbf v_{\mathrm{I}}) - \mathcal J_{\mathbf{F}_{\mathrm{I}}, \mathbf v_{\mathrm{E}}}\cdot\mathcal J_{\mathbf{F}_{\mathrm{E}}, \mathbf v_{\mathrm{E}}}^{-1} \cdot \mathcal J_{\mathbf{G}_{\mathrm{E}}, \mathbf v_{\mathrm{I}}}(\mathbf v_{\mathrm{I}}). $$

Then, the resolution algorithm for the system (2) is

1. $\mathbf x^k \leftarrow \mathbf x^{-1} + \delta\mathbf x^{k-1}$,
2. $\mathbf c^k \leftarrow (\mathbf x^k, \mathbf u^k, \mathbf p^k, \mathbf o^k)$,
3. $\mathbf v_{\mathrm{I}}^k\leftarrow \mathbf v_{\mathrm{I}}^*$ such that $\mathbf{F}_{\mathrm{I}}(\mathbf v_{\mathrm{I}}^*, \mathbf c^k)=0$,
4. $\mathbf v_{\mathrm{E}}^k\leftarrow-\mathcal J_{\mathbf{F}_{\mathrm{E}}, \mathbf v_{\mathrm{E}}}^{-1}(\mathbf c^k) \cdot \mathbf{G}_{\mathrm{E}}(\mathbf v_{\mathrm{I}}^k, \mathbf c^k)$,
5. $\mathbf y^k\leftarrow\mathbf M_{\mathrm{yv}}\cdot\mathbf{f}_{\mathrm{}}(\mathbf{v}^k, \mathbf c^k) + \mathbf M_{\mathrm{yy}}\cdot\mathbf u^k$.