[ type = algorithm , float , name = Algorithm , ]algorithm

# Operator Splitting Based Dynamic Iteration for Linear Port-Hamiltonian
Systems

Andreas Bartel1 Michael Günther1 Birgit Jacob1 and Timo Reis2
1 IMACM  Bergische Unviersität Wuppertal
Gaußstr. 20  42119 Wuppertal  Germany
2 Institut für Mathematik  Technische Universität Ilmenau
Weimarer Str. 25  98693 Ilmenau  Germany

Abstract. A dynamic iteration scheme for linear differential-algebraic port-
Hamiltonian systems based on Lions-Mercier-type operator splitting methods is
developed. The dynamic iteration is monotone in the sense that the error is
decreasing and no stability conditions are required. The developed iteration
scheme is even new for linear port-Hamiltonian systems. The obtained algorithm
is applied to multibody systems and electrical networks.

Keywords.  Port-Hamiltonian Systems, DAE, Operator Splitting, Dynamic
Iteration.

MSC Classification.  37Jxx, 34A09

##  1 Introduction

Dynamic iteration schemes for differential-algebraic equations (DAE) have been
widely used and discussed for multiphysics problems, which allow for an easy
exploitation of the different properties of the subystems [4, 1, 8, 7, 6, 3].
In general, we consider an initial value problem for a (semi-explicit index-1)
DAE on $t\in[0,T]$ of the form

$\displaystyle\dot{v}$ $\displaystyle=f(t,v,w),\quad v(0)=v_{0},$ (1a)
$\displaystyle 0$ $\displaystyle=g(t,v,w),$ (1b)

where $\partial g(t,v,w)/\partial w$ is regular in a neighborhood of the
solution $(v(t),w(t))^{\top}$, $t\in[0,T]$. A general dynamic iteration scheme
is given by an initial guess $(v^{(0)}(t),w^{(0)}(t))^{\top}$ ($t\in[0,T]$,
i.e., an initial waveform) and the solution $(v^{(k)}(t),w^{(k)}(t))^{\top}$
of the DAE for integer $k\geq 1$

$\displaystyle\dot{v}^{(k)}$
$\displaystyle=F(t,v^{(k)},v^{(k-1)},w^{(k)},w^{(k-1)}),\quad
v^{(k)}(0)=v_{0},$ $\displaystyle 0$
$\displaystyle=G(t,v^{(k)},v^{(k-1)},w^{(k)},w^{(k-1)})$

with arbitrary splitting functions $F,G$ fulfilling the compatibility
conditions

$F(t,v,v,w,w)=f(t,v,w)\quad\text{and}\quad G(t,v,v,w,w)=g(t,v,w).$

Hence a dynamic iteration defines as a mapping from the $(k-1)$-th iterate to
the $k$-th iterate:

$(v^{(k-1)}(t),\,w^{(k-1)}(t))\mapsto(v^{(k)}(t),w^{(k)}(t)).$

One drawback of such schemes for DAE is given by the fact that stability
conditions for a guaranteed convergence must hold in a twofold manner:

  1. (a)

convergence within one time window has to be present [4, 15], and

  2. (b)

stable error propagation from window to window is needed [4, 1].

If one of these conditions is not met choosing arbitrary small time windows
will not help out (in general).

In the case of initial value problems of ordinary differential equations

$\displaystyle\dot{x}$ $\displaystyle=f(t,x),\quad x(0)=x_{0},\quad
t\in[0,T],$

convergence of dynamic iteration schemes (analogously defined) is given for
arbitrary long time windows $[0,T]$. However, the convergence might not be
monotone, and the error might increase drastically in the beginning, causing
an overflow in computer implementations, see chapter 3 for an example.

In this work, we show that an alternative iterative approach for linear DAE
systems can avoid both flaws described above if the DAE system (1) is composed
of coupled port-Hamiltonian systems. This type of systems is motivated by
physics, provides an energy balance and gives a simple framework for coupled
systems [23, 16, 20, 21, 24, 25, 12, 13, 9]. Now, the alternative iterative
approach is based on an iteration scheme from Lions and Mercier [19]. The
convergence will be monotone, and no stability conditions appear.

The paper is organized as follows: in the subsequent chapter 2, the framework
of linear coupled port-Hamiltonian DAE systems is set. Both, the perspective
of an overall coupled system and the perspective of a system composed of
coupled port-Hamiltonian DAE subsystems is given. Chapter 3 illustrates the
limits of classical dynamic iteration schemes such as the Jacobi iteration,
which motivates the derivation and analysis of operator splitting based
dynamic iteration schemes in the following. Chapter 4 recapitulates some basic
facts on maximal monotone operators, which are used in chapter 5 to derive the
monotone convergence results of operator splitting based dynamic iteration
schemes. An estimate for the convergence rate can be given in the case of ODE
systems and DAE systems with enough dissipation. Numerical results are
discussed in chapter 6. The paper finishes with some concluding remarks and an
outlook.

##  2 The setting

We consider linear time-invariant initial value problems of DAEs, which have
the following structure:

######  Definition 1.

For a final time $T>0$, the descriptor variable $x(t)\in{}^{n}$ and
$t\in[0,T]$, we consider

$\displaystyle\tfrac{\mathrm{d}}{\mathrm{d}t}Ex(t)$
$\displaystyle=(J-R)Qx(t)+Bu(t),\quad t\in[0,T],\quad Ex(0)=Ex_{0},$ (2a)
$\displaystyle y(t)$ $\displaystyle=B^{\top}Qx(t),$ (2b)

where $E$, $J$, $R$ and $Q$ are real $n\times n$-matrices, $B$ is a $n\times
m$-matrix, $Ex_{0}\in{}^{n}$ is a given initial value, $u$ is a given input
and $y$ is referred to as output. (The matrix $E$ may have no full rank.)
$\Box$

Throughout this article, we are working with weak solutions.

######  Definition 2 (Solution).

Let $x_{0}\in{}^{n}$ and $u\in L^{2}([0,T],{}^{m})$. A function $x\in
L^{2}([0,T],{}^{n})$ is called (weak) solution of (2) if $Ex\in
H^{1}([0,T],{}^{n})$, $Ex(0)=Ex_{0}$ and (2) is satisfied for $t\in[0,T]$
almost everywhere. $\Box$

######  Proposition 3.

Assume that a differential-algebraic system (2) is given with
$E,J,R\in{}^{n\times n}$, $B\in{}^{n\times m}$ satisfying
$E^{\top}Q=Q^{\top}E\geq 0$, $R=R^{\top}\geq 0$ and $J^{\top}=-J$.
Furthermore, let $u\in L^{2}([0,T],{}^{m})$, $x_{0}\in{}^{n}$, and let $x\in
L^{2}([0,T],{}^{n})$ be a solution of (2). Then the following statements hold:

  1. (a)

The function $t\mapsto x(t)^{\top}Q^{\top}Ex(t)$ is weakly differentiable with

$\tfrac{\mathrm{d}}{\mathrm{d}t}\big{(}x^{\top}Q^{\top}Ex\big{)}=2x^{\top}Q^{\top}\big{(}\tfrac{\mathrm{d}}{\mathrm{d}t}Ex)\,\in
L^{1}([0,T],\real).$

  2. (b)

The following dissipation inequality holds for all $t\in[0,T]$:

$\tfrac{1}{2}\big{(}x(t)^{\top}Q^{\top}Ex(t)\big{)}-\tfrac{1}{2}\big{(}x_{0}^{\top}Q^{\top}Ex_{0}\big{)}\\\
=-\int_{0}^{t}x(\tau)^{\top}Q^{\top}RQx(\tau)d\tau+\int_{0}^{t}u(\tau)^{\top}y(\tau)d\tau\leq\int_{0}^{t}u(\tau)^{\top}y(\tau)d\tau.$

###### Proof.

  1. (a)

We have $E^{\top}Q=Q^{\top}E$, and thus, for the Moore-Penrose inverse $E^{+}$
of $E$, we have

$E^{\top}Q=E^{\top}(E^{+})^{\top}E^{\top}Q=E^{\top}(E^{+})^{\top}Q^{\top}E.$

Consequently, for $X=(E^{+})^{\top}Q^{\top}\in\mathbb{R}^{n\times n}$, we have
$Q^{\top}E=E^{\top}XE$. By using that $Ex\in H^{1}([0,T],\mathbb{R}^{n})$, the
product rule for weak derivatives [2, Thm. 4.25] implies that

$x^{\top}Q^{\top}Ex=(Ex)^{\top}X(Ex)\in W^{1,1}([0,T],\mathbb{R})$

with

$\displaystyle\tfrac{\mathrm{d}}{\mathrm{d}t}(x^{\top}Q^{\top}Ex)$
$\displaystyle=\tfrac{\mathrm{d}}{\mathrm{d}t}\left((Ex)^{\top}X(Ex)\right)=2x^{\top}E^{\top}X\tfrac{\mathrm{d}}{\mathrm{d}t}\left(Ex\right)=2x^{\top}\tfrac{\mathrm{d}}{\mathrm{d}t}\left(E^{\top}XEx\right)$
$\displaystyle=2x^{\top}\tfrac{\mathrm{d}}{\mathrm{d}t}\left(Q^{\top}Ex\right)=2x^{\top}Q^{\top}\tfrac{\mathrm{d}}{\mathrm{d}t}\left(Ex\right).$

  2. (b)

The previous statement yields

$\displaystyle\tfrac{\mathrm{d}}{\mathrm{d}t}\tfrac{1}{2}(x^{\top}Q^{\top}Ex)$
$\displaystyle=x^{\top}Q^{\top}\tfrac{\mathrm{d}}{\mathrm{d}t}(Ex)=x^{\top}Q^{\top}((J-R)Qx+Bu)$
$\displaystyle=-(Qx)^{\top}R(Qx)+x^{\top}QBu\leq(B^{\top}Qx)^{\top}u=y^{\top}u.$

Now an integration on $[0,t]$ leads to the dissipation inequality.

∎

Thus, we investigate the following class of port-Hamiltonian DAEs:

######  Definition 4 (PH-DAE).

The system (2) with the assumptions of Prop. 3 is here referred as port-
Hamiltonian DAE (PH-DAE). $\Box$

###  2.1 Perspective as overall coupled system

We further assume that (2) is composed of several port-Hamiltonian systems
which are coupled in an energy-preserving way. The assumptions are collected
in the following.

######  Assumption 5.

(Overall coupled PH-DAE) We consider the setting (2).

  1. (a)

The matrices in (2) are, for some $s,n_{1},\ldots,n_{s}\in\mathbb{N}$,
structured as

$\displaystyle E=\begin{bmatrix}E_{1}&&\\\ &\ddots\\\
&&E_{s}\end{bmatrix}\\!,\;\;R=\begin{bmatrix}R_{1}&&\\\ &\ddots\\\
&&R_{s}\end{bmatrix}\\!,\;\;Q=\begin{bmatrix}Q_{1}&&\\\ &\ddots\\\
&&Q_{s}\end{bmatrix}\\!,$ (3) $\displaystyle
J=\begin{bmatrix}J_{1}&J_{12}&\cdots&J_{1s}\\\
-J_{12}^{\top}&\ddots&\ddots&\vdots\\\ \vdots&\ddots&\ddots&J_{s-1,s}\\\
-J_{1s}^{\top}&\cdots&-J_{s-1,s}^{\top}&J_{s}\end{bmatrix}\\!,\;\;B=\begin{bmatrix}B_{1}\\\
\vdots\\\ \vdots\\\ B_{s}\end{bmatrix},$ (4)

with $n_{1}+n_{2}+\dotsc+n_{s}=n$, where

    1. (i)

for $i=1,\dotsc,\,s$, $E_{i},Q_{i},J_{i},R_{i}\in\mathbb{R}^{n_{i}\times
n_{i}}$ with $R_{i}=R_{i}^{\top}\geq 0$, $J_{i}^{\top}=-J_{i}$,
$E_{i}^{\top}Q_{i}=Q_{i}^{\top}E_{i}\geq 0$, $Q_{i}$ invertible and
$B_{i}\in\mathbb{R}^{n_{i}\times m}$;

    2. (ii)

$J_{ij}\in\mathbb{R}^{n_{i}\times n_{j}}$ for $i,j\in\\{1,\ldots,s\\}$ with
$i<j$.

  2. (b)

$\mathrm{rk}\begin{bmatrix}E&R&J\end{bmatrix}=n$.

  3. (c)

$T>0$, $x_{0}\in\mathbb{R}^{n}$ and $u\in L^{2}([0,T];\mathbb{R}^{m})$.
Furthermore, for matrices $Z\in\mathbb{R}^{n\times r}$,
$Z_{1}\in\mathbb{R}^{r\times r_{1}}$, with full column rank and

$\displaystyle\operatorname{im}Z=$ $\displaystyle\,\ker E,$
$\displaystyle\operatorname{im}Z_{1}=$ $\displaystyle\,\ker RQZ\cap\ker
Z^{\top}Q^{\top}JQZ,$

the input $u$ and the initial value $x_{0}$ fulfill

$\displaystyle Z_{1}^{\top}Z^{\top}Q^{\top}Bu\in$
$\displaystyle\,H^{1}([0,T];\mathbb{R}^{r_{1}})\quad\text{ with
}Z_{1}^{\top}Z^{\top}Q^{\top}(JQx_{0}+Bu(0))=0.$ $\Box$

