The Hamiltonian equations of motion for a system with N particles consist of 2 N first-order equations ∂ t q n = ∂ p n H , (1a) ∂ t p n = − ∂ q n H , (1b) where ∂ x is the partial derivative with respect to a general variable x . Each particle in the system has positions q 1 , q 2 , … , q N and momenta p 1 , p 2 , … , p n (Goldstein, 1980). H defines the Hamiltonian—the total energy of the system H p , q = ∑ n T p n + V q n , (2) where T is the component defining kinetic energy and V is the component defining potential energy. These components are dependent on the momentum and position of each particle. The system described by the Hamiltonian equations is not general. Notably, damping is omitted. This motivates the use of the more general PHS framework.

Port-Hamiltonian systems describe the interconnection of energy storage elements, energy dissipating elements, and power conserving elements using general power-conjugate variables of effort, e , and flow, f (Duindam et al., 2009). The product of effort and flow variables is equal to the power of the system, that is, the change in the total energy over time ( ∂ t H ) . Effort and flow can be given physical meaning in electrical systems as voltage and current (Section 4), in mechanical systems as force and velocity (Section 5), and in acoustic systems as pressure and volume velocity (Section 6). The effort and flow variables are related to one another by a set of internal state variables which is the integral of the effort with respect to time. Taking the system in Equations 1a, 1b as an example, the state variables are the positions and momenta of each particle. The efforts are the gradient of the Hamiltonian with respect to the state variables and the flow is the time-derivative of each state variable. By the chain rule, the product of the effort and flows corresponds to power.

A useful representation of a dynamical PHS with state x , input u , and output y is the input-state-output port-Hamiltonian system with direct feed-through (Duindam et al., 2009) ∂ t x ⏟ f x = J − R ⏟ A ∇ H x ⏟ e x + G − P ⏟ B u , (3a) y = G + P T ⏟ C ∇ H x + M + S ⏟ D u . (3b) e and f are the state dependent efforts and flows. ∇ is the gradient operator and J is a skew-symmetric matrix ( J = − J T ) representing the energy-preserving interconnection between storage components. R is a matrix representing dissipative components and must be positive-definite in a passive system. G encapsulates the relationship between external ports and the internal states with P representing the dissipation between external ports and the storage components. M , also a skew-symmetric matrix, is the energy-preserving interconnection and S the dissipation between the input and output ports. The resulting PHS structure consists of a mixed set of differential and algebraic equations (DAE) (Duindam et al., 2009). The system in Equations 3a, 3b will also be familiar as the state-space equations from dynamical systems theory and signal processing (Kailath, 1980; Smith, 2010). Equation 3a is the state update equation and Equation 3b is the output equation. If damping, input, and output variables are omitted, and we define x = q p T , f ( x ) = ∂ t x , and e ( x ) = ∇ H ( x ) then Equation 3a is equivalent to Equation 1b.

Defining the total effort E = ∇ H x u , and flow variables F = ∂ t x − y , we can rewrite Equation 3 as F = J G − G T − M ⏟ J E − R P P T S ⏟ R E , (4) where we have the block skew-symmetric matrix J and the block dissipative matrix R . Due to the property for any skew-symmetric matrix x T J x = 0 , the inner product of efforts and flows expresses the power balance E , F = ∇ H T x ∂ t x ⏟ = ∂ t H Energy variation = u T y ⏟ B External power − E T R P P T S E ⏟ Q > 0 Dissipated power . (5)

Equation 5 states that the variation in the total stored energy is equal to the energy flowing in and out of the system, B , minus the energy that is dissipated over the dissipative elements, Q . The system is passive if the change in stored energy is less than the energy flowing into and out of the system ∂ t H ≤ y T u . This holds true if the dissipative elements satisfy R = R P P T S ≥ 0 . (6)

