11institutetext: Andreas Frommer and Michael Günther 22institutetext:
Bergische Universität Wuppertal, IMACM, Gaußstraße 20, D-42119 Wuppertal,
22email<EMAIL_ADDRESS>
Björn Liljegren-Sailer and Nicole Marheineke 33institutetext: Universität
Trier, FB IV – Mathematik, Universitätsring 15, D-54296 Trier,
33email<EMAIL_ADDRESS>
# Operator splitting for port-Hamiltonian systems

Andreas Frommer  Michael Günther  Björn Liljegren-Sailer
Nicole Marheineke

###### Abstract

The port-Hamiltonian approach presents an energy-based modeling of dynamical
systems with energy-conservative and energy-dissipative parts as well as an
interconnection over the so-called ports. In this paper, we apply an operator
splitting that treats the energy-conservative and energy-dissipative parts
separately. This paves the way for linear equation solvers to exploit the
respective special structures of the iteration matrices as well as the
multirate potential in the different right-hand sides. We illustrate the
approach using test examples from coupled multibody system dynamics.

##  1 Introduction

Operator splitting mclachlan_quispel_2002  is an efficient tool to numerically
solve initial-value problems of ordinary differential equations (ODE-IVPs)

$\dot{x}=f(x,t)=f_{1}(x,t)+f_{2}(x,t),\qquad x(t_{0})=x_{0},$ (1)

that allow for a splitting of the right hand side $f(x,t)$ into two parts
$f_{1}(x,t)$ and $f_{2}(x,t)$ of profoundly different behaviour, e.g., with
respect to stiffness, computational costs, dynamics etc. Rewriting the system
in the homogeneous from

$\displaystyle\begin{pmatrix}\dot{x}\\\ \dot{s}\end{pmatrix}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}f(x,s)\\\
1\end{pmatrix}=\underbrace{\begin{pmatrix}f_{1}(x,s)\\\
1\end{pmatrix}}_{\displaystyle\tilde{f}_{1}(x,s):=}+\underbrace{\begin{pmatrix}f_{2}(x,s)\\\
0\end{pmatrix}}_{\displaystyle\tilde{f}_{2}(x,s):=},\quad\begin{pmatrix}x(t_{0})\\\
s(t_{0})\end{pmatrix}=\begin{pmatrix}x_{0}\\\ 0\end{pmatrix},$ (2)

the idea is to alternately solve dynamical systems driven by $\tilde{f}_{1}$
and $\tilde{f}_{2}$, respectively.

Strang splitting strang1968  solves the ODE system with respect to the first,
second and again first split right-hand side with step sizes $h/2$, $h$ and
$h/2$, where the initial values are given by the respective final values,
i.e.,

$\displaystyle\begin{pmatrix}\dot{x}_{1/3}\\\ \dot{s}_{1/3}\end{pmatrix}$
$\displaystyle=\tilde{f}_{1}(x_{1/3},s_{1/3}),\quad$
$\displaystyle\begin{pmatrix}x_{1/3}(t_{0})\\\
s_{1/3}(t_{0})\end{pmatrix}=\begin{pmatrix}x_{0}\\\ 0\end{pmatrix},$ (3)
$\displaystyle\begin{pmatrix}\dot{x}_{2/3}\\\ \dot{s}_{2/3}\end{pmatrix}$
$\displaystyle=\tilde{f}_{2}(x_{2/3},s_{2/3})\quad$
$\displaystyle\begin{pmatrix}x_{2/3}(t_{0})\\\
s_{2/3}(t_{0})\end{pmatrix}=\begin{pmatrix}x_{1/3}(h/2)\\\ h/2\end{pmatrix},$
(4) $\displaystyle\begin{pmatrix}\dot{x}_{1}\\\ \dot{s}_{1}\end{pmatrix}$
$\displaystyle=\tilde{f}_{1}(x_{1},s_{1})\quad$
$\displaystyle\begin{pmatrix}x_{1}(t_{0})\\\
s_{1}(t_{0})\end{pmatrix}=\begin{pmatrix}x_{2/3}(h)\\\ h\end{pmatrix}.$ (5)