######  Remark 6.

Notice that with Assumption 5, we have for (2):

  1. (a)

$Q^{\top}E=E^{\top}Q\geq 0,\quad R=R^{\top}\geq 0,\quad\mathrm{and}\quad Q$ is
invertible. Hence, the overall coupled PH-DAE of Assumption 5 is a PH-DAE (in
the sense of Definition 4).

  2. (b)

As $Q^{\top}E\geq 0$, we have that $EQ^{-1}=Q^{-\top}(Q^{T}E)Q^{-1}$ is
positive semi-definite as well, and therefore possesses a matrix square root
$(EQ^{-1})^{1/2}$.

  3. (c)

The $i$th subsystem for the variable $x_{i}$ reads:

$E_{i}x_{i}=(J_{i}-R_{i})Q_{i}x_{i}+B_{i}u\;-\sum_{j<i}J_{ji}^{\top}Q_{j}x_{j}+\sum_{i<j}J_{ij}Q_{j}x_{j}.$

The coupling is facilitated by the off-diagonal terms within the matrix $J$.
$\Box$

######  Remark 7 (Coupled linear PH-DAE systems).

A PH-DAE system of type (2) with coupling structure (3) and (4) is naturally
given in the case of $s$ coupled linear PH-DAE systems as follows: consider
$s$ PH-DAE subsystems

$\displaystyle\tfrac{\mathrm{d}}{\mathrm{d}t}E_{i}x_{i}(t)=$
$\displaystyle(\tilde{J}_{i}-R_{i})Q_{i}x_{i}(t)+B_{i}u_{i}(t),$
$\displaystyle y_{i}(t)=B_{i}^{\top}Q_{i}x_{i}(t),$

$i=1,\ldots,s.$ with $\tilde{J}_{i}=-\tilde{J}_{i}^{\top}$,
$R_{i}=R_{i}^{\top}\geq 0$, and $Q_{i}^{\top}E_{i}=E_{i}^{\top}Q_{i}\geq 0$.
The input $u_{i}$, the output $y_{i}$ and $B_{i}$ are split into

$u_{i}(t)=\begin{pmatrix}\hat{u}_{i}(t)\\\ \bar{u}_{i}(t)\end{pmatrix},\quad
y_{i}(t)=\begin{pmatrix}\hat{y}_{i}(t)\\\ \bar{y}_{i}(t)\end{pmatrix},\quad
B_{i}=\begin{pmatrix}\hat{B}_{i}&\bar{B}_{i}\end{pmatrix}$

according to external and coupling quantities. The subsystems are coupled via
external inputs and outputs by

$\displaystyle\begin{pmatrix}\hat{u}_{1}\\\ \vdots\\\
\hat{u}_{k}\end{pmatrix}+\hat{C}\begin{pmatrix}\hat{y}_{1}\\\ \vdots\\\
\hat{y}_{k}\end{pmatrix}=0,\qquad\hat{C}=-\hat{C}^{\top}.$

These $s$ systems can be condensed to one large PH-DAE system [8]

$\displaystyle\tfrac{\mathrm{d}}{\mathrm{d}t}Ex$
$\displaystyle=(J-R)Qx+\bar{B}\bar{u},$ (5a) $\displaystyle\bar{y}$
$\displaystyle=\bar{B}^{\top}Qx$ (5b)

with $J=\tilde{J}-\hat{B}\hat{C}\hat{B}^{\top}$ and the condensed quantities

$\displaystyle v^{\top}=(v_{1}^{\top},\ldots,v_{s}^{\top})\quad\text{ for
}v\in\\{x,\,\bar{u},\bar{y}\\},$ $\displaystyle
F=\operatorname{diag}\,(F_{1},\ldots,F_{s})\quad\text{ for
}F\in\\{E,\,Q,\,\tilde{J},\,R,\,\hat{B},\,\bar{B}\\}.$

Equation (5) defines now a PH-DAE system of type (2) with structure given by
(3) and (4): the matrices $\tilde{J}$ and $R$ are block-diagonal, and the
Schur complement type matrix $-\hat{B}\hat{C}\hat{B}^{\top}$ has only off-
block diagonal entries. Note that the coupling has been shifted from the port
matrices $B_{i}$ to the off-block diagonal part of $J$. $\Box$

Now we collect some properties of systems (2) with properties as in Assumption
5. To this end, we note that a matrix pencil $sE-A\in\mathbb{R}[s]^{n\times
n}$ is called regualar, if $\det(sE-A)$ is not the zero polynomial.
Furthermore, for a regular pencil $sE-A$, the index  $\nu\in\mathbb{N}_{0}$ of
the DAE $\tfrac{\mathrm{d}}{\mathrm{d}t}Ex(t)=Ax(t)+f(t)$ is the smallest
number for which the mapping $\lambda\mapsto\lambda^{1-\nu}(\nu E-A)^{-1}$ is
bounded outside a compact subset of $\mathbb{C}$. Note that the index equals
to the number of differentiations of the DAE needed to obtain an ordinary
differential equation [17, Chap. 2].

######  Proposition 8.

Any system (2) fulfilling Assumption 5 has the following properties.

  1. (a)

The pencil $sE-(J-R)Q\in\mathbb{R}[s]^{n\times n}$ is regular.

  2. (b)

There exists a unique solution $x$ of (2).

  3. (c)

The index of the DAE (2) is at most two.

###### Proof.

  1. (a)

For $\lambda\in\mathbb{C}$ with $\operatorname{Re}\lambda>0$, we show that
$C:=\lambda E-(J-R)Q$ has a trivial kernel. Let $z\in\mathbb{C}^{n}$ with
$Cz=0$. Then

$0=\operatorname{Re}(z^{*}Q^{\top}Cz)=\operatorname{Re}(\lambda
z^{*}Q^{\top}Ez)-\operatorname{Re}(z^{*}Q^{\top}JQx)+\operatorname{Re}(z^{*}Q^{\top}RQz).$

Invoking that $J=-J^{\top}$, $Q^{\top}E=E^{\top}Q\geq 0$ and $R=R^{\top}\geq
0$, we obtain

$\displaystyle\operatorname{Re}(\lambda z^{*}Q^{\top}Ez)=$
$\displaystyle\,\operatorname{Re}(\lambda)z^{*}Q^{\top}Ez\geq 0,$
$\displaystyle\operatorname{Re}(z^{*}Q^{\top}JQz)=$
$\displaystyle\,\operatorname{Re}((Qz)^{*}J(Qz))=0,$
$\displaystyle\operatorname{Re}(z^{*}Q^{\top}RQz)=$
$\displaystyle\,\operatorname{Re}((Qz)^{*}R(Qz))=(Qz)^{*}R(Qx)\geq 0,$

which leads to $z^{*}Q^{\top}Ez=0$ and $(Qz)^{*}R(Qz)=0$. The positive semi-
definiteness of $Q^{\top}E$ and $R$ implies that $E^{\top}Qz=Q^{\top}Ez=0$ and
$RQz=0$. Using $(\lambda E-(J-R)Q)z=0$, this yields $JQz=0$. Hence, by
$R=R^{\top}$ and $J=-J^{\top}$,

$Qz\in\ker\left[\begin{matrix}E&R&J\end{matrix}\right]^{\top}.$

Assumption 5 (b) now gives $Qz=0$. Since $Q$ is invertible, we are led to
$z=0$.

  2. (b)

By using that the pencil $sE-(J-R)Q$ is regular, uniqueness of solutions of
the initial value problem (2) follows from the considerations in [18, Sec.
1.1].
It remains to prove that a solution of (2) exists. Let
$Z\in\mathbb{R}^{n\times r}$, $Z_{1}\in\mathbb{R}^{r\times r_{1}}$ be as in
Assumption 5 (c). Furthermore, let $Z^{\prime}\in\mathbb{R}^{n\times(n-r)}$,
$Z_{1}^{\prime}\in\mathbb{R}^{r\times(r-r_{1})}$ be matrices with
$\operatorname{im}Z^{\prime}=(\operatorname{im}Z)^{\bot}$ and
$\operatorname{im}Z_{1}^{\prime}=(\operatorname{im}Z_{1})^{\bot}$. Moreover,
we introduce matrices $Z_{2}\in\mathbb{R}^{(n-r)\times r_{1}}$,
$Z_{2}^{\prime}\in\mathbb{R}^{(n-r)\times(n-r-r_{1})}$ with full column rank
and

$\operatorname{im}Z_{2}=\ker
Z_{1}^{\top}Z^{\top}Q^{\top}(J-R)Z^{\prime},\;\;\operatorname{im}Z_{2}^{\prime}=(\operatorname{im}Z_{2})^{\bot}.$

The construction of $Z,Z^{\prime},Z_{1},Z_{1}^{\prime},Z_{2},Z_{2}^{\prime}$
yields that the columns of the matrix

$V:=\big{[}Z^{\prime}Z_{2}^{\prime},Z^{\prime}Z_{2},ZZ_{1}^{\prime},ZZ_{1}\big{]}$
(6)

are linearly independent, and thus, $V,V^{\top}Q^{\top}\in\mathbb{R}^{n\times
n}$ are invertible matrices. Then we obtain

$V^{\top}Q^{\top}(sE-(J-R)Q)V=\begin{bmatrix}sE_{11}-A_{11}&sE_{12}-A_{12}&-A_{13}&-A_{14}\\\
sE_{12}^{\top}-A_{21}&sE_{22}-A_{22}&-A_{23}&0\\\ -A_{31}&-A_{32}&-A_{33}&0\\\
A_{14}^{\top}&0&0&0\end{bmatrix},$ (7)

where

$\displaystyle E_{11}=$
$\displaystyle\,(Z_{2}^{\prime})^{\top}(Z^{\prime})^{\top}Q^{\top}EZ^{\prime}Z_{2}^{\prime},$
$\displaystyle E_{12}=$
$\displaystyle\,(Z_{2}^{\prime})^{\top}(Z^{\prime})^{\top}Q^{\top}EZ^{\prime}Z_{2},$
$\displaystyle E_{22}=$
$\displaystyle\,Z_{2}^{\top}(Z^{\prime})^{\top}Q^{\top}EZ^{\prime}Z_{2},$
$\displaystyle A_{11}=$
$\displaystyle\,(Z_{2}^{\prime})^{\top}(Z^{\prime})^{\top}Q^{\top}(J-R)QZ^{\prime}Z_{2}^{\prime},$
$\displaystyle A_{12}=$
$\displaystyle\,(Z_{2}^{\prime})^{\top}(Z^{\prime})^{\top}Q^{\top}(J-R)QZ^{\prime}Z_{2},$
$\displaystyle A_{21}=$
$\displaystyle\,Z_{2}^{\top}(Z^{\prime})^{\top}Q^{\top}(J-R)QZ^{\prime}Z_{2}^{\prime},$
$\displaystyle A_{13}=$
$\displaystyle\,(Z_{2}^{\prime})^{\top}(Z^{\prime})^{\top}Q^{\top}(J-R)QZZ_{1}^{\prime},$
$\displaystyle A_{31}=$
$\displaystyle\,(Z_{1})^{\prime\top}Z^{\top}Q^{\top}(J-R)QZ^{\prime}Z_{2}^{\prime},$
$\displaystyle A_{22}=$
$\displaystyle\,Z_{2}^{\top}(Z^{\prime})^{\top}Q^{\top}(J-R)QZ^{\prime}Z_{2},$
$\displaystyle A_{33}=$
$\displaystyle\,(Z_{1})^{\prime\top}Z^{\top}Q^{\top}(J-R)QZZ_{1}^{\prime},$
$\displaystyle A_{14}=$
$\displaystyle\,(Z_{2}^{\prime})^{\top}(Z^{\prime})^{\top}Q^{\top}JQZZ_{1},$
$\displaystyle B_{1}=$
$\displaystyle\,(Z_{2}^{\prime})^{\top}(Z^{\prime})^{\top}Q^{\top}B,$
$\displaystyle B_{2}=$
$\displaystyle\,Z_{2}^{\top}(Z^{\prime})^{\top}Q^{\top}B,$ $\displaystyle
B_{3}=$ $\displaystyle\,(Z_{1}^{\prime})^{\top}Z^{\top}Q^{\top}B,$
$\displaystyle B_{4}=$ $\displaystyle\,Z_{1}^{\top}Z^{\top}Q^{\top}B.$

The construction of the matrices
$Z,Z^{\prime},Z_{1},Z_{1}^{\prime},Z_{2},Z_{2}^{\prime}$ yields that

     * •

$\left[\begin{smallmatrix}E_{11}&E_{12}\\\
E_{12}^{\top}&E_{22}\end{smallmatrix}\right]$ is positive definite. In
particular, $E_{22}$ is invertible,

     * •

$A_{33}$ is invertible, and

     * •

$A_{14}$ has full row rank.

