A Generalized Scalar Auxiliary Variable Method for Structure-Preserving and Efficient Integration of Nonlinear port-Hamiltonian DAEs

Motivation

• Nonlinear DAE time discretization generally requires solving a nonlinear algebraic system at every time step.For large-scale problems, repeated residual and Jacobian evaluations, matrix assembly, and factorizations can dominate time-integration cost.Source: §1 · p. 1

Related approaches

• Existing SAV and GSAV approaches enable linearly implicit treatment of nonlinear contributions while targeting unconditional energy stability.Relaxed and energy-optimal corrections address the fact that auxiliary energy need not coincide with the true energy.Source: §1 · pp. 1–2

Gap

• A GSAV formulation for nonlinear pH-DAEs must incorporate external ports, energy-induced geometry, degenerate dissipation, and descriptor structure.The paper identifies this formulation as previously unavailable for the considered pH-DAE type.Source: §1 · p. 2

Contribution

• The proposed EOP-GSAV framework gathers state dependence in one explicitly evaluated term, leaving a constant implicit core and requiring one linear solve per time step.The constant matrix is invertible for every positive step size, so one factorization can serve the whole fixed-step integration.Source: §1 · p. 2
• The method uses an open-power passivity budget and provides order-k convergence for BDF-1 and BDF-2 in the index-one descriptor setting.The analysis also establishes consistency and positivity of the auxiliary variable.Source: §1 · p. 2

Port-Hamiltonian setting

• The pH-DAE class imposes skew-symmetric interconnection, positive semidefinite dissipation, regular index-one structure, input regularity, and Hamiltonian compatibility.The formulation also assumes consistent initial data, solution regularity, and a Hamiltonian bounded from below.Source: §2.1 · pp. 2–3

Constant-core splitting

• The discretization separates linear effort, structure, and dissipation components from a constant-core operator, gathering all state dependence into g(x).The constant operator is generally nonsymmetric, with Q-weighted accretivity replacing the positive self-adjoint structure of standard GSAV.Source: §2.1 · p. 3

Energy balance

• The Hamiltonian obeys an open power balance in which dissipation is offset by external port power.Because ports can inject or remove energy, the Hamiltonian need not decrease monotonically.Source: §2.1 · p. 4

Passivity projection

• The open-system EOP projection replaces a closed-system decay bound with a discrete passivity budget that accounts for port energy exchange.The resulting construction yields a discrete passivity inequality.Source: §2.2 · pp. 4–5

Temporal accuracy

• BDF-1 shows first-order convergence, while BDF-2 and implicit midpoint show second-order convergence in both state and Hamiltonian errors.The proposed schemes’ order-k convergence is also established rigorously.Source: §2.2 · p. 5

Computational properties

• The implementation uses a linear system with a reusable factorization for fixed order and step size, while explicitly evaluating physical dissipation and nonlinear terms.The extrapolation residual is of order O(δt^k) under the stated assumptions.Source: §2.2 · p. 5

• The projected GSAV scheme is analyzed for k=1 and k=2, proving order-k convergence, consistency of the relaxation scalar, and positivity and energy consistency of the projected auxiliary variable.Source: §3 · p. 5
• The index-one descriptor structure separates differential and algebraic components, with the algebraic variable recovered through a nonsingular algebraic block despite degeneracy of the energy semi-norm.Source: §3.1 · p. 6, §3.2 · p. 8
• Under the stated regularity and nonsingularity assumptions, confinement arguments keep pre-relaxation and projected states in bounded sets where nonlinear evaluations remain uniformly controlled.Source: §3.1 · pp. 6–7, §3.3 · p. 8
• For every δt > 0, the constant implicit matrix is invertible, so each fixed-order, fixed-step integration uses a time-independent factorization.Source: §3.1 · pp. 6–7
• Theorem 3.14 establishes convergence with consistent initialization, including the stated BDF-2 initial error bounds for the pre-relaxed state, projected state, Q-norm, relaxation scalar, and auxiliary energy.Source: §3.4 · p. 9
• The projected state satisfies the nonlinear algebraic constraint with the same order as the state approximation, while the EOP projection is energy-non-worsening and exactly tracks the Hamiltonian when the passivity budget is inactive.Source: §3.4 · pp. 11–12

Low-dimensional nonlinear stress test

• At δt = 0.1, BDF-2-EOP completes every stress-test step, while modified and simple Newton hit the 20-correction limit at t = 2.6.Full Newton converges for all steps at this coarsest tested step.Source: §4.2 · p. 14
• Under refinement, the frozen-iteration contraction quantity decreases to 0.11 at δt = 6.25 × 10−3, allowing simple Newton to recover despite Jacobian variation of approximately 23–27.The increasing E/δt contribution dominates the midpoint Jacobian.Source: §4.2 · p. 14
• BDF-2-EOP needs two factorizations versus 120 for full Newton in a near-matched state-accuracy stress-test comparison.The comparison is reported for median laptop runtimes.Source: §4.2 · p. 14