We will also define a conjugate PHS where the variables ∂ t x and ∇ H x are interchanged. The total effort will then be E = ∂ t x u and the total flow will be F = ∇ H x − y . The conjugate can be seen as a particular case of a hybrid PHS (van der Schaft and Jeltsema, 2014, Chapter 5) where all the storage efforts and flows have been interchanged. This form will be seen in Section 6.2. A hybrid PHS, where one set of efforts and flows have been interchanged, is used to model the switching PHS tonehole in Section 7. These alternative representations still maintain the power balance in Equation 5 as each row in the total effort and flow corresponding to storage elements is still a set of power-conjugated variables. Hybrid Dirac structures are used in Müller and Hélie (2021) along with dedicated solvers based on projection methods and a Newton iteration.

An important class of Hamiltonian functions are quadratic Hamiltonians. Quadratic Hamiltonian functions appear frequently in the systems we are interested in and have the following form and gradient, H x = 1 2 x T L x ,   ∇ H x = L x , (7) where L is a symmetric matrix. Passivity is guaranteed if L ≥ 0 . Even if the Hamiltonian is not quadratic it can be turned into a quadratic function by a change of variables using energy quadratization methods (Lopes et al., 2015; Yang, 2016).

We define the following general shift operators applied to a function of time, t , dimension z , and a general variable s e s + u t , z , s = u t , z , s + Δ s = u l n s + Δ s , e s − u t , z , s = u t , z , s − Δ s = u l n s + Δ s . (8)

Because FDTD methods represent variables at regular discrete time, t = n Δ t , and spatial coordinates, z = l Δ z , a shorthand is used where the temporal and spatial coordinates are absorbed into a superscript and subscript, respectively. Here, we describe these operators on a general variable s to highlight that the operators are not limited to just time and spatial dimensions.

Elemental finite difference operators, represented by the operator δ approximate continuous-time partial derivatives ∂ s . Some elementary operators are the forward ( δ s + ) , backward ( δ s − ) , and centered ( δ s ⋅ ) difference operators ∂ s ≈ δ s + = 1 Δ s e s + − 1 δ s − = 1 Δ s 1 − e s − δ s ⋅ = 1 2 Δ s e s + − e s − . (9) Of use are the forward ( μ s + ) , backward ( μ s − ) and centered ( μ s ⋅ ) , which approximate unity, 1 ≈ μ s + = 1 2 e s + + 1 μ s − = 1 2 1 + e s − μ s ⋅ = 1 2 e s + + e s − . (10)

These operators are useful for collocating variables in discrete time. Higher-order derivatives can be approximated by combining these elementary operators.

The discrete gradient method is an energy-conserving discretization method for Hamiltonian systems. It involves the approximation of the gradient operator by a forward difference with respect to the state ∇ H x n ≈ δ x + H x = H x n + Δ x − H x n ⊘ Δ x , (19) where δ x + is the forward difference with respect to the state vector x and ⊘ defines the Hadamard or element-wise division between vectors. The key to the discrete gradient method is defining Δ x as the difference between the state at successive time steps Δ x = x n + 1 − x n . (20)

Defining the inner product as x , y = x T y for two vectors x and y ∈ R N , the discrete variation in energy is equivalent to the inner product of the discrete gradient and the forward difference of the state in time δ t + h = δ x + H x , δ t + x , (21) and the discrete gradient method observes a discrete chain rule. The presented definition of the discrete gradient is only one particular definition and other definitions exist (Yalçin et al., 2015).

For a quadratic Hamiltonian, the discrete gradient approximates the continuous-time gradient in Equation 7 via the midpoint rule. This corresponds to the forward average operator δ x + H x n = 1 2 L x n + 1 + x n = μ t + L x n . (22)

Based on Equation 11, the discrete energy is conserved as δ t + x n , μ t + , L x n = δ t + 1 2 x n T L x n . (23)

Any PHS discretized with the discrete gradient method produces an implicit and unconditionally stable scheme. However, if the Hamiltonian is quadratic and the resulting scheme is linearly implicit, an explicit update form can be derived via a linear system solution. To derive the scheme approximating Equations 3a, 3b, we approximate ∂ t with the forward difference in time Equation 9 and aprapproximate ∇ H as described in Equation 22. The resulting linearly implicit scheme has the following explicit update and output equations x n + 1 = A d x n + B d u n , (24a) y n = C d x n + D d u n , (24b)