This scheme provides an order two approximation $x_{1}(h/2)$ for the exact
solution $\varphi(h)\circ\begin{pmatrix}x_{0}\\\ 0\end{pmatrix}$ of (1) at
time point $t_{0}+h$, starting at the initial value $x_{0}$ at $t_{0}$. In
short hand:

$\begin{pmatrix}x_{1}(h/2)\\\
s_{1}(h/2)\end{pmatrix}=\varphi_{\tilde{f}_{1}}(\frac{h}{2})\circ\varphi_{\tilde{f}_{2}}(h)\circ\varphi_{\tilde{f}_{1}}(\frac{h}{2})\circ\begin{pmatrix}x_{0}\\\
0\end{pmatrix}=\varphi(h)\circ\begin{pmatrix}x_{0}\\\
0\end{pmatrix}+\begin{pmatrix}{\mathcal{O}}(h^{3})\\\ 0\end{pmatrix}$ (6)

where $\varphi_{\tilde{f}_{1}}(h)$ and $\varphi_{\tilde{f}_{2}}(h)$ denote the
solution of the flow $\tilde{f}_{1}$ and $\tilde{f}_{2}$, respectively, at
time point $t_{0}+h$, starting at $t_{0}$.

If the three initial-value problems (3)-(5) are solved by consistent
approximation schemes $\Phi_{\tilde{f}_{1}}(h/2)$, $\Phi_{\tilde{f}_{2}}(h)$
and $\Phi_{\tilde{f}_{1}}(h/2)$ with step sizes $h/2$, $h$ and $h/2$, the
corresponding numerical approximation

$\tilde{\Phi}(h)\circ\begin{pmatrix}x_{0}\\\
0\end{pmatrix}:=\Phi_{\tilde{f}_{1}}(h/2)\circ\Phi_{\tilde{f}_{2}}(h)\circ\Phi_{\tilde{f}_{1}}(h/2)\circ\begin{pmatrix}x_{0}\\\
0\end{pmatrix}$

is of second order due to the symmetry of the approach. We will apply this
idea of operator splitting to the special case of port-Hamiltonian ODE
systems.

The paper is organized as follows: Port-Hamiltonian ODE systems are briefly
introduced in Section 2. The following sections shows how the operator
splitting approach can be applied to port-Hamiltonian systems, with discrete
gradient schemes being the method-of-choice. Sections 5 and 6 show how this
approach can exploit both the multirate behaviour and system structure at the
level of nested integration and linear solvers. Numerical results for test
example from coupled multibody systems are discussed in Section 7. Finally, we
give a summary and outlock to future work.

##  2 Port-Hamiltonian ODE systems

An ODE-IVP of port-Hamiltonian structure vanderschaft2006  is given by

$\displaystyle\dot{x}$ $\displaystyle=(J(x)-R(x))\nabla H(x)+B(x)u(t),$ (7)
$\displaystyle y$ $\displaystyle=B(x)^{\top}\nabla H(x)$ (8)

with state variable $x:[t_{0},t_{end}]\rightarrow\mathbb{R}^{n}$, a twice
continuously differentiable Hamiltonian
$H:\mathbb{R}^{n}\rightarrow\mathbb{R}$, $n\times n$ matrices
$J(x)=-J(x)^{\top}$ skew-symmetric and $R(x)\geq 0$ positive semi-definite, an
input signal $u:[t_{0},t_{end}]\rightarrow\mathbb{R}^{d}$, an output
$y:[t_{0},t_{end}]\rightarrow\mathbb{R}^{d}$ and a $n\times d$-dimensional
port matrix $B(x)$.

A fundamental property for the solution $x(t)$ of a port-Hamiltonian system is
the dissipativity inequality

$\displaystyle\frac{d}{dt}H(x(t))$ $\displaystyle=\nabla
H(x(t))^{\top}\dot{x}(t)=-\nabla H(x(t))^{\top}R(x(t))\nabla
H(x(t))+y(t)^{\top}u(t)$ $\displaystyle\leq y(t)^{\top}u(t).$

This can also be written in the integral form

$\displaystyle H(x(t_{0}+h))$ $\displaystyle-H(x(t_{0}))=$
$\displaystyle=-\int_{t_{0}}^{t_{0}+h}\nabla H(x(s))^{\top}R(x(s))\nabla
H(x(s))\,ds+\int_{t_{0}}^{t_{0}+h}y(s)^{\top}u(s)\,ds$ (9)
$\displaystyle\leq\int_{t_{0}}^{t_{0}+h}y(s)^{\top}u(s)\,ds$