The already proved regularity of $sE-(J-R)Q$, which implies regularity of the
pencil $V^{\top}Q^{\top}(sE-(J-R)Q)QV$, which yields that $A_{14}$ has
moreover full column rank. As a consequence, $A_{14}$ is square (i.e.,
$2r_{1}=n-r$) and invertible.
Set $z_{0}=(z_{10}^{\top},z_{20}^{\top},z_{30}^{\top},z_{40}^{\top})$,
$z_{10}\in\mathbb{R}^{n-r-r_{1}}$, $z_{20}\in\mathbb{R}^{r_{1}}$,
$z_{30}\in\mathbb{R}^{r-r_{1}}$, $z_{40}\in\mathbb{R}^{r_{1}}$. By Assumption
5 (c), we have $B_{4}u\in H^{1}([0,T];\mathbb{R}^{r_{1}})$ and

$0=Z_{1}^{\top}Z^{\top}Q^{\top}(JQx_{0}+Bu(0))=Z_{1}^{\top}Z^{\top}Q^{\top}JQx_{0}+B_{4}u(0)\\\
=Z_{1}^{\top}Z^{\top}Q^{\top}JQ(Z^{\prime}Z_{2}^{\prime}z_{10}+Z^{\prime}Z_{2}z_{20}+ZZ_{1}^{\prime}z_{30}+ZZ_{1}z_{40})+B_{4}u(0)\\\
=Z_{1}^{\top}Z^{\top}Q^{\top}JQZ^{\prime}Z_{2}^{\prime}z_{10}+B_{4}u(0)=-Z_{1}^{\top}Z^{\top}Q^{\top}J^{\top}QZ^{\prime}Z_{2}^{\prime}z_{10}+B_{4}u(0)\\\
=-A_{14}^{\top}z_{10}+B_{4}u(0).$

Consequently, $z_{1}=(A_{14}^{\top})^{-1}B_{4}u$ fulfills $z_{1}\in
H^{1}([0,T];\mathbb{R}^{r})$, $-A_{14}^{\top}z_{1}+B_{4}u(t)=0$ and
$z_{1}(0)=z_{10}$. Moreover, let $z_{2}\in H^{1}([0,T];\mathbb{R}^{r_{1}})$ be
the solution of the following ODE with initial value $z_{2}(0)=z_{20}$

$\tfrac{\mathrm{d}}{\mathrm{d}t}z_{2}(t)=E_{22}^{-1}(A_{22}-A_{23}A_{33}^{-1}A_{32})z_{2}(t)\\\
\quad+E_{22}^{-1}\big{(}-\tfrac{\mathrm{d}}{\mathrm{d}t}E_{12}^{\top}z_{1}(t)+A_{21}z_{1}(t)+B_{2}u(t)-A_{23}A_{33}^{-1}z_{1}(t)-A_{23}A_{33}^{-1}B_{3}u(t)\big{)},$

and we successively define $z_{3}\in L^{2}([0,T];\mathbb{R}^{r-r_{1}})$,
$z_{4}\in L^{2}([0,T];\mathbb{R}^{r_{1}})$ by

$\displaystyle z_{3}(t)=$
$\displaystyle\,-A_{33}^{-1}A_{23}z_{2}(t)-A_{33}^{-1}A_{13}z_{1}(t),$
$\displaystyle z_{4}(t)=$
$\displaystyle\,-A_{14}^{-1}\big{(}\tfrac{\mathrm{d}}{\mathrm{d}t}E_{11}z_{1}(t)+\tfrac{\mathrm{d}}{\mathrm{d}t}E_{12}z_{2}(t)-A_{11}z_{1}(t)-A_{12}z_{1}(t)-A_{13}z_{3}(t)\big{)}.$

Altogether, we have that
$z:=(z_{1}^{\top},z_{2}^{\top},z_{3}^{\top},z_{4}^{\top})^{\top}$ is a
solution of the differential-algebraic equation

$\displaystyle\tfrac{\mathrm{d}}{\mathrm{d}t}\begin{bmatrix}E_{11}&E_{12}&0&0\\\
E_{12}^{\top}&E_{22}&0&0\\\ 0&0&0&0\\\
0&0&0&0\end{bmatrix}\begin{pmatrix}z_{1}(t)\\\ z_{2}(t)\\\ z_{3}(t)\\\
z_{4}(t)\end{pmatrix}=$
$\displaystyle\begin{bmatrix}A_{11}&A_{12}&A_{13}&A_{14}\\\
A_{21}&A_{22}&A_{23}&0\\\ A_{31}&A_{32}&A_{33}&0\\\
-A_{14}^{\top}&0&0&0\end{bmatrix}\begin{pmatrix}z_{1}(t)\\\ z_{2}(t)\\\
z_{3}(t)\\\ z_{4}(t)\end{pmatrix}+\begin{bmatrix}B_{1}\\\ B_{2}\\\ B_{3}\\\
B_{4}\end{bmatrix}u(t),$ (8) $\displaystyle\begin{bmatrix}E_{11}&E_{12}&0&0\\\
E_{12}^{\top}&E_{22}&0&0\\\ 0&0&0&0\\\
0&0&0&0\end{bmatrix}\begin{pmatrix}z_{1}(0)\\\ z_{2}(0)\\\ z_{3}(0)\\\
z_{4}(0)\end{pmatrix}=$ $\displaystyle\begin{bmatrix}E_{11}&E_{12}&0&0\\\
E_{12}^{\top}&E_{22}&0&0\\\ 0&0&0&0\\\
0&0&0&0\end{bmatrix}\begin{pmatrix}z_{10}\\\ z_{20}\\\ z_{30}\\\
z_{40}\end{pmatrix}.$

Now setting

$x(t):=(Q^{\top})^{-1}(V^{\top})^{-1}z(t),$

we obtain from (7) and (8) that $x$ is a solution of (2).

  3. (c)

The result has been shown in [22, Thm. 4.3] and [12, Thm. 6.6]. Notice, this
can be also seen from the proof of (b). ∎

######  Remark 9.

Assume that the DAE (2) fulfills Assumption 5 (a) & (b). By using
$V\in\mathbb{R}^{n\times n}$ as defined as in (6), we can transform (2) to an
equivalent DAE (8). Consequently, the existence of a solution of (2) implies
$z_{1}\in H^{1}([0,T];\mathbb{R}^{r_{1}})$ with, and thus
$B_{4}u=(A_{14}^{\top})^{-1}\in H^{1}([0,T];\mathbb{R}^{r_{1}})$. Invoking
that $B_{4}=Z_{1}^{\top}Z^{\top}Q^{\top}B$, we obtain that
$Z_{1}^{\top}Z^{\top}Q^{\top}Bu\in H^{1}([0,T];\mathbb{R}^{r_{1}})$ with
$Z_{1}^{\top}Z^{\top}Q^{\top}(JQx_{0}+Bu(0))=0$ has to be fulfilled.
Consequently, under Assumption 5 (a) & (b), the conditions in Assumption 5 (c)
are also necessary for the existence of a solution of the DAE (2). $\Box$

###  2.2 Perspective as coupled PH-subsystems

We take the perspective of the paper [14] for the system (2) and its
partitioning given in Rem. 6. To this end, we need to introduce the internal
inputs $u_{ij}$, i.e., data stemming from the $j$th system being input for the
$i$th system, and outputs $y_{ij}$ (i.e., data from the $i$ system to be
transferred to the $j$th system). This read for the $i$th subsystems
($i=1,\ldots,s$):

$\displaystyle\tfrac{\mathrm{d}}{\mathrm{d}t}E_{i}x_{i}(t)$
$\displaystyle=(J_{i}-R_{i})Q_{i}x_{i}(t)+\sum_{\begin{smallmatrix}j=1\\\
j\neq i\end{smallmatrix}}^{s}B_{ij}{u}_{ij}(t)+B_{i}u_{i}(t),\qquad
E_{i}x_{i}(0)=Ex_{i0},$ $\displaystyle y_{ij}(t)$
$\displaystyle=B_{ij}^{\top}Q_{i}x_{i}(t),\qquad j=1,\ldots,s,\,j\neq i,$
$\displaystyle y_{i}(t)$ $\displaystyle=B_{i}^{\top}Q_{i}x_{i}(t)$

with matrices $E_{i},Q_{i},J_{i},R_{i}\in\mathbb{R}^{n_{i}\times n_{i}}$, and
$B_{i}\in\mathbb{R}^{n_{i}\times m}$ fulfilling the properties of Assumption 5
(a), and $B_{ij}\in\mathbb{R}^{n_{i}\times m_{ij}}$ for $i=1,\ldots,s$,
$j=1,\ldots,s$, $i\neq j$ with $m_{ij}=m_{ji}$. The dissipation inequality for
this subsystem reads

$\tfrac{1}{2}\big{(}x_{i}(t)^{\top}Q_{i}^{\top}E_{i}x_{i}(t)\big{)}-\tfrac{1}{2}\big{(}x_{i}(0)^{\top}Q_{i}^{\top}E_{i}x_{i}(0)\big{)}\\\
\leq\int_{0}^{t}u_{i}(\tau)^{\top}y_{i}(\tau)d\tau+\sum_{\begin{smallmatrix}j=1\\\
j\neq
i\end{smallmatrix}}^{s}\int_{0}^{t}u_{ij}(\tau)^{\top}y_{ij}(\tau)d\tau\qquad\forall\,t\in[0,T]$

The overall input is $u=u_{1}=\ldots=u_{s}$, whereas the overall output is
given by the sum

$y=\sum_{j=1}^{s}y_{i}.$

Coupling is done such that $x\mapsto x^{\top}Q^{\top}Ex$ with $E$ and $Q$ as
in (3) is a storage function for the overall system [11]. This is for instance
achieved by

$\displaystyle\forall\,i=1,\ldots,s:$ $\displaystyle u_{ij}+y_{ji}=0,\qquad
j=1,\ldots,i-1,$ $\displaystyle u_{ji}-y_{ij}=0,\qquad j=i+1,\ldots,s,$

since

$\displaystyle\,\tfrac{1}{2}\big{(}x(t)^{\top}Q^{\top}Ex(t)\big{)}-\tfrac{1}{2}\big{(}x(0)^{\top}Q^{\top}Ex(0)\big{)}$
$\displaystyle=$
$\displaystyle\,\sum_{i=1}^{s}\tfrac{1}{2}\Big{(}\big{(}x_{i}(t)^{\top}Q_{i}^{\top}E_{i}x_{i}(t)\big{)}-\tfrac{1}{2}\big{(}x_{i}(0)^{\top}Q_{i}^{\top}E_{i}x_{i}(0)\big{)}\Big{)}$
$\displaystyle\leq$
$\displaystyle\sum_{i=1}^{s}\int_{0}^{t}u_{i}(\tau)^{\top}y_{i}(\tau)d\tau+\sum_{i=1}^{s}\sum_{\begin{smallmatrix}j=1\\\
j\neq i\end{smallmatrix}}^{s}\int_{0}^{t}u_{ij}(\tau)^{\top}y_{ij}(\tau)d\tau$
$\displaystyle=$
$\displaystyle\int_{0}^{t}\\!\\!u(\tau)^{\\!\top}\\!\sum_{i=1}^{s}y_{i}(\tau)d\tau\\!+\\!\sum_{i=1}^{s}\sum_{j=1}^{i-1}\\!\int_{0}^{t}\\!\\!u_{ij}(\tau)^{\\!\top}\\!y_{ij}(\tau)d\tau\\!+\\!\sum_{i=1}^{s}\sum_{j=i+1}^{s}\\!\int_{0}^{t}\\!\\!u_{ij}(\tau)^{\\!\top}y_{ij}(\tau)d\tau$
$\displaystyle=$
$\displaystyle\int_{0}^{t}\\!\\!u(\tau)^{\\!\top}\\!y(\tau)d\tau+\sum_{i=1}^{s}\sum_{j=1}^{i-1}\\!\int_{0}^{t}\\!\\!\\!-y_{ji}(\tau)^{\\!\top}u_{ji}(\tau)d\tau+\sum_{i=1}^{s}\sum_{j=i+1}^{s}\\!\int_{0}^{t}\\!\\!u_{ij}(\tau)^{\\!\top}\\!y_{ij}(\tau)d\tau$
$\displaystyle=$ $\displaystyle\,\int_{0}^{t}u(\tau)^{\top}y(\tau)d\tau.$

The overall system is given by (2) with $E$, $Q$ and $R$ as in (3), and $J$
structured as in (4) with

$\displaystyle\forall\,i=1,\ldots,s-1,\;j=i+1,\ldots,s:$ $\displaystyle\quad
J_{ij}=B_{ij}B_{ji}^{\top}.$

######  Remark 10.

This coupling has the drawback that all variables enter the definition of the
coupling. And we have as many outputs of a subsystem as the other subsystem
has variables times the number of other subsystems. This will be treated
differently in practice. A ’sparse’ coupling has to be introduced.

Before we investigate the Lions/Mercier based algorithm, we revisit the Jacobi
dynamic iteration and discuss some limitations.

##  3 Failure of Jacobi dynamic iteration

To motivate the development of monotone operator splitting based dynamic
iteration schemes in chapter 5, we remind that classical approaches such as
the Jacobi or Gauss-Seidel iteration scheme may not always be feasible for a
computer implementation. To this end, we address system (2) in the following
ODE setting: assigning $E=Q=I$, the initial value problem (2) reads

$\displaystyle\dot{x}$ $\displaystyle=(J-R)x+Bu(t),\quad x(0)=x_{0},\quad
t\in[0,T].$ (9)