with H = I − Δ t 2 A L − 1 ,   A d = H I + Δ t 2 A L ,   B d = Δ t ⋅ H B , (25a) C d = 1 2 C L A d + I ,   D d = 1 2 C L B d + D , (25b)

defined based on the matrices in Equations 3a, 3b. Unless the system parameters are time-varying, the discrete state-space matrices (with subscript d ) can be computed and stored prior to the simulation. In the case that we have a conjugate PHS system where ∂ t x and ∇ H are interchanged, the discrete gradient method results in the following scheme L μ t + x n = A δ t + x n + B u n , (26a) y n = C δ t + x n + D u n , (26b) and the explicit update matrices are then H = A − Δ t 2 L − 1 ,   A d = H A + Δ t 2 L ,   B d = − Δ t ⋅ H B , (27a) C d = 1 Δ t C A d − I ,   D d = 1 Δ t C B d + D . (27b)

Symplectic methods such as the Störmer-Verlet method are specifically used to discretize partitioned systems such as those describing Hamiltonian dynamics. Consider a general non-autonomous system defined by two variables ( q , p ) and the input variable u . q and p are both dependent on a variable x and related to one another by functions f , g , and the partial derivative with respect to x ∂ x p = f q , p , u , ∂ x q = g q , p , u . (28)

Normally p and q correspond to momentum and position, and ∂ t q = p / m . This is not necessarily the case in the general system. Hairer et al. (2003) describe a Störmer-Verlet scheme for discretizing general autonomous partitioned systems. We propose an extension of the method to non-autonomous systems δ x − p x + 1 2 Δ x = 1 2 f q x , p x + 1 2 Δ x , u x + f q x , p x − 1 2 Δ x , u x , (29a) δ x + q x = 1 2 g q x + Δ x , p x + 1 2 Δ x , u x + 1 2 Δ x + g q x , p x + 1 2 Δ x , u x + 1 2 Δ x . (29b)

For linearly damped autonomous systems, the method in Equations 29a, 29b results in a contraction of symplectic area. In Chatziioannou (2019) it was found for a mass-spring-damper system that the contraction of the Störmer-Verlet scheme is a Padé approximation to the continuous-time contraction. Determining if this property holds for the class of all linear non-autonomous PHSs discretized with the proposed method and if the proposed method is symplectic in the non-autonomous case is outside the scope of this article.

The Störmer-Verlet method involves placing one variable, here p , on an interleaved grid at indices x + 1 2 Δ x and keeping the other variable, here q , on the normal grid. Positions between the interleaved and normal grid are related by the forward or backward averaging operators. For instance, in the case of the input variable u x = μ x − u x + 1 2 Δ x , u x + 1 2 Δ x = μ x + u x , (30)

Equation 29a can be understood to be centered on the normal grid. As such, the right-hand-side of Equation 29a averages f evaluated at p x + 1 2 Δ x and p x − 1 2 Δ x with q and u fixed at x . Conversely, Equation 29b is centered on the interleaved grid. The traditional three-step form of the Störmer-Verlet scheme is given for the general method as p x + 1 2 Δ x = p x + Δ x 2 f q x , p x + 1 2 Δ x , u x , (31a) q x + Δ x = q x + Δ x 2 g q x + Δ x , p x + 1 2 Δ x , u x + 1 2 Δ x + g q x , p x + 1 2 , u x + 1 2 Δ x , (31b) p x + Δ x = p x + 1 2 Δ x + Δ x 2 f q x + Δ x , p x + 1 2 Δ x , u x + Δ x . (31c)

The expression in Equations 29a, 29b can be retrieved by shifting Equation 31c back by Δ x and substituting the result into Equation 31a. Computationally, the form in Equations 29a, 29b is more efficient than Equations 31a, 31b, 31c (Hairer et al., 2003). If regular grid values are required, Equation 31c can be applied ad hoc. Unlike the non-general Störmer-Verlet scheme, the general scheme is not guaranteed to be explicit. If f and g are linear functions, then the resulting linearly implicit scheme can be resolved into an explicit update form.