Without input and dissipation $R\equiv 0$, the Hamiltonian is an invariant of
the port-Hamiltonian ODE (7).

If port-Hamiltonian systems are solved numerically, the aim is not only
efficiency, but also structure preservation. We aim for schemes that preserve
the dissipativity inequality (2) at a discrete level.

##  3 Operator splitting for port-Hamiltonian systems

The splitting of the right-hand side of the port-Hamiltonian ODE (7) into an
energy-preserving part $f_{1}(x,t):=J(x)\nabla H(x)$ and a dissipative energy-
coupling part $f_{2}(x,t):=-R(x)\nabla H(x)+B(x)u(t)$ is a natural choice.

The approximation $z(h)$, given by the Strang splitting approximation

$\displaystyle\begin{pmatrix}\dot{v}\\\ \dot{s}_{v}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}-R(v)\nabla H(v)+B(v)u(s_{v})\\\
1\end{pmatrix},\qquad$ $\displaystyle\begin{pmatrix}v(t_{0})\\\
s_{v}(t_{0})\end{pmatrix}=\begin{pmatrix}x_{0}\\\ 0\end{pmatrix},$ (10)
$\displaystyle\begin{pmatrix}\dot{w}\\\ \dot{s}_{w}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}J(w)\nabla H(w)\\\ 0\end{pmatrix},\qquad$
$\displaystyle\begin{pmatrix}w(0)\\\
s_{w}(t_{0})\end{pmatrix}=\begin{pmatrix}v(h/2)\\\ h/2\end{pmatrix},$ (11)
$\displaystyle\begin{pmatrix}\dot{z}\\\ \dot{s}_{z}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}-R(z)\nabla H(z)+B(z)u(s_{z})\\\
1\end{pmatrix},\qquad$ $\displaystyle\begin{pmatrix}z(t_{0})\\\
s_{z}(t_{0})\end{pmatrix}=\begin{pmatrix}w(h)\\\ h\end{pmatrix}$ (12)
$\displaystyle y_{z}$ $\displaystyle=B(z)^{\top}\nabla H(z),$

provides an order two approximation $z(h/2)$ for the exact solution
$\varphi(h)\circ\begin{pmatrix}x_{0}\\\ 0\end{pmatrix}$of (7) at time point
$t_{0}+h$ and the output $y$, respectively, starting at the initial value
$x_{0}$ at $t_{0}$. In short hand:

$\displaystyle\begin{pmatrix}z(h/2)\\\ s_{z}(h/2)\end{pmatrix}$
$\displaystyle=$
$\displaystyle\varphi_{\tilde{f}_{1}}(\frac{h}{2})\circ\varphi_{\tilde{f}_{2}}(h)\circ\varphi_{\tilde{f}_{1}}(\frac{h}{2})\circ\begin{pmatrix}x_{0}\\\
0\end{pmatrix}=\varphi(h)\circ\begin{pmatrix}x_{0}\\\
0\end{pmatrix}+\begin{pmatrix}{\mathcal{O}}(h^{3})\\\ 0\end{pmatrix}$
$\displaystyle y_{z}(h/2)$ $\displaystyle=$ $\displaystyle
y(h)+{\mathcal{O}}(h^{3})$

The dissipativity inequality is fulfilled exactly at a discrete level, i.e.,

$\displaystyle H(z(h/2))-H(x_{0})$
$\displaystyle=(H(z(h/2))-H(z(0)))+(H(w(h))-H(w(0)))+$
$\displaystyle\quad+(H(v(h/2))-H(x_{0}))$
$\displaystyle\leq\int_{0}^{h/2}y_{v}(\tau)^{\top}u(\tau)\,d\tau+\int_{0}^{h/2}y_{z}(\tau)^{\top}u(\tau+h/2)\,d\tau.$

##  4 Discrete gradient methods

Applying operator splitting to the port-Hamiltonian system (7), (8) requires
the following demands:

  * •

the numerical approximations of (10) and (12) are dissipative at a discrete
level for vanishing $u\equiv 0$, i.e., it holds in this case

