Source-linked AI summary
A Generalized Scalar Auxiliary Variable Method for Structure-Preserving and Efficient Integration of Nonlinear port-Hamiltonian DAEs
Aashutosh Sharma, Andreas Bartel, Manuel Schaller
TL;DR
Nonlinear pH-DAE integration must preserve energetic structure without making nonlinear solves prohibitively expensive. The paper develops an EOP-GSAV framework with constant-core BDF schemes, obtaining robust matched-accuracy performance while variable-order adaptive IDA remains faster at high accuracy.
Problem
Nonlinear pH-DAE time integration requires preserving energy and passivity while avoiding costly repeated nonlinear solves, and a GSAV algorithm for this setting was previously unavailable.
Method
The EOP-GSAV framework gathers state-dependent effort, interconnection, and dissipation terms into an explicitly evaluated contribution and applies relaxed BDF-1 and BDF-2 schemes with a passivity budget.
Results
The schemes require one linear solve per step with reusable factorizations, outperform tested implicit-midpoint realizations over matched accuracies, and show comparable work–precision to equal-order IDA.
Takeaways & Limitations
For the considered pH-DAE class, structure-preserving nonlinear integration can avoid repeated nonlinear solves while retaining order-k convergence and discrete auxiliary passivity.
Takeaways & Limitations
The analysis uses sufficient regularity assumptions and structurally requires the nondegeneracy condition Mdd > 0 for its energy argument.
Abstract
from arXiv · showhide
We develop an energy-optimal generalized scalar auxiliary variable (EOP-GSAV) framework for nonlinear index-one port-Hamiltonian differential-algebraic equations (pH-DAEs). Exploiting the port-Hamiltonian structure, we separate the nonlinear effort, interconnection, and dissipation terms from a constant implicit core. The resulting BDF-1 and BDF-2 schemes require one linear solve per time step with a reusable factorization, while retaining discrete passivity and accurate tracking of the Hamiltonian. The schemes are compared with the implicit midpoint method equipped with full, modified, and frozen-Jacobian Newton iterations. Numerical experiments ranging from a strongly state-dependent nonlinear stress test to large-scale benchmarks demonstrate robust and competitive performance, with substantial efficiency gains in matched-accuracy regimes. A comparison with SUNDIALS IDA shows comparable work-precision behavior at equal order despite a non-specialized Python/SciPy implementation, while unrestricted variable-order adaptive IDA is faster in the high-accuracy regime.
1 Introduction
The paper addresses the cost of preserving energetic structure when time-stepping nonlinear pH-DAEs, where nonlinear algebra can dominate integration. It develops an EOP-GSAV framework that gathers state dependence into an explicitly evaluated term while retaining passivity and efficient linear solves.
- 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.
- 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.
- 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.
- 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.
- Contribution: 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.
2 Adapting GSAV to pH-DAEs
The pH-DAE adaptation splits the model into a constant core and explicitly evaluated state-dependent terms, then applies relaxed IMEX BDF schemes with an open-system passivity projection. The resulting methods recover the observed first- and second-order temporal behavior while preserving the targeted auxiliary energy structure.
- 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.
- 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.
- 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.
- 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.
- 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.
- 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.
3 Analysis of the pH-DAE-GSAV scheme
The analysis establishes solvability, convergence, confinement, constraint accuracy, passivity, and energy fidelity for projected BDF-1 and BDF-2 EOP-GSAV discretizations of index-one pH-DAEs.
- 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.
- 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.
- Under the stated regularity and nonsingularity assumptions, confinement arguments keep pre-relaxation and projected states in bounded sets where nonlinear evaluations remain uniformly controlled.
- For every δt > 0, the constant implicit matrix is invertible, so each fixed-order, fixed-step integration uses a time-independent factorization.
- 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.
- 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.
4 Numerical results
The numerical studies test robustness, structure preservation, accuracy, scaling, and work–precision against implicit midpoint/Newton variants and SUNDIALS IDA. BDF-k-EOP remains robust and efficient across stress tests and large benchmarks, while unrestricted adaptive IDA leads at high accuracy.
- 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.
- Low-dimensional nonlinear stress test: 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.
- Low-dimensional nonlinear stress test: 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.
- 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.
- 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.
- Time-integrator comparison: 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.
- 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.
- 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.
5 Conclusions
The EOP-GSAV framework gives nonlinear index-one pH-DAEs fixed-core BDF schemes with reusable linear algebra, passivity, and Hamiltonian fidelity. Experiments show robustness and matched-accuracy advantages, while unrestricted variable-order IDA is faster at high accuracy.
- 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.
- 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.
- Numerical results: 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.
- 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.
- 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.
A Nonlinear mass-spring-damper benchmark system
The benchmark is a one-dimensional alternating-mass chain whose neighboring masses couple through nonlinear FPU-β springs and Maxwell branches. Generalized strain coordinates and state-dependent interconnection and dissipation produce an index-one, passive nonlinear pH-DAE.
- 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.
- 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.
- 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.
- 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).
- 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.
- 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.