When using the Störmer-Verlet method to discretize a PHS with a quadratic gradient, we no longer approximate the gradient. The output equation is untouched save for the addition of temporal indexes and averaging operators. Consider the quadratic PHS defined by a generally partitioned version of Equation 4 ∂ t p ∂ t q − y = J p K p q G p − K p q T J q G q − G p T − G q T − M − R p R p q P p R q p R q P q P p T P q T S L p p L q q u , (32) where J p and J q are skew-symmetric matrices. All resistive matrices are ≥ 0 and L p and L q are symmetric matrices arising from the quadratic Hamiltonian definition. A discretization with the Störmer-Verlet method produces the following system of equations. For brevity, we omit temporal indexes. p is located at times n + 1 2 and q and u are located at times n δ t − p = J p − R p L p μ t − p + K p q − R p q L q q + G p − P p u , (33a) δ t + q = J q − R q L q μ t + q + − K p q T − R q p L p p + G q − P q μ t + u , (33b) − y = − G p T − P p T μ t − L p p + − G q T − P q T L q q + − M − S u . (33c)

Any system discretized with the Störmer-Verlet method will be conditionally stable. We will not prove the general case and leave proofs for specific systems.

Energy analysis is an informative tool that can aid the understanding of any FDTD scheme. In particular, it can help one determine if a scheme is conditionally or unconditionally stable, and—if a scheme is conditionally stable—what are the necessary criteria for stability. We evaluate the energy balance using the accumulated energy error, h err , described in (Harrison-Harsley, 2018; Equation 3.75) h err n = h n + 1 − h 0 + Δ t ∑ m = 0 n q m − b m ⌊ h 0 ⌋ 2 , (34) where h , q and b are the discrete time correlates to the energetic quantities in Equation 5 and Equation 34 is a discrete-time integration of said equation. ⌊ ⋅ ⌋ 2 denotes the rounding down to the nearest power of two used to reduce variations to machine precision. The stability of a scheme is established by deriving the conditions where h n , q n ≥ 0 , ∀ n . In an unconditionally stable scheme this will always be true, but in conditionally stable schemes the limiting factors must be derived.

Frequency domain analysis is also a helpful tool and is equivalent to a z -Transform analysis (Bilbao, 2009). Frequency domain analysis is typically carried out using a frequency domain ansatz or single frequency solution ( u n = z = e j ω Δ t ) and is useful for evaluating the effects of numerical dispersion. Frequency domain analysis, in this manner, is limited to linear systems.

We now apply the numerical methods discussed in Section 3 to the RLC PHS in Equation 4. Two criteria that are compared here are numerical stability and numerical dispersion. The accuracy of the two discretization methods is shown in Figure 2 in comparison to the analytical admittance given in Equation 37. Additionally, the accumulated energy error Equation 34 and a comparison of numerical dispersion in the two schemes is shown.

A lumped reed model aims to capture the behavior at the tip of the reed. We model the reed as a mass-spring-damper system which is a mechanical analogy of the RLC circuit presented in Section 4. The system is based on the position y as seen in Figure 3. y = 0 is the equilibrium position of the reed once the player has positioned their embouchure (the position and mechanical properties of the player’s lips), y l is the displacement from equilibrium to the mouthpiece lay, and y c is the displacement at which the Hunt-Crossley contact force begins to take effect. The lumped reed system is governed by the following equation: m d t 2 y + m γ d t y + k y − f c h = S r p Δ , p Δ = p m − p i n . (52) m , γ , and k are the effective mass, damping, and stiffness of the reed, which are affected by the player’s embouchure. S r is the effective reed area, also determined by the player’s embouchure, and S r p Δ is the external force due to the pressure difference across the reed channel (van Walstijn and Avanzini, 2007). p m is the upstream pressure supplied by the player and p i n is the pressure at the interface with the instrument bore. f c is the Hunt-Crossley contact force and is a function of compression: h = y − y c . It is important to note that y c is not equivalent to y l , as shown in Figure 3. y c has been empirically found to be a factor around 0.4 − 0.6 times the displacement y l (Chatziioannou and van Walstijn, 2012; Chatziioannou et al., 2019). The contact force is defined as, f c h = − ∂ h V c h , t − γ c ∂ t V c h , t = − ∂ h V c 1 + γ c ∂ t h , (53) with γ c the contact damping and V c the contact potential energy, V c h , t = k c α + 1 h t + α + 1 . (54) k c represents the contact stiffness and α ≥ 1 is the contact power-law exponential. These parameters may be determined via the material parameters of the colliding objects (van Walstijn et al., 2024b) or through inverse modeling (Chatziioannou et al., 2019). The binary relationship between contact and no-contact is determined by the operator ⋅ + where h + = max 0 , h . (55)