$\displaystyle H(v(h/2))\leq H(x_{0})\hskip 22.76228pt\text{and}\hskip
22.76228pt{H(z(h/2))\leq H(w(h))}$

holds for all step sizes $h>0$.

  * •

the numerical approximation $w(h)$ of (11) is energy preserving at a discrete
level, i.e.,

$\displaystyle{H(w(h))=H(v(h/2))}$

holds for all step sizes $h>0$.

Though (11) defines a symplectic flow for $J$ regular, symplectic schemes are
only energy preserving for quadratic Hamiltonians, but not in general. Hence
methods of choice for solving (7)-(8) here are discrete gradient methods
GM_Go96 , defined by

$\displaystyle\frac{x_{1}-x_{0}}{h}$
$\displaystyle=(\bar{J}(x_{0},x_{1},h)-\bar{R}(x_{0},x_{1},h))\bar{\nabla}H(x_{0},x_{1})+\bar{B}(x_{0},x_{1},h)\bar{u}(t_{0},t_{1},h),$
(13) $\displaystyle y_{1}$
$\displaystyle=\bar{B}(x_{0},x_{1},h)^{\top}\bar{\nabla}H(x_{0},x_{1}),$ (14)

with the discrete gradient
$\bar{\nabla}H:\mathbb{R}^{n}\times\mathbb{R}^{n}\rightarrow\mathbb{R}$
fulfilling the two conditions

$\displaystyle i)$ $\displaystyle\hskip
14.22636pt\bar{\nabla}H(x,y)^{\top}(y-x)=H(y)-H(x),$ $\displaystyle ii)$
$\displaystyle\hskip 14.22636pt\bar{\nabla}H(x,x)=\nabla H(x)$

for all $x,y\in\mathbb{R}^{n}$, and the skew-symmetric matrices
$\bar{J}(x,y,h)$, the positive-semidefinite matrices $\bar{R}(x,y,h)$, the
matrices $\bar{B}(x,y,h)$ and $\bar{u}(t_{0},t_{1},h)$ fulfilling the
compatability conditions
$\bar{J}(x,x,0)=J(x),\bar{R}(x,x,0)=R(x),\bar{B}(x,x,h)=B(x)$ and
$\bar{u}(t_{0},t_{0},0)=u(t_{0})$.

By construction, the discrete gradient method satisfies a discrete version of
the dissipativity inequality mclachlan1999geometric :

$\displaystyle H(x_{1})-H(x_{0})$ $\displaystyle=$
$\displaystyle-h\bar{\nabla}H(x_{0},x_{1},h)^{\top}\bar{R}(x_{0},x_{1},h)\bar{\nabla}H(x_{0},x_{1},h)+hy_{1}^{\top}\bar{u}(t_{0},t_{1},h).$

One notes that discrete gradient methods are implicit and limited to order two
for non-vanishing $R$, whereas order one would have been sufficient for
deriving an overall order two operator splitting scheme.

##  5 Multirate potential

A typical situation is that the energy preserving part defines a highly
oscillatory behaviour, which requires small stepsizes, whereas the dissipative
part may be characterized by a very slow dynamics, allowing for large step
sizes.

One idea to exploit this multrate behaviour is to use low (order two) schemes
for the dissipative part, as these are sufficient for obtaining accurate
results also for large stepsizes, but to use highly accurate numerical schemes
for the energy preserving part, which then also allows for larger step sizes.

Another way of exploiting this multirate behaviour given by different parts of
the right-hand side with a time constants ratio of $m\in\mathbb{N}$ is to
combine operator splitting with nested integration SEXTON1992665 , i.e., to
replace the step (4) with step size $h$ by $m$ subsequent steps of step size
$h/m$:

> $\dot{w}_{1}=J(w_{1})\nabla H(w_{1}),\quad w_{1}(0)=v(h/2)$
> for $i=2,\ldots,m:$ $\dot{w}_{i}=J(w)\nabla H(w),\qquad w_{i}(0)=w_{i-1}(h)$
> $w(h):=w_{m}(h)$