For this system, we recall that the dynamic iteration is based on a splitting
$J-R=M-N$ and an iteration count $k\in\mathbb{N}_{0}$. Now, let be an initial
guess $k=0$ be given $x^{(0)}\in C^{0}([0,T],{}^{n})$ with $x^{(0)}(0)=x_{0}$
(e.g., $x^{(0)}(t)=x_{0}$) for all $t\in[0,T]$, the dynamic iteration scheme
reads:

$\displaystyle\dot{x}^{(k+1)}$ $\displaystyle=Mx^{(k+1)}-Nx^{(k)}+Bu(t),\quad
x^{(k+1)}(0)=x_{0}.$ (10)

######  Remark 11.

Note that setting $M=J_{d}-R$ and $N=-J_{o}$ in the case of the pHS-ODE system
(9) defines a block-Jacobi dynamic iteration scheme. $\Box$

If we restrict to $u\in C^{0}([0,T],{}^{m})$, we obtain from (10)
$x^{(k+1)}\in C^{0}([0,T],{}^{n})$. Let $x\in C^{0}([0,T],{}^{n})$ denote the
analytic solution of (9). Then, the recursion error is given by
$\epsilon^{(k)}:=x^{(k)}-x$. Following Burrage [10], we can derive an exact
error recursion for the dynamic iteration (10). Doing this, for the system (9)
with

$R=\operatorname{diag}(\tau,\dotsc,\tau),\quad B=0,$ (11)

and setting $M=-R,N=-J_{o}=-J$, we get the error recursion on $t\in[0,\,T]$:

$\displaystyle\epsilon^{(k)}\\!(t)=\int_{0}^{t}\\!\\!e^{-(t-s_{k\\!-\\!1})R}J\epsilon^{(k-1)}(s_{k\\!-\\!1})\text{d}s_{k-1}=J\\!\\!\int_{0}^{t}\\!\\!e^{-\tau(t-s_{k\\!-\\!1})}\epsilon^{(k-1)}(s_{k-1})\text{d}s_{k-1}$
$\displaystyle=J^{k}\\!\\!\int_{0}^{t}\\!\\!e^{-\tau(t-s_{k-1})}\\!\\!\int_{0}^{s_{k\\!-\\!1}}\\!\\!\\!e^{-\tau(s_{k\\!-\\!1}\\!-\\!s_{k\\!-\\!2})}\\!\\!\dotsc\\!\int_{0}^{s_{1}}\\!\\!\\!e^{-\tau(s_{1}\\!-\\!s_{0})}\epsilon^{(0)}(s_{0})\text{d}s_{0}\dotsc\text{d}s_{k-2}\text{d}s_{k-1}$
$\displaystyle=e^{-\tau
t}J^{k}\int_{0}^{t}\int_{0}^{s_{k-1}}\dotsi\int_{0}^{s_{1}}e^{\tau
s_{0}}\epsilon^{(0)}(s_{0})\text{d}s_{0}\dotsi\text{d}s_{k-2}\text{d}s_{k-1}.$

Now, just for demonstration purposes and simplicity, let us assume that
initial error has the format:

$\epsilon^{(0)}(s_{0})=Ds_{0}e^{-\tau s_{0}}\mathbbm{1}\;\text{ with
}\mathbbm{1}=[1,\,1]^{\top}$

and $D>0$. Using the maximum norm $\|f\|_{T}:=\max_{t\in[0,T]}\|f(t)\|_{2}$
for $f\in C^{0}([0,T],{}^{n})$, the error reads for all $k\geq\tau T-1$
(maximal error at $T$):

$\|\epsilon^{(k)}\|_{T}=\sqrt{2}\,D\|J\|_{2}^{k}\frac{e^{-\tau
T}T^{k+1}}{(k+1)!}\quad\Rightarrow\quad\frac{\|\epsilon^{(k+1)}\|_{T}}{\|\epsilon^{(k)}\|_{T}}=\frac{\|J\|_{2}\cdot
T}{k+2}.$

######  Example 12.

We analyse a $2\times 2$ version of (9) with the restriction (11), which
reads:

$\displaystyle x^{\prime}=(J-R)x,\quad x(0)=x_{0}=\begin{bmatrix}2\\\
2\end{bmatrix},\qquad\text{with }J=\begin{bmatrix}0&\nu\\\
-\nu&0\end{bmatrix}\\!,\;R=\begin{bmatrix}\tau&0\\\ 0&\tau\end{bmatrix}$ (12)

and employing $\tau,\nu>0$. The analytic solution and error evolution are
given by:

$x(t)=e^{-\tau t}\begin{bmatrix}\cos(\nu t)&\sin(\nu t)\\\ -\sin(\nu
t)&\cos(\nu
t)\end{bmatrix}x_{0},\qquad\frac{\|\epsilon^{(k+1)}\|_{T}}{\|\epsilon^{(k)}\|_{T}}=\frac{\nu\,T}{k+2}\quad\forall
k\geq\tau T-1.$

The error amplification factor $\nu T/(k+2)$ is bounded by one for all $k>\nu
T-2$. For all $k<\nu T-2$ the error is increasing. Hence, the convergence is
monotone (for all $k$) only if $\nu T-2\leq 0$. We give results for both
cases:

  1. (a)

For $\nu T>2$, we obtain a decrease of the error only if $k$ is large enough.
On the other hand, we may trigger overflows if $\nu T$ is large enough, see
Fig. 1, left.

  2. (b)

If $\nu T\leq 2$, we have monotone convergence. Hence windowing is always a
solution to obtain (good) convergence.

An example for $T=0.5$ is shown in Fig. 1, right. Thus using Jacobi with
windows, we can obtain convergences. Note that the behavior in the window
$[0,0.5]$ yields a kind of worst case for the convergence. $\Box$

Figure 1: Jacobi iteration for system (12) with $\nu=15$, $\tau=0.01$, $D=1$.
Left: $T=10$. The maximum error would be obtained at iteration $k=148$:
$\sqrt{2}\cdot 150^{148}/148!$ and error decrease from $k=149$ on. Overflow in
single precision at $k=40$.
Right: $T=0.5$. Maximum error for $k=5$ ($\nu T-2=5.5$), then decrease. After
$k=29$, we achieve the solution (in single precision).

##  4 Recap on maximal monotone operators

For $\omega\geq 0$, we introduce the weighted $L^{2}$-space

$L^{2}_{\omega}([0,T],\mathbb{R}^{n}):=\\{f:[0,T]\rightarrow\mathbb{R}^{n}\mid
t\mapsto e^{-t\omega}f(t)\in L^{2}([0,T],\mathbb{R}^{n})\\}$

equipped with the norm $\|\cdot\|_{2,\omega}$ defined as

$\|f\|_{2,\omega}^{2}=\int_{0}^{T}e^{-2t\omega}\|f(t)\|^{2}dt.$

Clearly, as $T\in(0,\infty)$, we have
$L^{2}_{\omega}([0,T],\mathbb{R}^{n})=L^{2}([0,T],\mathbb{R}^{n})$ with
equivalent norms. Furthermore, if $\omega=0$, then the norms in
$L^{2}_{\omega}([0,T],\mathbb{R}^{n})$ and $L^{2}([0,T],\mathbb{R}^{n})$
coincide, i.e.,

$\forall\,f\in
L^{2}([0,T],\mathbb{R}^{n}):\quad\|f\|_{2,0}=\|f\|_{2}=\left(\int_{0}^{T}\|f(t)\|^{2}dt\right)^{1/2}.$

######  Definition 13 (Contraction).

Let $X$ be a Banach space space with norm $\|\cdot\|$. A (possibly nonlinear)
operator $A:D(A)\subset X\to X$ is called contractive, if

$\|Ax-Ay\|\leq\|x-y\|,\qquad\forall\,x,y\in X,$

and strictly contractive, if there exists some $L<1$, such that

$\|Ax-Ay\|\leq L\,\|x-y\|,\qquad\forall\,x,y\in X.$ $\Box$

######  Definition 14 (Maximally monotone operator).

Let $X$ be a real Hilbert space with inner product
$\langle\cdot,\cdot\rangle$. A set $Y\subset X\times X$ is called monotone, if

$\langle x-u,y-v\rangle\geq 0,\qquad(x,y),(u,v)\in Y.$

Further, $Y\subset X\times X$ is called maximally monotone, if it is monotone
and not a proper subset of a monotone subset of $X\times X$.
A (possibly nonlinear) operator $A:D(A)\subset X\rightarrow X$ is called
(maximally) monotone, if the graph of $A$, i.e., $\\{(x,Ax):\,x\in D(A)\\}$,
is (maximally) monotone. $\Box$

######  Remark 15.

Let $A:D(A)\subset X\rightarrow X$ be a monotone operator. It follows from the
definition of monotonicity that $I+\lambda A$ is injective for all
$\lambda>0$.
Moreover, by [5, Theorem 2.2], the following three statements are equivalent:

  1. (i)

$A$ is maximally monotone,

  2. (ii)

$I+\lambda A$ is surjective for some $\lambda>0$,

  3. (iii)

$I+\lambda A$ is surjective for all $\lambda>0$.

Consequently, if $A$ is maximally monotone, then $I+\lambda A$ is bijective
for all $\lambda>0$. The Cauchy-Schwarz inequality yields that

$\|(I+\lambda A)x-(I+\lambda A)y\|\geq\|x-y\|\quad\forall x,y\in D(A),$

whence $(I+\lambda A)^{-1}:X\rightarrow X$ is contractive.

Furthermore, $(I-\lambda A)(I+\lambda A)^{-1}$ is contractive. This follows
with $\tilde{x}:=(I+\lambda A)^{-1}x$ and $\tilde{y}:=(I+\lambda A)^{-1}y$
from

$\displaystyle\|x-y\|^{2}-\|(I-\lambda A)$ $\displaystyle(I+\lambda
A)^{-1}x-(I-\lambda A)(I+\lambda A)^{-1}y\|^{2}$
$\displaystyle=\|\tilde{x}-\tilde{y}+\lambda A\tilde{x}-\lambda
A\tilde{y}\|^{2}-\|\tilde{x}-\tilde{y}-(\lambda A\tilde{x}-\lambda
A\tilde{y})\|^{2}$
$\displaystyle=4\lambda\langle\tilde{x}-\tilde{y},A\tilde{x}-A\tilde{y}\rangle\geq
0.$ $\Box$

##  5 Dynamic iteration scheme

Now we develop a dynamic iteration scheme for the DAE (2) fulfilling
Assumption 5. To facilitate the decoupling, we split $J$ as follows

$J=\,\underbrace{\begin{bmatrix}J_{1}&0&\cdots&0\\\ 0&\ddots&\ddots&\vdots\\\
\vdots&\ddots&\ddots&0\\\
0&\cdots&0&J_{s}\end{bmatrix}}_{=:J_{d}}+\,\underbrace{\begin{bmatrix}0&J_{12}&\cdots&J_{1s}\\\
-J_{12}^{\top}&\ddots&\ddots&\vdots\\\ \vdots&\ddots&\ddots&J_{s-1,s}\\\
-J_{1s}^{\top}&\cdots&-J_{s-1,s}^{\top}&0\end{bmatrix}}_{=:J_{o}}.$ (13)

such that $J_{d}=-J_{d}^{\top}$ and $J_{o}=-J_{o}^{\top}$ (compare with Rem. 6
c).

Let $\alpha\in[0,1]$, $\mu,\omega\in\mathbb{R}$ with $0\leq\mu\leq\omega$. We
introduce the affine-linear operator $M$ by:

$\displaystyle M\\!:$ $\displaystyle\,D(M)\subset
L^{2}_{\omega}([0,T],\mathbb{R}^{n})$ $\displaystyle\rightarrow
L^{2}_{\omega}([0,T],\mathbb{R}^{n}),$ (14a) $\displaystyle x$
$\displaystyle\mapsto\tfrac{\mathrm{d}}{\mathrm{d}t}(E{x})-(J_{d}-\alpha R+\mu
EQ^{-1})Qx\\!-\\!Bu$ (14b) with domain $D(M)=\left\\{x\in
L^{2}_{\omega}([0,T],\mathbb{R}^{n})\mid Ex\in
H^{1}([0,T],\mathbb{R}^{n}),Ex(0)=Ex_{0}\right\\},$ (14c)

and the linear bounded multiplication operator $N$ by

$\displaystyle N\\!:$ $\displaystyle L^{2}_{\omega}([0,T],\mathbb{R}^{n})$
$\displaystyle\,\rightarrow L^{2}_{\omega}([0,T],\mathbb{R}^{n}),$ (15a)
$\displaystyle x$ $\displaystyle\,\mapsto((1-\alpha)R+\mu EQ^{-1}-J_{o})Qx.$
(15b)

######  Remark 16.

Notice that $x$ is a solution of the DAE (2) if, and only if,

$Mx+Nx=0$ (16)

holds.  $\Box$

######  Lemma 17.

Let $\omega,\mu\in\mathbb{R}$ with $0\leq\mu\leq\omega$, and let
$\alpha\in[0,1]$. Assume that the DAE (2) fulfills Assumption 5, and let the
operator $M$ be defined as in (14). Then $MQ^{-1}$ is maximally monotone on
$L^{2}_{\omega}([0,T],\mathbb{R}^{n})$. More precisely, $MQ^{-1}$ satisfies