The Hunt-Crossley model can be interpreted as a parallel combination of a nonlinear spring and damper that resists the incursion of one object into another.

A second nonlinearity within the reed model is the flow through the reed channel. We assume a steady, simple, and quasi-stationary Bernoulli flow model (Chaigne and Kergomard, 2016; Sec. 10.3.1.1) u f = sign p Δ S j y 2 p Δ ρ , S j y = w y l − y + . (56) S j ( y ) is the variable jet-area approximated by a rectangular surface with jet-width w . The steady Bernoulli flow model is found in Scavone (1997); van Walstijn and Avanzini (2007); Chatziioannou et al. (2019), among many others. It is likely pertinent to consider an unsteady flow model as was done by Lopes and Hélie (2016). However, we compensate for the simplified model by including a pumping flow, u r = S r d t y , (57) that considers the contribution of reed motion to the flow. It was observed in Lopes and Hélie (2016) that the inclusion of a pumping flow is necessary for the model to be physically realistic and energy conserving. We now present a complete PHS description of the reed system based on Equations 52, 53, 56, 57, and denote momentum ν , reserving p for pressure variables d t ν y h ⏟ f r = ( 0 − 1 − 1 1 0 0 1 0 0 ⏟ J r − m γ + γ c ∂ h V c 0 0 0 0 0 0 0 0 ⏟ R r ) ∂ ν T r = d t y ∂ y V r = k y ∂ h V c ⏟ e r + S r − S r 0 0 0 0 ⏟ G r p m p i n , (58a) u m u i n = G r T e r + u f − u f . (58b) u m and u i n are the upstream and downstream volume velocities associated with p m and p i n , respectively. The total stored energy of the system in Equations 58a, 58b is H r = ν 2 2 m ⏟ T r + k y 2 2 ⏟ V r + V c . (59)

It is worth noting that the flow out of the reed model u i n = − ( u r + u f ) is negative in contrast to earlier literature due to sign conventions in defining the PHS. This is rectified when coupling the system with the woodwind bore. The pumping flow produced by G r T is a direct result of the PHS framework as the pumping flow complements the external forces due to p Δ . It is not immediately clear that the system in Equation 58 maintains an energy balance due to the nonlinear contact damping and Bernoulli flow. Following the derivation in Bilbao et al. (2015b) and Chatziioannou et al. (2019), we can demonstrate passivity by taking the inner product of the flows and efforts in Equation 58a e r , f r = ∂ t H r = − m γ ∂ ν T 2 ⏟ Q r − γ c ∂ h V c ∂ ν T 2 ⏟ Q c + e r , G r p m p i n T , (60) and identifying Q r as the damping due to the lumped reed model and Q c as the damping due to the contact force. Q c is non-negative as ∂ h V c is non-negative by definition. To simplify the last product in Equation 60, the inner product of Equation 58b is taken with p m p i n T resulting in the following energetic quantities e r , G r p m p i n T = p i n u i n + p m u m − p m u f + p i n u f , = p i n u i n + p m u m ⏟ B r − u f p Δ ⏟ Q f . (61) B r is the power flowing through the boundaries and Q f is the dissipation due to the Bernoulli flow through the reed channel equal to Q f = S j sign p Δ 2 ρ p Δ 1 / 2 ⋅ p Δ = S j 2 ρ p Δ 3 / 2 ≥ 0 . (62) The dissipation due to the Bernoulli flow is guaranteed to be non-negative by p Δ = sign p Δ p Δ and sign x 2 = 1 . Restating the energy balance of the system, ∂ t H r = B r − Q r + Q c + Q f ⏟ ≥ 0 , (63) the system in Equations 58a, 58b maintains an energy balance as all the dissipative terms and H r are all guaranteed to be non-negative.