A straightforward way to derive higher order methods are composition methods.
Using the operator splitting approach based on the discrete gradient method as
the base scheme of order two, a composition of three base schemes with three
different step sizes yields an order four scheme. However, composition schemes
with order higher than two demand negative time steps suzuki1991 , which
contradicts the discrete dissipativity condition in the case of $R\neq 0$.

In lattice quantom chromodynamics, where the gauge action plays the role of
the fast and cheap part, and the fermionic action plays the role of the slow
and expensive part this problem could be circumvented by the following idea:
if one uses an integrator of second order for the slow action with step size
$h$, and approximates the fast action by $m$ steps of the second order scheme
with step size $h/(2m)$, the error of the overall multirate scheme will be of
order
${\mathcal{O}}(h^{2})+{\mathcal{O}}((\frac{h}{m})^{2})+{\mathcal{O}}(h^{4})$.
With the use of force gradient OMELYAN2003272  information only at the slowest
level it is possible to cancel the leading error term of order
${\mathcal{O}(h^{2})}$. If the mutirate factor $m$ between the time constants
of both schemes is high enough, one gets already for quite small step size
$h\leq\frac{1}{m}$ the overall order is then given by the leading error term
of order ${\mathcal{O}}(h^{4})$, i.e., the scheme has an effective order of
four GM_sh18 . Whether this approach can be applied to the port-Hamiltonian
setting is an open question.

Another idea is to use highly accurate numerical schemes for the energy
preserving part, which allows for larger step sizes.

##  6 Linear solvers

Operator splitting also allows for exploiting the special structure one
obtains when solving the energy-preserving and dissipative subsystems. For
simplicity of discussion, we assume constant matrices $J,R$ and $B$ and a
quadratic Hamiltonian $H(x)=x^{\top}Qx$ here. Note that the discrete gradient
method is equivalent to the implicit midpoint rule in case of linear systems.

Applying the standard average vector field method as a classical discrete
gradient method of order two, one yields the following linear system

$\displaystyle\frac{x_{1}-x_{0}}{h}=JQ\frac{x_{1}+x_{0}}{2}$
$\displaystyle\Rightarrow$
$\displaystyle\left(I-\frac{h}{2}\tilde{J}\right)\tilde{x}_{1}=\left(I+\frac{h}{2}\tilde{J}\right)\tilde{x}_{0}$
(15)

for the energy preserving part (11), where we have applied a congruency with
$Q^{1/2}$ with $\tilde{x}_{1}=Q^{1/2}x_{1}$, $\tilde{x}_{0}=Q^{1/2}x_{0}$ and
$\tilde{J}=Q^{1/2}JQ^{1/2}$.

Computing approximations for $\tilde{x}_{1}$ in (15) by applying a linear
solver to the matrix $I-\frac{h}{2}\tilde{J}$ is not a structure-aware
approach, since the computed approximations will, typically, not reflect
energy preservation, at least if, for efficiency reasons, we do not aim at a
very accurate solution. A structure aware approach to compute approximations
for (15) arises, if we use a matrix function approach, i.e. we consider

$\tilde{x}_{1}=C(\tilde{J},\tfrac{h}{2})\tilde{x_{0}},$

where $C(\tilde{J},\frac{h}{2})$ is the matrix function evaluation of the
Cayley transform

$C(z,a)=\frac{z+a}{z-a},\enspace z\in\mathbb{C}\setminus\\{a\\}.$

We can now use the Arnoldi method FrSi06  to evaluate the action of
$C(\tilde{J},\frac{h}{2})$ on the vector $\tilde{x}_{0}$. Each iteration of
this method requires one matrix-vector multiplication with $\tilde{J}$, and
unlike the general Arnoldi method it relies on a short recurrence because
$\tilde{J}$ is skew-hermitian. This approach is structure preserving, since
all iterates it produces will have the same 2-norm than $\tilde{x}_{0}$. This
is because the Arnoldi method obtains its $k$-th iterate as the Cayley
transform with a matrix which is the orthogonal projection of $\tilde{J}$ onto
the $k$-th Krylov subspace.

We can work in a similar manner with the original matrix $J$ by using the
inner product defined by $Q$ rather than the Euclidian inner product; see
ConGol76 ; FroKah2022 ; Gueetal2022 .

For the dissipative energy-coupling part (10) one gets the linear system