$\displaystyle\langle MQ^{-1}v$ $\displaystyle-
MQ^{-1}w,\;v-w\rangle_{2,\omega}$
$\displaystyle=\tfrac{1}{2}e^{-2T\omega}\left(\|(EQ^{-1})^{1/2}(v(T)-w(T))\|^{2}\right)$
$\displaystyle\qquad(\omega-\mu)\|(EQ^{-1})^{1/2}(v-w)\|^{2}_{2,\omega}+\alpha\|R^{1/2}(v-w)\|^{2}_{2,\omega}\geq
0$

for $v,w\in QD(M)$, and $I+\lambda MQ^{-1}$ is surjective for every
$\lambda>0$.

###### Proof.

First of all, we have for $v,w\in QD(M)$

$\displaystyle\langle MQ^{-1}v-MQ^{-1}w,v-w\rangle_{2,\omega}$
$\displaystyle=\langle\tfrac{\mathrm{d}}{\mathrm{d}t}EQ^{-1}(v-w),\,v-w\rangle_{2,\omega}-\langle(J_{d}-\alpha
R+\mu EQ^{-1})(v-w),v-w\rangle_{2,\omega}$
$\displaystyle=\langle\tfrac{\mathrm{d}}{\mathrm{d}t}EQ^{-1}(v\\!-\\!w),v-w\rangle_{2,\omega}-\mu\langle(EQ^{-1}(v\\!-\\!w),v\\!-\\!w\rangle_{2,\omega}+\langle\alpha
R(v\\!-\\!w),v\\!-\\!w\rangle_{2,\omega}$
$\displaystyle=\tfrac{1}{2}e^{-2T\omega}\left(\|(EQ^{-1})^{1/2}(v(T)-w(T))\|^{2}\right)$
$\displaystyle\qquad+(\omega-\mu)\|(EQ^{-1})^{1/2}(v-w)\|^{2}_{2,\omega}+\langle\alpha
R(v-w),v-w\rangle_{2,\omega}\geq 0,$

which shows that $MQ^{-1}$ is monotone.
To prove that $MQ^{-1}$ is maximally monotone, it suffices, by Remark 15, to
prove that $I+\lambda MQ^{-1}$ is surjective for some $\lambda>0$. Assume that
$y\in L^{2}_{\omega}([0,T],\mathbb{R}^{n})$. Consider the matrices
$\tilde{E}=EQ^{-1}$, $\tilde{R}=\alpha R+\frac{1}{\lambda}I$ and
$\tilde{Q}=I$. Then the DAE

$\tfrac{\mathrm{d}}{\mathrm{d}t}\tilde{E}\tilde{v}(t)=(J_{d}-\tilde{R})\tilde{Q}\tilde{v}(t)+Be^{-\mu
t}u(t)-\tfrac{1}{\lambda}e^{-\mu
t}y(t),\quad\tilde{E}\tilde{v}(0)=\tilde{E}Qx_{0},$ (17)

This DAE clearly fulfills Assumption 5 (a)&(b). Moreover, by positive
definiteness of $\tilde{R}$, we obtain that any matrix
$Z\in\mathbb{R}^{n\times r}$ with full column rank and
$\operatorname{im}Z=\ker E$ fulfills $\ker\tilde{R}QZ=\\{0\\}$. This means
that Assumption 5 (c) is trivially fulfilled by (17), and Proposition 8
implies that (17) has a unique solution $\tilde{v}\in
L^{2}_{\omega}([0,T],\mathbb{R}^{n})$ with $\tilde{E}\tilde{v}\in
H^{1}([0,T],\mathbb{R}^{n})$. Now consider ${v}\in
L^{2}_{\omega}([0,T],\mathbb{R}^{n})$ with $v(t)=e^{\mu t}\tilde{v}(t)$. Then
$\tilde{E}v(0)=\tilde{E}\tilde{v}(0)=\tilde{E}Qx_{0}$, and the product rule
for weak derivatives [2, Thm. 4.25] yields $\tilde{E}{v}\in
H^{1}([0,T],\mathbb{R}^{n})$ with

$\displaystyle\tfrac{\mathrm{d}}{\mathrm{d}t}\tilde{E}v=$
$\displaystyle\,\tfrac{\mathrm{d}}{\mathrm{d}t}e^{\mu t}\tilde{E}\tilde{v}$
$\displaystyle=$ $\displaystyle\,e^{\mu
t}\tfrac{\mathrm{d}}{\mathrm{d}t}\tilde{E}\tilde{v}+\mu e^{\mu
t}\tilde{E}\tilde{v}$ $\displaystyle=$ $\displaystyle\,e^{\mu
t}(J_{d}-\tilde{R})\tilde{Q}\tilde{v}+e^{\mu t}Be^{-\mu t}u(t)-e^{\mu
t}\tfrac{1}{\lambda}e^{-\mu t}y(t)+\mu e^{\mu t}\tilde{E}\tilde{v}$
$\displaystyle=$
$\displaystyle\,(J_{d}-\tilde{R})\tilde{Q}{v}+Bu(t)-\tfrac{1}{\lambda}y(t)+\mu\tilde{E}{v},$

and thus

$v+\lambda\big{(}\tfrac{\mathrm{d}}{\mathrm{d}t}{E}Q^{-1}v-(J_{d}-\alpha R+\mu
EQ^{-1})QQ^{-1}v+Bu)=y.$

Since this means that $(I+\lambda MQ^{-1})v=y$, the result is shown. ∎

######  Lemma 18.

Let $\lambda>0$, $\mu\geq 0$, $\alpha\in[0,1]$, and let
$R,E,Q\in\mathbb{R}^{n\times n}$ with the properties as in Assumption 5. Then,
for $K:=((1-\alpha)R+\mu EQ^{-1}-J_{o})$

$\|(I-\lambda K)(I+\lambda K)^{-1}x\|^{2}-\|x\|^{2}\\\
=-4{\lambda}\|((1-\alpha)R+\mu EQ^{-1})^{1/2}(I+\lambda
K)^{-1}x\|^{2},\quad\forall\,x\in\mathbb{R}^{n}.$ (18)

Moreover, if

$\mathrm{rk}\begin{bmatrix}\mu E&(1-\alpha)R\end{bmatrix}=n,$ (19)

then $(1-\alpha)R+\mu EQ^{-1}=K+J_{o}\in\mathbb{R}^{n\times n}$ is invertible,
and

$\|(I-\lambda K)(I+\lambda K)^{-1}x\|\leq
q\,\|x\|\quad\forall\,x\in\mathbb{R}^{n},$ (20)

where $q\geq 0$ with

$q^{2}=1-\frac{4\lambda}{(1+\lambda\|K\|)^{2}\cdot{\|(K+J_{o})^{-1}\|}}<1.$
(21)

###### Proof.

By using that $K+K^{\top}\geq 0$, we see that $I+\lambda K$ is an invertible
matrix. Let $x\in\mathbb{R}^{n}$ and $y:=(I+\lambda K)^{-1}x$. Then

$\|x\|^{2}-\|(I-\lambda K)(I+\lambda K)^{-1}x\|^{2}\\\ =\|y+\lambda
Ky\|^{2}-\|y-\lambda Ky\|^{2}=4\lambda\langle y,(K+K^{\top})y\rangle\\\
=4\lambda\langle y,\big{(}(1-\alpha)R+\mu EQ^{-1}\big{)}y\rangle$

shows (18). To prove the remaining result, assume that (19) is fulfilled. Then
it can be concluded from Assumption 5 (a) (i) that $(1-\alpha)R+\mu EQ^{-1}$
is positive definite. Furthermore, by

$\|x\|=\|(I+\lambda
K)y\|\leq(1+\lambda\|K\|)\|y\|=(1+\lambda\|K\|)\|(I+\lambda K)^{-1}x\|,$

we have

$\|(I+\lambda
K)^{-1}x\|\geq\frac{1}{1+\lambda\|K\|}\,\|x\|\quad\forall\,x\in\mathbb{R}^{n}.$

Plugging this into (18), we obtain for all $x\in\mathbb{R}^{n}$

$\displaystyle\|(I-\lambda K)(I+\lambda
K)^{-1}x\|^{2}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!$
$\displaystyle=$ $\displaystyle\,\|x\|^{2}-4{\lambda}\|((1-\alpha)R+\mu
EQ^{-1})^{1/2}(I+\lambda K)^{-1}x\|^{2}$ $\displaystyle\leq$
$\displaystyle\,\|x\|^{2}-4\lambda\|((1-\alpha)R+\mu
EQ^{-1})^{-1}\|^{-1}\|(I+\lambda K)^{-1}x\|^{2}$ $\displaystyle\leq$
$\displaystyle\,\|x\|^{2}-4\lambda\|((1-\alpha)R+\mu
EQ^{-1})^{-1}\|^{-1}\frac{1}{(1+\lambda\|K\|)^{2}}\,\|x\|^{2},$
$\displaystyle=$ $\displaystyle\,\left(1-\frac{4\lambda}{\|((1-\alpha)R+\mu
EQ^{-1})^{-1}\|\cdot(1+\lambda\|K\|)^{2}}\right)\,\|x\|^{2},$

which completes the proof. ∎

######  Remark 19.

Condition (19) is equivalent to at least one of the following three statements
being fulfilled:

  1. (i)

$\mathrm{rk}\begin{bmatrix}E&R\end{bmatrix}=n$, $\mu>0$ and $\alpha<1$,

  2. (ii)

$\mathrm{rk}E=n$ and $\mu>0$, or

  3. (iii)

$\mathrm{rk}R=n$ and $\alpha<1$ $\Box$.

As a direct conclusion of Lemma 18, we can draw the following result.

######  Lemma 20.

Let $\lambda,\mu>0$, $\alpha\in[0,1]$, let $R,E,Q\in\mathbb{R}^{n\times n}$
with the properties as in Assumption 5, and let the operator $N$ be defined as
in (15). The multiplication operator $NQ^{-1}$ on
$L^{2}_{\omega}([0,T],\mathbb{R}^{n})$ is linear, bounded and satisfies

$\langle x,NQ^{-1}x\rangle_{2,\omega}\leq 0,$

i.e., $NQ^{-1}$ is maximally monotone. Furthermore, the multiplication
operator $x\mapsto(I-\lambda NQ^{-1})(I+\lambda NQ^{-1})^{-1}x$ is a
contraction on $L^{2}_{\omega}([0,T],\mathbb{R}^{n})$ for $\lambda>0$
satisfying

$\displaystyle\|(I-\lambda NQ^{-1})$ $\displaystyle(I+\lambda
NQ^{-1})^{-1}x\|_{2,\omega}\leq\|x\|_{2,\omega}\quad\forall\,x\in
L^{2}_{\omega}([0,T],\mathbb{R}^{n}).$

Moreover, if $\operatorname{im}E+\operatorname{im}R=\mathbb{R}^{n}$ and
$\alpha<1$, then

$\displaystyle\|(I-\lambda NQ^{-1})$ $\displaystyle(I+\lambda
NQ^{-1})^{-1}x\|_{2,\omega}\leq q\|x\|_{2,\omega}\quad\forall\,x\in
L^{2}_{\omega}([0,T],\mathbb{R}^{n}).$

for $q<1$ as in (21).

###### Proof.

The first statement follows directly and the remaining statements are implied
by Lemma 18. ∎

Our dynamic iteration scheme is built on Rem. 16. Following Lions and Mercier
[19], we study the algorithms

$\displaystyle x^{k+1}:=$ $\displaystyle Q^{-1}(I+\lambda
MQ^{-1})^{-1}(I-\lambda NQ^{-1})(I+\lambda NQ^{-1})^{-1}(I-\lambda
MQ^{-1})Qx^{k}.$ (22)

with $x^{0}\in D(M)$ arbitrary.

######  Remark 21.

Note that equation (16) is equivalent to

$(I+\lambda NQ^{-1})(I+\lambda MQ^{-1})Qx=(I-\lambda NQ^{-1})(I-\lambda
MQ^{-1})Qx,$

or equivalently,

$x=Q^{-1}(I+\lambda MQ^{-1})^{-1}(I-\lambda NQ^{-1})(I+\lambda
NQ^{-1})^{-1}(I-\lambda MQ^{-1})Qx.$

######  Remark 22.

The iteration (22) is equivalent to

$\displaystyle x^{k+1}:=$ $\displaystyle(Q+\lambda M)^{-1}(Q-\lambda
N)(Q+\lambda N)^{-1}(Q-\lambda M)x^{k}.$ (23)

######  Theorem 23.

Let $\lambda>0$, $\mu\geq 0$, and let $\alpha\in[0,1]$. Assume that the DAE
(2) fulfills Assumption 5, and let the operators $M$, $N$ be defined as in
(14) and (15). Let $\lambda>0$ and $x^{0}\in D(M)$. Let the sequence
$(x^{k})_{k}$ be defined by (22), and let $x$ be the unique solution of (2).
Then the following statements are fulfilled:

  1. (a)

$(Ex^{k})_{k}$ converges to $Ex$ in $L^{2}([0,T],\mathbb{R}^{n})$,

  2. (b)

$(Ex^{k})_{k}$ converges pointwise to $Ex$ on $(0,T]$,

  3. (c)

If $\alpha>0$, then $(RQx^{k})_{k}$ converges to $RQx$ in
$L^{2}([0,T],\mathbb{R}^{n})$,

  4. (d)

For the sequence $(z^{k})_{k}$ given by

$\forall\,k\in\mathbb{N}:\quad z^{k}:=(I+\lambda MQ^{-1})Qx^{k},$

and the function $z:=(I+\lambda MQ^{-1})Qx$, the real sequence
$(\|z^{k}-z\|_{2,\omega})_{k}$ is monotonically decreasing, and

$\forall\,k\in\mathbb{N}:\quad\|x^{k}-x\|_{2,\omega}\leq\|Q^{-1}\|\|z^{k}-z\|_{2,\omega}.$
(24)

###### Proof.

Let $\omega\in\mathbb{R}$ with $\mu<\omega$. We have

$\displaystyle z^{k+1}$ $\displaystyle=(I-\lambda NQ^{-1})(I+\lambda
NQ^{-1})^{-1}(2I-(I+\lambda MQ^{-1}))(I+\lambda MQ^{-1})^{-1}{z}^{k}$
$\displaystyle=(I-\lambda NQ^{-1})(I+\lambda NQ^{-1})^{-1}[2Qx^{k}-z^{k}]$
(25)

and

$z=(I-\lambda NQ^{-1})(I+\lambda NQ^{-1})^{-1}(2Qx-z).$

We set

$\Delta x^{k}:=x-x^{k},\qquad\Delta z^{k}:=z-z^{k}.$

Then using Lemma 20 and Lemma 17 we calculate

$\displaystyle\|\Delta$ $\displaystyle z^{k+1}\|_{2,\omega}^{2}-\|\Delta
z^{k}\|_{2,\omega}^{2}=\|z-z^{k+1}\|_{2,\omega}^{2}-\|\Delta
z^{k}\|_{2,\omega}^{2}$ $\displaystyle=\|(I-\lambda NQ^{-1})(I+\lambda
NQ^{-1})^{-1}[(2Qx-z)-(2Qx^{k}-z^{k})]\|_{2,\omega}^{2}-\|\Delta
z^{k}\|_{2,\omega}^{2}$ $\displaystyle\leq\|2Q\Delta x^{k}-\Delta
z^{k}\|_{2,\omega}^{2}-\|\Delta z^{k}\|_{2,\omega}^{2}$
$\displaystyle=4\|Q\Delta x^{k}\|_{2,\omega}^{2}-4\langle Q\Delta x^{k},\Delta
z^{k}\rangle_{2,\omega}$ $\displaystyle=4\|Q\Delta
x^{k}\|_{2,\omega}^{2}-4\langle Qx-Qx^{k},(I+\lambda MQ^{-1})Qx-(I+\lambda
MQ^{-1})Qx^{k}\rangle_{2,\omega}$ $\displaystyle=-4\lambda\langle Qx-
Qx^{k},Mx-Mx^{k}\rangle_{2,\omega}$ $\displaystyle=-2\lambda
e^{-2T\omega}\|(EQ^{-1})^{1/2}Q\Delta x^{k}(T)\|^{2}$
$\displaystyle\quad-4\lambda(\omega-\mu)\|(EQ^{-1})^{1/2}Q\Delta
x^{k}\|_{2,\omega}^{2}-4\lambda\alpha\|R^{1/2}Q\Delta x^{k}\|_{2,\omega}^{2}.$
(26)

Thus the sequence $(\|\Delta z^{k}\|_{2,\omega})_{k}$ monotonically
decreasing, and therefore it converges. Inequality (24) follows from Remark
15, which completes the proof of d. Furthermore, the latter inequality gives

$\displaystyle\|(EQ^{-1})^{1/2}Q\Delta x^{k}\|_{2,\omega}^{2}\leq$
$\displaystyle\,\frac{1}{4\lambda(\omega-\lambda)}(\|\Delta
z^{k}\|_{2,\omega}^{2}-\|\Delta z^{k+1}\|_{2,\omega}^{2}),$ (27)
$\displaystyle\|R^{1/2}Q\Delta x^{k}\|_{2,\omega}^{2}\leq$
$\displaystyle\,\frac{1}{4\lambda\alpha}(\|\Delta
z^{k}\|_{2,\omega}^{2}-\|\Delta z^{k+1}\|_{2,\omega}^{2}).$ (28)

Combining (28) with the fact that the norms $\|\cdot\|_{\omega,2}$ and
$\|\cdot\|_{2}$ are equivalent, we are led to (c). Further, (27) together with
$E=(EQ^{-1})^{1/2}(EQ^{-1})^{1/2}Q$ gives $(\|E\Delta
x^{k}\|_{2,\omega}^{2})_{k}$ converges to zero, and the equivalence between
the norms $\|\cdot\|_{\omega,2}$ and $\|\cdot\|_{2}$ now implies (a).
It remains to prove (b), using (26) for $t\in(0,T]$ instead of $T$ we obtain

$\|(EQ^{-1})^{1/2}Q\Delta
x^{k}(t)\|^{2}\leq\frac{e^{2t\omega}}{2\lambda}(\|\Delta
z^{k}\|_{2,\omega}^{2}-\|\Delta z^{k+1}\|_{2,\omega}^{2}),$

which shows that, for all $t\in(0,T]$, $Ex^{k}(t)\rightarrow Ex(t).$ ∎

######  Remark 24.

Though the sequence $(\|\Delta x^{k}\|_{2,\omega})_{k}$ is usually not
monotone decreasing, it is bounded by the monotone decreasing sequence
$(\|Q^{-1}\|\|z^{k}-z\|_{2,\omega})_{k}$ due to (24). $\Box$

######  Remark 25.

As the sequences $(x^{k})_{k}$, $(z^{k})_{k}$ are bounded in $L^{2}$, each of
them has a weakly convergent subsequence. $\Box$

######  Theorem 26.

Let $\lambda>0$, $\mu\geq 0$, and let $\alpha\in[0,1]$. Assume that the DAE
(2) fulfills Assumption 5, and let the operators $M$, $N$ be defined as in
(14) and (15). Let $\lambda>0$ and $x^{0}\in D(M)$. Let the sequence
$(x^{k})_{k}$ be defined by (22), and let $x$ be the unique solution of (2).
Then the following statements are fulfilled:

  1. (a)

If (19) holds, then the sequence $(z^{k})$ as defined in Theorem 23 (d)
converges to $z:=(I+\lambda MQ^{-1})Qx$ in $L^{2}([0,T],\mathbb{R}^{n})$.
Further, $(Ex^{k})_{k}$ converges to $Ex$ in $H^{1}([0,T],\mathbb{R}^{n})$.

  2. (b)

If (19) or

$\mathrm{rk}\begin{bmatrix}E&\alpha R\end{bmatrix}=n,$ (29)

is fulfilled, then $(x^{k})_{k}$ converges to $x$ in
$L^{2}([0,T],\mathbb{R}^{n})$,

###### Proof.

We start with the proof of (a): By (22) and Remark 21, we have

$\displaystyle z^{k+1}:=$ $\displaystyle\,(I-\lambda MQ^{-1})(I+\lambda
MQ^{-1})^{-1}(I-\lambda NQ^{-1})(I+\lambda NQ^{-1})^{-1}z^{k},$ $\displaystyle
z:=$ $\displaystyle\,(I-\lambda MQ^{-1})Q(I+\lambda MQ^{-1})^{-1}(I-\lambda
NQ^{-1})(I+\lambda NQ^{-1})^{-1}z,$

and we define $\Delta z^{k}:=z-z^{k}$. Assume that (19) holds. By invoking
Remark 15, Lemma 17 and Lemma 20, we obtain, for $q<1$ as in Lemma 18,

$\displaystyle\|\Delta z^{k+1}\|_{2,\omega}$ $\displaystyle=\|(I-\lambda
MQ^{-1})(I+\lambda MQ^{-1})^{-1}(I-\lambda NQ^{-1})(I+\lambda NQ^{-1})^{-1}z$
$\displaystyle\quad-(I-\lambda MQ^{-1})(I+\lambda MQ^{-1})^{-1}(I-\lambda
NQ^{-1})(I+\lambda NQ^{-1})^{-1}z^{k}\|_{2,\omega}$
$\displaystyle\leq\|(I-\lambda NQ^{-1})(I+\lambda NQ^{-1})^{-1}\Delta
z^{k}\|_{2,\omega}$ $\displaystyle\leq q\|\Delta z^{k}\|_{2,\omega}.$

This implies that $z^{k}$ converges in $L^{2}([0,T],\mathbb{R}^{n})$ to $z$.
Next we prove that $(Ex^{k})_{k}$ converges to $Ex$ in
$H^{1}([0,T],\mathbb{R}^{n})$. We already know from Theorem 23 (a) that
$(Ex^{k})_{k}$ converges to $Ex$ in $L^{2}([0,T],\mathbb{R}^{n})$. Hence, it
suffices to prove that $(\tfrac{\mathrm{d}}{\mathrm{d}t}Ex^{k})_{k}$ converges
to $\tfrac{\mathrm{d}}{\mathrm{d}t}Ex$ in $L^{2}([0,T],\mathbb{R}^{n})$. By
invoking that $(z^{k})$ converges to $z$ in $L^{2}([0,T],\mathbb{R}^{n})$,
(24) yields that $(x^{k})$ converges to $x$ in $L^{2}([0,T],\mathbb{R}^{n})$.
Now using the definition of $x$, $z$ and $M$, we have

$\displaystyle z^{k}=$
$\displaystyle\,Qx^{k}+\lambda\big{(}\tfrac{\mathrm{d}}{\mathrm{d}t}(E{x}^{k})-(J_{d}-\alpha
R+\mu EQ^{-1})Qx^{k}-Bu),$ $\displaystyle z=$
$\displaystyle\,Qx+\lambda\big{(}\tfrac{\mathrm{d}}{\mathrm{d}t}(E{x})-(J_{d}-\alpha
R+\mu EQ^{-1})Qx^{k}-Bu),$

which gives

$\tfrac{\mathrm{d}}{\mathrm{d}t}(E{x})-\tfrac{\mathrm{d}}{\mathrm{d}t}(E{x}^{k})=\tfrac{1}{\lambda}(z-z^{k})+\tfrac{1}{\lambda}Q(x-x^{k})+(J_{d}-\alpha
R+\mu EQ^{-1})Q(x-x^{k}).$

Since $(x^{k})$ converges to $x$ in $L^{2}([0,T],\mathbb{R}^{n})$, and
further, $(z^{k})$ converges to $z$ in $L^{2}([0,T],\mathbb{R}^{n})$, the
above equation implies that $\tfrac{\mathrm{d}}{\mathrm{d}t}(E{x}^{k})$
converges in $L^{2}([0,T],\mathbb{R}^{n})$ to
$\tfrac{\mathrm{d}}{\mathrm{d}t}(E{x})$. Altogether, this means that
$(Ex^{k})$ converges to $Ex$ in $H^{1}([0,T],\mathbb{R}^{n})$.

Next we prove (b). The case where (19) is fulfilled, the result follows by a
combination of (a) with (24).
Assume that (29) holds: Then, by invoking Assumption 5 (a) (i), the matrix
$Q^{\top}E+\alpha Q^{\top}RQ$ is positive definite, and thus invertible. On
the other hand, $(((Q^{\top}E+Q^{\top}RQ)x^{k})_{k}$ converges to
$(Q^{\top}E+\alpha Q^{\top}RQ)x$ in $L^{2}([0,T],\mathbb{R}^{n})$ by Theorem
23. Consequently, $(x^{k})_{k}$ converges to $x$ in
$L^{2}([0,T],\mathbb{R}^{n})$. ∎

Based on the results of the previous section we develop the following
iteration scheme. We recall that we choose $\omega>0$ and we define

$x^{k+1}:=(Q+\lambda M)^{-1}(Q-\lambda N)(Q+\lambda N)^{-1}(Q-\lambda
M)x^{k},$

with $x^{0}\in D(M)$ arbitrary, and

$z^{k}:=(I+\lambda MQ^{-1})Qx^{k}.$

Using (25) we obtain

$z^{k+1}=(I-\lambda NQ^{-1})(I+\lambda NQ^{-1})^{-1}[2Qx^{k}-z^{k}].$

Furthermore, relation $z^{k}=(Q+\lambda M)x^{k}$ corresponds to the DAEs

$\tfrac{d}{dt}E_{i}x_{i}^{k}=(J_{i}-R_{i}+\mu
E_{i}Q_{i}^{-1})Q_{i}x_{i}^{k}-\tfrac{1}{\lambda}Q_{i}x_{i}^{k}+B_{i}u+\tfrac{1}{\lambda}z_{i}^{k},\quad
i=1,\ldots,s,$

with initial conditions $E_{i}x_{i}^{n}(0)=E_{i}x_{0,i}$ ($i=1,\ldots,s$).

######  Remark 27.

  1. (a)

If $\mathrm{rk}\begin{bmatrix}E&R\end{bmatrix}<n$, then it is not possible to
formulate any convergence results for $(x^{k})$, in general. As an example,
consider a system (2) with

$E=R=0_{2\times 2},\;Q=I_{2},\;J=\left[\begin{smallmatrix}0&1\\\
-1&0\end{smallmatrix}\right],\;u=0\in L^{2}([0,T];\mathbb{R}^{2}),$

which belongs to the class specified in Assumption 5 with $s=2$ and
$n_{1}=n_{2}=1$. The unique solution of (2) is clearly given by $x=0$ (the
specification of the initial value $x_{0}$ is obsolete by $E=0$).
The above choice of $E$, $R$ and $Q$ leads to $M$ in (14) being the zero
operator, whereas $N$ as in (15) corresponds to the multiplication with the
skew-Hermitian matrix $-J$. Consequently, for $\lambda>0$, the iteration (23)
now reads

$x^{k+1}:=Ux^{k},\quad\text{where}\quad U=(I+\lambda J)(I-\lambda J)^{-1}.$

By invoking Lemma 18, we see that $U$ is a unitary matrix with $U\neq I_{2}$.
Therefore, $\|x^{k+1}(t)\|=\|x^{k}(t)\|$ for all $k\in\mathbb{N}$,
$t\in[0,T]$, and thus $\|x^{k+1}\|_{2}=\|x^{k}\|_{2}$. It can be concluded
that the sequence $(x^{k})$ does not converge in $L^{2}([0,T],\mathbb{R}^{2})$
to $x=0$. A closer look to the iteration yields that $(x^{k})$ does not even
have a subsequence which weakly converges in $L^{2}([0,T],\mathbb{R}^{2})$ to
$x=0$.

  2. (b)

Under Assumption 5, (19) is equivalent to $E$ and $R$ fulfilling
$\mathrm{rk}E=n$, or $\mathrm{rk}\begin{bmatrix}E&R\end{bmatrix}=n$ and
$\alpha>0$.

  3. (c)

In the case of ordinary differential equations, i.e., $E=I$, then Theorem 26
simplifies to the following:

    1. (i)

If $\mu>0$, or $\mathrm{rk}R=n$ and $\alpha<1$, then $(z^{k})$ converges to
$z$ in $L^{2}([0,T],\mathbb{R}^{n})$, and $(x^{k})_{k}$ converges to $x$ in
$H^{1}([0,T],\mathbb{R}^{n})$.

    2. (ii)

$(x^{k})_{k}$ converges to $x$ in $L^{2}([0,T],\mathbb{R}^{n})$. $\Box$

The above discussion results in a Lions-Mercier-type dynamic iteration for
PHS, which give summarize in Alg. 5.

Lions-Mercier-type Dynamic Iteration for PHS 1:problem data: $J_{i}$, $J_{o}$,
$R_{i}$, $Q_{i}$, $E_{i}$, $B_{i}$, $u_{i}$, $x_{i,0}$ ($i=1,\dotsc,s$)
2:choose for $i=1,\ldots,s$ $z^{0}_{i}\in L^{2}([0,T];\,\mathbb{R}^{n_{i}})$ %
initial guess  3:choose parameters $\alpha$, $\mu$, $\lambda$
4:$R\,\leftarrow\operatorname{diag}(R_{1},\dotsc,R_{s})$, $E,Q$ analogously
5:$K\,\leftarrow\,(1-\alpha)R+\mu EQ^{-1}-J_{o}$
6:$z^{0}\,\leftarrow\,\begin{bmatrix}(z_{1}^{0})^{\top},\,\dotsc,\,(z_{s}^{0})^{\top}\end{bmatrix}^{\top}$
7:for $k=0,1,2,\ldots$ do 8: for $i=1,2,\dotsc,\,s$ do 9: solve for
$x_{i}^{k}$ $\displaystyle\tfrac{d}{dt}E_{i}x_{i}^{k}=(J_{i}-\alpha R_{i}+\mu
E_{i}Q_{i}^{-1})Q_{i}x_{i}^{k}-\tfrac{1}{\lambda}Q_{i}x_{i}^{k}+B_{i}u+\tfrac{1}{\lambda}z_{i}^{k}$
$\displaystyle\text{with initial value }E_{i}x_{i}^{k}(0)=E_{i}x_{0,i}$ 10:
end for 11:
$x^{k}\,\leftarrow\,\begin{bmatrix}(x_{1}^{k})^{\top},\,\dotsc,\,(x_{s}^{k})^{\top}\end{bmatrix}^{\top}$
12: $z^{k+1}\,\leftarrow\,(I-\lambda K)(I+\lambda
K)^{-1}\left(2Qx^{k}-z^{k}\right)$ 13: partition $z^{k+1}$ into components
$z_{1}^{k+1},\,\dotsc,\,z_{s}^{k+1}$ 14: if error small then return  15: end
if 16:end for

##  6 Numerical results

We analyse the convergence rates for the new type of monotonic dynamic
iterations (Lions-Mercier-type iteration). Then, we show numerical convergence
results for the promising monotone convergence algorithms developed here.

###  6.1 Convergence rate discussion

For $0<\alpha<1$ and $K=((1-\alpha)R+\mu EQ^{-1}-J_{o})$ the convergence rate
is given by the contraction factor (21)

$q^{2}=1-\frac{4\lambda}{(1+\lambda\|K\|)^{2}\|(K+J_{o})^{-1}\|}.$

The optimal $\lambda^{*}$ will have the smallest error reduction factor. We
find:

$\lambda^{*}(\mu)=\frac{1}{\|K\|}\quad\;\text{and}\quad\;q^{\ast}\bigl{(}\mu,\,\lambda^{\ast}(\mu)\bigr{)}^{2}=1-\frac{1}{\|K\|\cdot\|(K+J_{o})^{-1}\|}.$

######  Remark 28 (Decoupled setting).

In the decoupled case, $J_{o}=0$, and hence the optimal error reduction reads:

$q^{\ast}\bigl{(}\mu,\,\lambda^{\ast}(\mu)\bigr{)}^{2}=1-\frac{1}{\text{cond}{(K)}},$

which is only small for nearly equilibrated matrices $K$. This is in contrast
to (e.g.) Jacobi dynamic iteration procedures, where we have immediate
convergence. $\Box$

Considering ODEs with $E=I$ and setting $\alpha=1$, we have $K=\mu
Q^{-1}-J_{o}$.

######  Remark 29 (ODE case).

  1. a)

Transformation to $Q=I$: using the Cholesky factorization $Q=V^{\top}V$, we
multiply (2) with $V$, which yields

$\displaystyle V\tfrac{\mathrm{d}}{\mathrm{d}t}x$
$\displaystyle=V(J-R)(V^{\top}V)x+VBu(t)\qquad\Leftrightarrow\quad\tfrac{\mathrm{d}}{\mathrm{d}t}\tilde{x}=(\tilde{J}-\tilde{R})\tilde{x}+\tilde{B}u(t)$

with $\tilde{x}=Vx$, $\tilde{J}=VJV^{\top},\tilde{R}=VRV^{\top}$ and
$\tilde{B}=VB$ provides a remedy. Correspondingly, $y=B^{\top}Qx$ is
transformed into $\tilde{y}={\tilde{B}}^{\top}\tilde{x}$.

  2. b)