We discretize the system in Equations 58a, 58b using the discrete gradient method and obtain Equation 64 by approximation of the time derivative and gradients. This results in a pair of implicit nonlinear equations in both the state-update equation (due to the Hunt-Crossley force) and the output equation (due to the Bernoulli flow). We will focus on the state-update equation and the discretization of the output equation will be left to Section 8. We are particularly concerned with the discrete-time transcendental equation: δ t + ν n = − k μ t + y n − δ h + V c n 1 + γ c μ t + ν n m − γ m μ t + ν n m + S r p Δ , (64) which corresponds to the discretized form of the first line in Equation 58a. It is helpful here to introduce an intermediary variable: ξ = h n + 1 − h n = y n + 1 − y n = Δ t δ t + y . (65)

We apply the property in Equation 13 to replace δ t + ν n and μ t + y n with factors of μ t + ν n and δ t + y n , respectively. Then, we utilize the discretized form of the second line in Equation 58a, δ t + y n = μ t + ν n m , (66) and Equation 65 to rewrite the transcendental equation as a function whose root is equal to ξ F ξ = ξ − ξ ∗ + 1 a 0 δ h + V c n 1 + γ c Δ t ξ = 0 , (67) with ξ ∗ = 1 a 0 S r p Δ + a 1 n , (68a) a 0 = 2 m Δ t 2 + k 2 + m γ Δ t , (68b) a 1 n = 2 Δ t ν n − k y n . (68c)

ξ ∗ represents the solution to ξ in the absence of the contact force. We approximate δ h + V c using the definition of the discrete gradient in Equation 19, F ξ = ξ − ξ ∗ + 1 a 0 V c ξ + h n − V c h n ξ 1 + γ c Δ t ξ . (69)

This equation can be solved using a Newton-Raphson iteration with a guaranteed solution due to the convexity of F ( ξ ) (Bilbao et al., 2015b). Division by ξ can prevent convergence and a more numerically robust definition of F ( ξ ) and ∂ ξ F ( ξ ) can be produced following the procedure in van Walstijn et al. (2024b). After ξ has been derived, the next states ν n + 1 , y n + 1 , and h n + 1 can be updated based on ξ via ν n + 1 = 2 m Δ t ξ − ν n , (70a) y n + 1 = ξ + y n , (70b) h n + 1 = ξ + h n . (70c)

The discrete Hamiltonian and dissipated energy terms associated with this scheme are h r = ν n 2 2 m + k y n 2 2 + h c , h c = V c h n , (71a) h r = m γ + γ c ∂ h V c h n μ t + ν n m 2 . (71b)

Explicit and linearly implicit schemes representing Hertzian contact dynamics have been derived using energy quadratization strategies in Ducceschi et al. (2021) and van Walstijn et al. (2024b). Explicit Hunt-Crossley contact dynamics were recently presented in van Walstijn et al. (2024a). Energy quadratization has also been used in the discretization of nonlinear port-Hamiltonian systems (Lopes et al., 2015). Invariant energy quadratization (IEQ) (Zhao et al., 2016) and scalar auxiliary variable (SAV) methods (Shen et al., 2018) apply a change of basis to the energy function by representing the energy with a quadratic variable. For single-variable problems, as treated here, the methods are equivalent. We will only quadratize the contact potential, leaving the other energies untouched. The contact potential is now defined based on an auxiliary state σ with energy V σ , V σ = 1 2 σ 2 , (72) and h is related to the new state by the function f ( ⋅ ) with σ = f ( h ) . This substitution is viable here as the given energy function is non-negative (Shen et al., 2018). ∂ h V c and ∂ t h are defined relative to the auxillary state by, ∂ h V c = g σ , ∂ t σ = g ∂ t h , (73) where g is the gradient variable equal to ∂ h f ( h ) , g = 1 2 k c α + 1 h + α − 1 . (74)