$\displaystyle\frac{x_{1}-x_{0}}{h/2}$ $\displaystyle=$ $\displaystyle-
RQ\frac{x_{1}+x_{0}}{2}+B\frac{u_{(}t_{0}+h/2)+u(t_{0})}{2}\Rightarrow$
$\displaystyle\left(I+\frac{h}{4}\tilde{R}\right)\tilde{x}_{1}$
$\displaystyle=$
$\displaystyle\left(I-\frac{h}{4}\tilde{R}\right)\tilde{x}_{0}+\tilde{B}\frac{u_{(}t_{0}+h/2)+u(t_{0})}{2},$

where we have applied again a congruency with $Q^{1/2}$ with
$\tilde{x}_{1}=Q^{1/2}x_{1}$, $\tilde{x}_{0}=Q^{1/2}x_{0}$ and
$\tilde{R}=Q^{1/2}RQ^{1/2}$, $\tilde{B}=Q^{1/2}B$. Now the iteration matrix
$I+\frac{h}{2}\tilde{R}$ is symmetric positive-definite. In a coupling
context, $R$ will be composed of the dissipation operators of the individual
systems, and we have the potential of using targeted preconditioners for each
of the systems to obtain particulalry efficient solvers for the composed
system, also for larger time steps. We do not go into further details here.

##  7 Numerical examples

Finally we discuss two numerical examples: a two mass oscillator with damping
GM_gu22b  for testing the different time integration approaches, and a single
mass-spring-damper (MSD) chain GPBS2012  available from the Port Hamiltonian
Benchmark
System111https://algopaul.github.io/PortHamiltonianBenchmarkSystems.jl to
validate the structure preserving properties of the the matrix Arnoldi
approach.

###  7.1 A two mass osciallator with damping

As an example to validate the different numerical time integration approaches
we consider two systems ($i=1,2$), each consisting of a mass $m_{i}>0$, which
is connected via a massless spring $K_{i}$ to walls with damping $r_{i}>0$,
which are coupled by a spring $K$, see Fig. 1.

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

To set up the coupled system, we have for the two systems the positions
$q_{1},\,q_{2}$ and momenta $p_{1},\,p_{2}$ for the masses $m_{1},\,m_{2}$,
respectively, as well as position $q=q_{2}$ as coupling variable between the
two systems. The port-Hamiltonian description of this coupled system (with the
coupling described by the off-diagonal elements of the skew-symmetric matrix
$J$) is given by

$\displaystyle\dot{x}$ $\displaystyle=(J-R)Qx,\qquad x(0)=x_{0}$

for the unknown $x:=(q_{1},q_{1}-q,q_{2},p_{1},p_{2})^{\top}$, with

$\displaystyle J=\begin{pmatrix}0&0&0&1&0\\\ 0&0&0&1&-1\\\ 0&0&0&0&1\\\
-1&-1&0&0&0\\\ 0&1&-1&0&0\end{pmatrix},\quad
R=\mbox{diag}(0,0,0,r_{1},r_{2}),\quad
Q=\mbox{diag}(K_{1}.K,K_{2},\frac{1}{m_{1}},\frac{1}{m_{2}}),$

and the Hamiltonian given by $H(x)=x^{\top}Qx$. With $\tilde{x}:=Q^{1/2}x$ one
gets

$\dot{\tilde{x}}=(\tilde{J}-\tilde{R})\tilde{x}$ (16)

with

$\displaystyle\tilde{J}=\begin{pmatrix}0&0&0&\sqrt{\frac{K_{1}}{m_{1}}}&0\\\
0&0&0&\sqrt{\frac{K}{m_{1}}}&-\sqrt{\frac{K}{m_{2}}}\\\
0&0&0&0&\sqrt{\frac{K_{2}}{m_{2}}}\\\
-\sqrt{\frac{K_{1}}{m_{1}}}&-\sqrt{\frac{K}{m_{1}}}&0&0&0\\\
0&\sqrt{\frac{K}{m_{2}}}&-\sqrt{\frac{K_{2}}{m_{2}}}&0&0\end{pmatrix},\quad
R=\mbox{diag}\left(0,0,0,\frac{r_{1}}{m_{1}},\frac{r_{2}}{m_{2}}\right),$