Decoupled setting ($J_{o}=0$): The reduction factor is
$q^{\ast}=1-\tfrac{1}{\text{cond}{(Q)}}$. Again, this is only small for
scaling matrices $Q$ with small condition number. The above transformation to
$Q=I$ yields $q^{*}=0$. An alternative is a shift with a matrix-valued
$\lambda$ in Alg. 5.

  3. c)

Coupled setting: for $Q=I$ and the Euclidean norm this gives

$\displaystyle\lambda^{*}(\mu)$
$\displaystyle=\tfrac{1}{\sqrt{\mu^{2}+\lambda_{\max}(J_{o}^{\top}J_{o})}}\quad\text{and}\quad
q^{*}\bigl{(}\mu,\lambda^{*}(\mu)\bigr{)}^{2}=1-\tfrac{1}{\sqrt{1+{\lambda_{\max}(J_{o}^{\top}J_{o})}/{\mu^{2}}}}.$
(30)

###  6.2 Results for monotone Lions-Mercier-type algorithm

####  6.2.1 Example with Jacobi failure

For system (12) holds $Q=I$ and $\lambda_{\max}(J_{o}^{\top}J_{o})=\nu^{2}$.
This gives the optimal contraction constant via (30):
$(q^{*})^{2}={1-\tfrac{1}{\sqrt{1+{\nu^{2}}/{\mu^{2}}}}}$. We see the
following:

######  Remark 30.

  1. a)

We have monotone convergence for the Lions-Mercier-type algorithm, Alg. 5.

  2. b)

In the limit of an undamped error bound $\omega\rightarrow 0^{+}$, the optimal
convergence rate goes to one (for $\nu\neq 0$). This indicates an extremely
slow convergence.

  3. c)

One can get arbitrarily small convergence rates $q^{*}$ for
$\mu\rightarrow\infty$. However, this corresponds to more and more damped
error norms, which do not allow to estimate the numerically important undamped
error. $\Box$

####  6.2.2 Simple 2x2 system with scaling

We introduce into the simple 2x2 example (12) a scaling matrix $Q$:

$x^{\prime}=(J\\!-\\!R)Qx,\quad\\!x(0)=x_{0},\quad\\!\text{with}\\!\\!\quad\\!J\\!-\\!R=\begin{bmatrix}-\tau&\nu\\\
-\nu&-\tau\end{bmatrix}\\!,\quad\\!Q=\begin{bmatrix}q_{1}&0\\\
0&q_{2}\end{bmatrix}\\!,$ (31)

where $\tau,\nu>0$ and $q_{1}\geq q_{2}>0$. In the decoupled case, $\nu=0$,
the optimal reduction factor is given by (see Rem. 29b):

$q(\lambda^{*})=1-\tfrac{1}{\text{cond}(Q)}=1-\tfrac{q_{2}}{q_{1}}$

Fig. 2 gives numerical results for the Lions-Mercier-type algorithm, showing
that $Q$ not being a multiple of $I$ destroys the convergence in one step,
although the system is fully decoupled.

|
---|---
(a) decoupled, $Q\\!=\\!I$, non-optimal $\lambda$ | (d) coupled, $Q\\!=\\!I$, non-optimal $\lambda$
Decoupled with optimal $\lambda$: the error reduced directly to zero (no logarithmic plot). |
(b) decoupled, $Q=I$, optimal $\lambda$ | (e) coupled, $Q=I$, optimal $\lambda$
|
(c) decoupled, $Q\neq I$, optimal $\lambda$ | (f) coupled, $Q\neq I$, optimal $\lambda$
$q_{1}=1.5$, $q_{2}=1$, $q(\lambda^{*})\approx 0.3333$ |  $q_{1}=1.5$, $q_{2}=1$, $q(\lambda^{*})\approx\frac{1}{3}$
Figure 2: Lions-Mercier-type iteration with (31): number of iterations ($x$)
versus scaled error in coupling term $v$/$z$ for various settings. First
column provide the decoupled case ($\nu=0$), and the right column the coupled
cases ($\nu=1.5$). Always $\tau=-0.5$ and $\omega=2.1$.

######  Remark 31.

The convergence rates for the coupled cases – even for the optimal choice of
$\lambda$ should be improved to have more competitive dynamic iteration
schemes. $\Box$

####  6.2.3 Two masses and three springs example with damping

Figure 3:  ODE two masses oscillator with damping. The coordinates
$q_{1},\,q_{2}$ describe the position of the masses.

Now, we consider two masses $m_{1},\,m_{2}>0$, which are connected via
massless springs $K_{1},\,K_{2},K_{3}>0$ to walls with damping
$r_{1},r_{2}>0$, see Fig. 3. To set up the system, we have positions
$q_{1},\,q_{2}$ and momenta $p_{1},\,p_{2}$ for the masses $m_{1},\,m_{2}$,
respectively. Then, the system can be modeled by the following Hamiltonian
equation of motion:

$\displaystyle\dot{q}_{1}$ $\displaystyle=\tfrac{1}{m_{1}}p_{1},$
$\displaystyle\qquad\dot{q}_{2}$ $\displaystyle=\tfrac{1}{m_{2}}p_{2},$
$\displaystyle\dot{p}_{1}$
$\displaystyle=-(K\\!+\\!K_{1})q_{1}+Kq_{2}\\!-\\!\tfrac{r_{1}}{m_{1}}p_{1},$
$\displaystyle\dot{p}_{2}$
$\displaystyle=-(K\\!+\\!K_{2})q_{2}+Kq_{1}\\!-\\!\tfrac{r_{2}}{m_{2}}p_{2}.$

Notice, this is not yet given in our PHS-format (2) for the overall systems.
Before we treat the overall system, we set up a description as coupled PH-
subsystem as in Sec. 2.2, let $q=q_{2}$ denote the coupling variables, which
we need for the first subsystem. Then, the two coupled PHS systems read:

$\displaystyle\dot{x}_{1}$ $\displaystyle=(J_{1}-R_{1})z_{1}+B_{1}u_{1},$
$\displaystyle\qquad\ \qquad\qquad\quad\dot{x}_{2}$
$\displaystyle=(J_{2}-R_{2})z_{2}+B_{2}u_{2},$ $\displaystyle y_{1}$
$\displaystyle=B_{1}^{\top}z_{1},$ $\displaystyle y_{2}$
$\displaystyle=B_{2}^{\top}z_{2},$ $\displaystyle u$ $\displaystyle=Cy,$

with quantities for the first subsystem:
$x_{1}=[p_{1},\,q_{1},\,q_{1}-q]^{\top}$, $z_{1}=Q_{1}x_{1}$,