Accuracy and structure verification

• BDF-1-EOP and BDF-2-EOP have zero auxiliary passivity violation at every tested step size, while physical passivity is observed empirically for all three methods.Only the auxiliary inequality is guaranteed by the EOP construction.Source: §4.3.1 · p. 15

Time-integrator comparison

• BDF-2-EOP uses one linear solve per step after two startup/update factorizations, whereas midpoint variants require substantially more linear work from nonlinear iterations.This comparison uses N = 2048, δt = 5 × 10−3, and 1000 time steps.Source: §4.3.2 · p. 16
• BDF-2-EOP is faster throughout the tested work–precision range, although its advantage over modified and simple Newton narrows at tighter tolerances.The same trend holds for both state and Hamiltonian errors.Source: §4.3.2 · p. 16

Problem-size scaling

• BDF-2-EOP remains faster than all midpoint realizations as system size grows from 253 to 131 069 unknowns, with its advantage over simple Newton increasing.For N ≥8192, runtime grows approximately linearly with problem size and normalized cost is nearly constant across the three largest systems.Source: §4.3.2 · p. 16

Comparison with SUNDIALS IDA

• At equal order, the finest IDA comparison is nearly balanced: BDF-2-EOP takes 15.70 s versus 15.45 s for order-two IDA, while unrestricted IDA increasingly leads as accuracy tightens.Unrestricted IDA reaches ρt = 9.49 for state error and ρt = 10.86 for Hamiltonian error at the finest measured points.Source: §4.3.3 · pp. 18–19

Framework

• BDF-1 and BDF-2 use one linear solve per step with reusable factorizations while retaining order-k convergence and discrete auxiliary passivity.When the EOP budget is inactive, the auxiliary energy exactly tracks the shifted physical Hamiltonian.Source: §5 · p. 19

Numerical results

• In the nonlinear stress test, BDF-2-EOP remains robust where Jacobian-reuse midpoint iterations deteriorate or fail.The test isolates robustness under strongly varying nonlinear Jacobians.Source: §5 · p. 19
• On the large pH-DAE benchmark, BDF-2-EOP outperforms tested implicit midpoint realizations over the matched-accuracy range.It avoids repeated nonlinear solves and state-dependent factorizations.Source: §5 · p. 19

Comparison with IDA

• Unrestricted variable-order IDA becomes substantially faster at high accuracy, while order-two comparisons show comparable state work–precision behavior.The order-two diagnostic remains within a small constant factor in the Hamiltonian metric.Source: §5 · p. 19

Scope and outlook

• The fixed-operator formulation supports future higher-order, adaptive, broader-structure, distributed, and sparse or iterative-linear-algebra extensions.These directions are identified as future work rather than established capabilities of the present schemes.Source: §5 · p. 19

Physical model

• The system contains alternating heavy and light masses connected by N−1 parallel links with nonlinear FPU-β springs and Maxwell branches.Each Maxwell branch combines a linear spring with a dashpot, and massless internal nodes create algebraic constraints.Source: A Nonlinear mass-spring-damper… · pp. 19–20, Figure 3 · p. 20

Nonlinear structure

• A monotone generalized strain coordinate creates nonlinear interconnection through a coordinate change without altering the underlying physical kinematics.The relations are ṙ = Γ(r)Dv and δ̇ = Dv.Source: A Nonlinear mass-spring-damper… · p. 20

State and Hamiltonian

• The state comprises momentum, generalized primary-spring strain, Maxwell-spring extension, and massless-node velocity.These variables define the benchmark’s physical and auxiliary components.Source: A Nonlinear mass-spring-damper… · p. 20

Dissipation

• State-dependent dissipation is exercised through a cubic relative-velocity damper alongside constant Maxwell dissipation.The dissipation is represented within the pH formulation through R(x) = R0 + R1(x).Source: A Nonlinear mass-spring-damper… · p. 21

pH-DAE formulation

• The resulting descriptor system is index one and passive, with the algebraic equation uniquely solvable under positive-definite Maxwell damping.The nonlinear pH-DAE is written using a constant-core decomposition plus an explicitly evaluated nonlinear term.Source: A Nonlinear mass-spring-damper… · p. 21

Scaling and forcing

• The benchmark scales dynamically by fixing spatial wavelengths so that increasing N adds active nonlinear components rather than quiescent degrees of freedom.The nonlinear algebraic component vanishes identically for this benchmark.Source: A Nonlinear mass-spring-damper… · p. 22