The waveforms of $x_{1},x_{2}$ and $x_{3}$ are depicted in Fig. 2.

Figure 2:  Exact solution of the first three components of $x$ for the
parameters $m_{1}=m_{2}=200,K=K_{1}=10,K_{2}=1000,r_{1}=5,r_{2}=2$

Fig. 3 shows the error plot for different operator splitting approaches
applied to the two mass oscillator. All schemes show nicely an order two
behaviour, they only differ in the error constant. Though both are of order
two. the exact operator splitting is about $10^{4}$ times more accurate than
operator splitting based on the discrete gradient method. The same applies to
the multirate approaches (large step size with highly accurate scheme vs.
nested integration with small time steps for the energy preserving part).

Figure 3:  Final error (Euclidean norm) for the exact operator splitting
solution (10)-(12) (x), operator splitting with discrete gradient method
(13)-(14) ($\cdot$), multirate approach with nested integration (o) and
multirate integration with higher order scheme ($\ast$).

###  7.2 A single mass-spring-damper chain

To highlight the structure preserving properties of the matrix Arnoldi
approach for system (15) we use the single MSD chain from the Port Hamiltonian
Benchmark
System222https://algopaul.github.io/PortHamiltonianBenchmarkSystems.jl with
vanishing inputs where we chose the size parameters such that we obtain a
dimension $N=10,000$. For this example, the largest eigenvalue of $J$ is $10i$
so that we chose $h=0.05\cdot\frac{1}{10}$ as a step size which should sample
the periods of all frequencies sufficiently well. The results are given in
Figure 4

Figure 4: Comparison of GMRES and the matrix Arnoldi method. Left: 2-norm of
the residual as a function of the iteration number. Right:
$|1-\|x^{(k)}\|/\|\tilde{x}_{0}|\|\,|$ as a measure for the deviation in
energy preservation.

The left part of Figure 4 reports the convergence history for two approaches
to solve (15). The first approach just uses the GMRES method for the matrix
$I-\frac{h}{2}\tilde{J}$, the second is the Arnoldi matrix function method. We
report the size of the residuals
$(I-\frac{h}{2}\tilde{J})x^{(k)}-(I+\frac{h}{2}\tilde{J})\tilde{x}_{0}$ for
the iterate $x^{(k)}$ for both methods. We see that when measuring accuracy
via the residual, the GMRES method is about 25% faster than the matrix Arnoldi
method. This corresponds to similar computational cost, since both methods
require one matrix-vector multiplication with $\tilde{J}$ per iteration. The
Arnoldi approach shows a decrease in quality in every other iteration. This
can be attributed to the fact that for $k$ odd the projected matrix cannot
well accomodate the symmetry of the spectrum of $\tilde{J}$ with respect to
the real axis.

The right part of Figure 4 shows that the matrix Arnoldi method does a perfect
job in energy preservation, since the 2-norm of all its iterates differ from
that of $\tilde{x}_{0}$ ,in a relative sense, by just machine precision. For
the GMRES approximations, the violation of energy preservation is quite
pronounced for the early iterates, and it becomes less as the iteration
proceeds.

Figure 5: Deterioration of energy conservation after several time steps. We
plot $|1-\|x_{i}\|/\|\tilde{x}_{0}\|\,|$, $x_{i}$ the numerical approximation
for $x(t_{0}+ih)$, as a measure for the deviation in energy preservation.

As a further illustration of the effects of the inherent structure
preservation of the matrix Arnoldi approach, Figure 5 reports a situation
occurring within a multi-rate setting, see Section 5. We now do 20 consecutive
steps of numerical integration, and in each step we stop the computation of
the numerical approximation to the solution of (15) once the residual is of
size $h^{2}$. This choice is motivated by the fact that the integration scheme
by itself has order 2. The figure again reports the quality of energy
preservation as the relative difference of the 2-norms of the computed
approximations for $x(t_{0}+ih),i=1,\ldots,20$ and the 2-norm of the initial
value. We see that the GMRES approach now presents a very severe violation of
energy preservation at later time points, whereas the matrix Arnoldi approach
again preserves energy perfectly up to machine precision. We note that we
required the same residual accuracy for both methods, which means that the
matrix Arnoldi approach takes about 25% more matrix vector products than
GMRES.