$\displaystyle Q_{1}\\!=\\!\begin{bmatrix}\frac{1}{m_{1}}&0&0\\\ 0&K_{1}&0\\\
0&0&K\end{bmatrix}\\!,\;J_{1}\\!=\\!\begin{bmatrix}0&-1&-1\\\ 1&0&0\\\
1&0&0\end{bmatrix}\\!,\;R_{1}\\!=\\!\begin{bmatrix}r_{1}&0&0\\\ 0&0&0\\\
0&0&0\end{bmatrix}\\!,\;B_{1}\\!=\\!\begin{bmatrix}0\\\ 0\\\ -1\end{bmatrix}$

and for the second subsystem we have: $x_{2}=[p_{2},\,q_{2}]^{\top}$,
$z_{2}=Q_{2}x_{2}$

$\displaystyle Q_{2}\\!=\\!\begin{bmatrix}\frac{1}{m_{2}}&0\\\
0&K_{2}\end{bmatrix}\\!,\;J_{2}\\!=\\!\begin{bmatrix}0&-1\\\
1&0\end{bmatrix}\\!,\quad\;R_{2}\\!=\\!\begin{bmatrix}r_{2}&0\\\
0&0\end{bmatrix}\\!,\quad B_{2}\\!=\\!\begin{bmatrix}1\\\ 0\end{bmatrix}\\!,$
and the coupling interface reads: $\displaystyle u$
$\displaystyle\\!=\\!\begin{bmatrix}u_{1}\\\ u_{2}\end{bmatrix},\quad
y=\begin{bmatrix}y_{1}\\\ y_{2}\end{bmatrix},\quad C=\begin{bmatrix}0&1\\\
-1&0\end{bmatrix}.$

Also the overall coupled system can be written in our pH-structure given in
(2), with $F=\operatorname{diag}(F_{1},F_{2})$ for $F\in\\{J,R,Q,B\\}$, cf.
Rem. 7:

$\displaystyle\begin{bmatrix}x_{1}\\\ x_{2}\end{bmatrix}^{\prime}$
$\displaystyle=\left([J-BCB^{\top}]-R\right)\cdot Q\begin{bmatrix}x_{1}\\\
x_{2}\end{bmatrix}.$ (32)

Numerical results for the Lions-Mercier-type Alg. 5 are given in Fig. 4. These
results show that the estimate $q^{\star}\approx 0.975$ for the convergence
rate (with an optimal choice leading to $q^{\star}\approx 0.944$) is quite
pessimistic in this case (for $\mu=2$).

Figure 4: Convergence analysis for Lions-Mercier-type iteration
($\lambda=1.5$, $\mu=2$ (first row) $\mu=0.1$ (second row), $\omega=2.2$,
$\alpha=0.5$) for system (32) with $K_{i}=2$, $m_{i}=2$ ($i=1,2$), $K=4$,
$r_{1}=0.5$, $r_{2}=0.75$, initial values are zero apart from $p_{2,0}=0.1$.

######  Remark 32.

Including coupling and output quantities as variables, i.e.,
$\tilde{x}_{i}=\begin{bmatrix}x_{i},\,u_{i},\,y_{i}\end{bmatrix}^{\top}$
yields a DAE, which has no damping ($R$) in the coupling or output equations.
Thus it does not fit into our framework. $\hfill\Box$

Finally, we like to discuss a DAE example. To have a DAE multibody system, a
simple choice are holonomic constraints. However, such a constraint comes
without dissipation. Hence $\text{rk}(E,R)$ is not full. Therefore, we change
the field of applications.

###  6.3 Circuit example (DAE)

The electric circuit in Fig. 5 can be modelled as follows: for
$x^{\top}=\begin{bmatrix}u_{1},\,u_{2},\,u_{3},\,\jmath_{1},\,u_{4},\,u_{5}\end{bmatrix}$
and $E=\text{diag}(0,C_{1},0,L_{1},0,C_{2})$
$B^{\top}=\begin{bmatrix}1,\,0,\,0,\,0,\,0,\,0\end{bmatrix}$

$\displaystyle(Ex)^{\prime}=B\imath(t)+$
$\displaystyle\left(\\!\begin{bmatrix}\begin{array}[]{cccc|cc}0&&&&&\\\
&0&&&&\\\ &&0&1&&\\\ &&-1&0&1&\\\ \hline\cr&&&-1&0&\\\
&&&&&0\end{array}\end{bmatrix}\\!-\\!\begin{bmatrix}\begin{array}[]{cccc|cc}\tfrac{1}{R_{1}}&-\tfrac{1}{R_{1}}&&&&\\\
-\tfrac{1}{R_{1}}&\tfrac{1}{R_{1}}\\!+\\!\tfrac{1}{R_{2}}&-\tfrac{1}{R_{2}}&&&\\\
&-\tfrac{1}{R_{2}}&\tfrac{1}{R_{2}}+\tfrac{1}{R_{3}}&&&\\\ &&&0&&\\\
\hline\cr&&&&\tfrac{1}{R_{4}}&-\tfrac{1}{R_{4}}\\\
&&&&-\tfrac{1}{R_{4}}&\tfrac{1}{R_{4}}\\!+\\!\tfrac{1}{R_{5}}\end{array}\end{bmatrix}\\!\right)x.$

Figure 5: Example of a simple electric circuit; $R_{i}>0$, $L_{i}>0$,
$C_{i}>0$ and given $\imath\in L^{2}[0,T]$.

This example fulfills the rank condition $\text{rk}(E,R)=n=6$ and has block-
diagonal $R$.

Numerical results for the Lions-Mercier-type Alg. 1 are given in Fig. 6. These
results show that the estimate $q\approx 0.9998$ for the convergence rate
(with an optimal choice leading to $q^{\star}\approx 0.9980$) is quite
pessimistic in this case, too.

(a) error reduction: achieved precision vs. iteration count
(left in $x$, center in $z$, right in $z$ with weighted $L^{2}$-norm);

(b) convergence in $x_{3}$: we give results for iteration $l=1,10,20,30$;

(c) convergence in $x_{5}$: we give results for iteration $l=1,10,20,30$.

Figure 6: Results for the academic circuit Fig. 5 using:
$R_{1}=R_{2}=R_{3}=R_{4}=0.5,\,R_{5}=5$, $C_{1}=C_{2}=5\cdot 10^{-4}$, $L=20$,
$i(t)=\sin(2\pi\cdot 50t)\cdot\sin(2\pi\cdot 500t)$. Lions-Mercier parameters:
$\lambda=1.2$, $\mu=1$, $\omega=2.2$, $\alpha=0.2$. (a) shows the $L^{2}$
error reduction for different variables, (b) and (c) depict the convergence in
a variable of the first subsystem $x_{3}$ and in the second subsystem $x_{5}$,
respectively.

######  Remark 33 (Rank condition for circuits).

We discuss circuit equation modelled by modified nodal analysis.
i) Let a circuit include one voltage sources. Then the model has a pure
algebraic equation, which has no dissipation. Thus the rank condition is not
satisfied.
ii) Electric circuits with current sources, capacitors, inductors and
resistors can be made to fulfill the rank condition by adding sufficiently
many resistors.
iii) The coupling—without auxiliary variables—can be facilitated via
inductors. $\Box$

##  7 Conclusions

For a general PHS setting, we have shown existence and uniqueness of weak
solutions. Furthermore, the perspective of coupled systems and coupling
variables are treated. Then, it is demonstrated that Jacobi-type of dynamic
iteration may fail in practice due to finite precision. To circumvent this, a
Mercier-Lions-type dynamic iteration scheme is developed, which guarantees
monotone convergence in a related variable. The proof is based on properties
of the Cayley transform. Numerical results demonstrate the monotone
convergence. As a drawback, we observe that the currently achieved optimal
convergence rates still need some improvements in order to pave the way for a
broader application.

## References

  * [1] Giuseppe Alì, Andreas Bartel, Michael Günther, Vittorio Romano, and Sebastian Schöps.  Simulation of coupled PDAEs: Dynamic iteration and multirate simulation.  In Michael Günther, editor, Coupled Multiscale Simulation and Optimization in Nanoelectronics, pages 103–156. Springer Berlin Heidelberg, Berlin, Heidelberg, 2015.
  * [2] Hans Wilhelm Alt.  Linear Functional Analysis.  Universitext. Springer-Verlag, London, 2016.
  * [3] Martin Arnold.  Modular time integration of coupled problems in system dynamics.  In Michael Günther and Wil Schilders, editors, Novel Mathematics Inspired by Industrial Challenges, pages 57–72. Springer International Publishing, Cham, 2022.
  * [4] Martin Arnold and Michael Günther.  Preconditioned dynamic iteration for coupled differential-algebraic systems.  BIT Numer. Math., 41:1–25, 2001.
  * [5] V. Barbu.  Nonlinear Differential Equations of Monotone Types in Banach Spaces.  Springer Monographs in Mathematics. Springer-Verlag, New York, 2010.
  * [6] Andreas Bartel, Markus Brunk, Michael Günther, and Sebastian Schöps.  Dynamic iteration for coupled problems of electric circuits and distributed devices.  SIAM Journal on Scientific Computing, 35(2):B315–B335, 2013.
  * [7] Andreas Bartel, Markus Brunk, and Sebastian Schöps.  On the convergence rate of dynamic iteration for coupled problems with multiple subsystems.  Journal of Computational and Applied Mathematics, 262:14–24, 2014\.
  * [8] Andreas Bartel and Michael Günther.  PDAEs in refined electrical network modeling.  SIAM Review, 60(1):56–91, 2018.
  * [9] C. Beattie, V. Mehrmann, H. Xu, and H. Zwart.  Linear port-hamiltonian descriptor systems.  Math. Control Signals Systems, 30(4):17, 2018.
  * [10] Kevin Burrage.  Parallel and Sequential Methods for Ordinary Differential Equations.  Clarendon Press, Oxford, 1995.
  * [11] J. Cervera, A. J. van der Schaft, and A. Baños.  Interconnection of port-Hamiltonian systems and composition of Dirac structures.  Automatica, 43(2):212–225, 2007.
  * [12] H. Gernandt, F.E. Haller, and T. Reis.  A linear relation approach to port-Hamiltonian differential-algebraic equations.  SIAM J. Matrix Anal. Appl., 42(2):1011–1044, 2020.
  * [13] H. Gernandt, F.E. Haller, T. Reis, and A.J. van der Schaft.  Port-Hamiltonian formulation of nonlinear electrical circuits.  J. Geom. Phys., 159:102959, 2020.
  * [14] Michal Günther, Andreas Bartel, Birgit Jacob, and Timo Reis.  Dynamic iteration schemes and port-Hamiltonian formulation in coupled DAE circuit simulation.  International Journal of Circuit Theory and Applications, 49:430–452, 2021.
  * [15] Z. Jackiewicz and M. Kwapisz.  Convergence of waveform relaxation methods for differential-algebraic systems.  SIAM Journal on Numerical Analysis, 33(6):2303–2317, 1996.
  * [16] D. Jeltsema and A.J. van der Schaft.  Port-Hamiltonian systems theory: An introductory overview.  Foundations and Trends in Systems and Control, 1(2-3):173–387, 2014\.
  * [17] Peter Kunkel and Volker Mehrmann.  Differential-Algebraic Equations. Analysis and Numerical Solution.  EMS Publishing House, Zürich, Switzerland, 2006.
  * [18] René Lamour, Roswitha März, and Caren Tischendorf.  Differential Algebraic Equations: A Projector Based Analysis, volume 1 of Differential-Algebraic Equations Forum.  Springer-Verlag, Heidelberg-Berlin, 2013.
  * [19] P. L. Lions and B. Mercier.  Splitting algorithms for the sum of two nonlinear operators.  SIAM Journal on Numerical Analysis, 16(6):964–979, 1979.
  * [20] B. Maschke and A.J. van der Schaft.  Generalized port-Hamiltonian DAE systems.  System & Control Letters, 121:31–37, 2018.
  * [21] B. Maschke and A.J. van der Schaft.  Dirac and Lagrange algebraic constraints in nonlinear port-Hamiltonian systems.  Vietnam J. Math., 48:929–939, 2020.
  * [22] C. Mehl, V. Mehrmann, and M. Wojtylak.  Linear algebra properties of dissipative Hamiltonian descriptor systems.  SIAM Journal on Matrix Analysis and Applications, 39(3):1489–1519, 2018.
  * [23] Volker Mehrmann and Riccardo Morandin.  Structure-preserving discretization for port-Hamiltonian descriptor systems.  In 2019 IEEE 58th Conference on Decision and Control (CDC), pages 6863–6868, 2019.
  * [24] A.J. van der Schaft.  Port-Hamiltonian systems: Network modeling and control of nonlinear physical systems.  In K. Schlacher and H. Irschnik, editors, Advanced Dynamics and Control of Structures and Machines, volume 444 of CISM courses and lectures, pages 127–167. Springer, Vienna, 2004.
  * [25] A.J. van der Schaft.  Port-Hamiltonian differential-algebraic systems.  In Achim Ilchmann and Timo Reis, editors, Surveys in Differential-Algebraic Equations I, Differential-Algebraic Equations Forum, pages 173–226. Springer, Berlin-Heidelberg, 2013.