A port‑Hamiltonian system [1,2] is defined by:
$$ \begin{aligned} \dot{\mathbf{x}} &= \bigl(\mathbf{J}(\mathbf{x}) - \mathbf{R}(\mathbf{x})\bigr) \nabla H(\mathbf{x}) + \mathbf{B}(\mathbf{x})\,\mathbf{u}_{\text{ext}}, \\ \mathbf{y} &= \mathbf{B}(\mathbf{x})^T \nabla H(\mathbf{x}) + \mathbf{D}(\mathbf{x})\,\mathbf{u}_{\text{ext}}. \end{aligned} $$The time derivative of the Hamiltonian satisfies:
$$ \frac{dH}{dt} = -(\nabla H)^T \mathbf{R} \,\nabla H \;+\; \mathbf{y}^T \mathbf{u}_{\text{ext}} \;\le\; \mathbf{y}^T \mathbf{u}_{\text{ext}}, $$which proves passivity: energy never increases beyond externally supplied power. This property is crucial for real‑time audio; it guarantees that no numerical instability can blow up the system.
Applying the implicit midpoint rule [2] yields a discrete system that preserves skew‑symmetry and passivity exactly for constant \(\mathbf{J}\). For state‑dependent \(\mathbf{J}(\mathbf{x})\) the Dirac structure is preserved to second order in the time step, and the energy balance \(dH/dt \le \mathbf{y}^T\mathbf{u}_{\text{ext}}\) remains satisfied because the discrete gradient is evaluated at the midpoint. This guarantees unconditional stability.
$$ \mathbf{x}^{n+1} = \mathbf{x}^n + \Delta t\,\bigl(\mathbf{J} - \mathbf{R}\bigr) \nabla H\!\left(\frac{\mathbf{x}^{n+1}+\mathbf{x}^n}{2}\right) + \Delta t\,\mathbf{B}\,\mathbf{u}_{\text{ext}}^{n+1/2}. $$This scheme is A‑stable, allowing large time steps without instability. The resulting sparse nonlinear system is solved via Newton‑Krylov iteration, exploiting the banded structure of the pipe equations.
The pipe is discretised with cell‑centred conserved variables: mass \(m_i = (\rho A \Delta x)_i\), momentum \(p_{\text{mom},i} = (\rho u A \Delta x)_i\), and total energy \(E_{\text{tot},i} = (\rho E A \Delta x)_i\). The Hamiltonian is \(H = \sum_i E_{\text{tot},i}\). The gradient \(\nabla H\) with respect to this state yields the intensive efforts (Chapter 3 shows how they are computed).
The discrete Dirac structure \(\mathbf{J}\) (skew‑symmetric) and dissipation matrix \(\mathbf{R}\) (positive semi‑definite) are derived from the staggered finite‑volume flux and the physical source terms. The construction follows the approach of Trenchant et al. (2018), where it is shown that the staggered scheme can be written in this form; explicit matrix blocks for a single cell are given in Chapter 3.