##  8 Summary

In this paper we have discussed how operator splitting methods at different
levels can be used for the numerical simulation of port-Hamiltonian systems
for both obtaining efficient and structure preserving methods: exact operator
splitting based on Strang splitting, the use of discrete gradient schemes,
exploiting the multirate behaviour in the splitting between structure-
preserving and dissipative part, and structure preserving numerical solution
of the respective linear systems by a matrix Arnoldi approach.

Open questions for future research comprise, amongst others, higher-order
schemes, tailored linear solvers for the dissipative part and generalization
to port-Hamiltonian DAE systems.

## References

  * (1) A. Bartel, M. Günther, B. Jacob, and T. Reis, Operator splitting based dynamic iteration for linear port-Hamiltonian systems, arXiv:2302.01195, (2022).
  * (2) P. Concus and G. H. Golub, A generalized conjugate gradient method for nonsymmetric systems of linear equations.  Comput. Meth. appl. Sci. Eng., 2nd int. Symp., Versailles 1975, Lect. Notes Econ. math. Syst. 134, 56-65 (1976)., 1976.
  * (3) M. Diab, A. Frommer, and K. Kahl, A flexible short recurrence Krylov subspace method for matrices arising in the time integration of port Hamiltonian systems and ODEs/DAEs with a dissipative Hamiltonian, arXiv:2205.13842, (2022).
  * (4) A. Frommer and V. Simoncini, Matrix Functions, in Model Order Reduction: Theory, Research Aspects and Applications, W. H. A. Schilders, H. A. van der Vorst, and J. Rommes, eds., Mathematics in Industry, Springer, Heidelberg, 2008, pp. 275–304.
  * (5) O. Gonzalez, Time integration and discrete Hamiltonian system, J. Nonlinear Sci., 6 (1996), pp. 449–467.
  * (6) C. Güdücü, J. Liesen, V. Mehrmann, and D. B. Szyld, On non-Hermitian positive (semi)definite linear algebraic systems arising from dissipative Hamiltonian DAEs, SIAM J. Sci. Comput., 44 (2022), pp. a2871–a2894.
  * (7) S. Gugercin, R. V. Polyuga, C. Beattie, and A. van der Schaft, Structure-preserving tangential interpolation for model reduction of port-Hamiltonian systems, Automatica J. IFAC, 48 (2012), pp. 1963–1974.
  * (8) R. I. McLachlan and G. R. W. Quispel, Splitting methods, Acta Numerica, 11 (2002), p. 341–434.
  * (9) R. I. McLachlan, G. R. W. Quispel, and N. Robidoux, Geometric integration using discrete gradients, Philos. Trans. R. Soc. Lond., B, Biol. Sci., 357 (1999), pp. 1021–1045.
  * (10) I. Omelyan, I. Mryglod, and R. Folk, Symplectic analytically integrable decomposition algorithms: classification, derivation, and application to molecular dynamics, quantum and celestial mechanics simulations, Computer Physics Communications, 151 (2003), pp. 272–314.
  * (11) J. Sexton and D. Weingarten, Hamiltonian evolution for the hybrid Monte Carlo algorithm, Nuclear Physics B, 380 (1992), pp. 665–677.
  * (12) D. Shcherbakov, M. Ehrhardt, J. Finkenrath, M. Günther, F. Knechtli, and M. Peardon, Adapted nested force-gradient integrators: The Schwinger model case, Communications in Computational Physics, 21 (2017), p. 1141–1153.
  * (13) G. Strang, On the construction and comparison of difference schemes, SIAM Journal on Numerical Analysis, 5 (1968), pp. 506–517.
  * (14) M. Suzuki, General theory of fractal path integrals with applications to many‐body theories and statistical physics, Journal of Mathematical Physics, 32 (1991), pp. 400–407.
  * (15) A. van der Schaft, Port-Hamiltonian systems: an introductory survey, in Proceedings of the International Congress of Mathematicians Vol. III, M. Sanz-Sole, J. Soria, J. Varona, and J. Verdera, eds., no. suppl 2, European Mathematical Society Publishing House (EMS Ph), 2006, pp. 1339–1365.  null ; Conference date: 22-08-2006 Through 30-08-2006.