Pith. sign in

REVIEW 4 major objections 3 minor 70 references

Differential-algebraic equations can be differentiated end-to-end: jaxdae computes the gradient of a full multi-physics DAE solve, algebraic constraints and adaptive stepping included, under a single jax.grad call.

Reviewed by Pith at T0; open to challenge. T0 means a machine referee read the full paper against a public rubric. the ladder, T0–T4 →

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗

T0 review · deepseek-v4-flash

2026-08-01 03:17 UTC pith:7ZLLJOAB

load-bearing objection Real JAX ecosystem gap-filler with honest caveats; abstract overclaims the replay adjoint but the engineering and validation are sound. the 4 major comments →

arxiv 2607.23202 v1 pith:7ZLLJOAB submitted 2026-07-25 cs.MS cs.SYeess.SY

jaxdae: A JAX-native Differentiable Solver for Differential-Algebraic Equations in Coupled Multi-physics

Chengyuan Li , Shanfang Huang , Jian Deng This is my paper
classification cs.MS cs.SYeess.SY MSC 65L8065Y20
keywords Differential-algebraic equationsAutomatic differentiationJAXIndex reductionParameter inversionDifferentiable simulationFrozen-grid replay adjointBatch gradients
verification ladder T0 review T1 audit T2 compute T3 formal T4 reserved

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The paper aims to close a gap: engineered models that start as PDEs become differential-algebraic equations after spatial discretization, and parameter inversion, uncertainty quantification, Bayesian inference, and control all need gradients of the DAE solve. Industrial acausal modeling tools simulate DAEs forward but stop at reverse-mode differentiation, while differentiable-physics frameworks handle explicit ODEs and PDEs but ignore the algebraic-constraint layer. jaxdae claims to supply the missing overlay: a JAX-native DAE solver combining adaptive BDF, Radau, and Rosenbrock integration with Pantelides index reduction, dummy derivatives, and Modelica-style connector assembly, all differentiable through one jax.grad call. The key trick for the adaptive path is to freeze the accepted step grid and re-solve a variable-step BDF-2 on it for the backward pass, so the adaptive controller never enters the gradient. If correct, DAE solves become composable primitives for inference, control, and design, with batch-parallel gradients that the compiler fuses into a single program whose wall time stays nearly flat from batch size 1 to 1000.

Core claim

The paper's central discovery is that the forward DAE solve and its reverse-mode sensitivity can be unified in one JAX-native suite, and specifically that the adaptive BDF path can be made differentiable without trying to differentiate the adaptive controller. The method records the accepted step grid during a normal adaptive BDF(1–5) forward solve, stops the gradient on the grid, and in the backward pass re-solves a variable-step BDF-2 problem on the frozen nodes using an implicit-function-theorem adjoint of the stage residual. This 'frozen-grid replay adjoint' follows the controller-decoupling idea in the literature; the contribution claimed here is the first correct, composable JAX realiz

What carries the argument

The central mechanism is the frozen-grid replay adjoint: the forward solve runs adaptive BDF (orders 1–5) and records the accepted step grid; the backward pass applies stop_gradient to the grid and re-solves a variable-step BDF-2 on those fixed nodes, then differentiates the stage residual via the implicit function theorem, accumulating the adjoint from final time back to initial time. This decouples the controller from cotangent propagation, because step accept/reject branches would otherwise break the chain rule and yield silently wrong gradients. Supporting machinery includes Pantelides index reduction with dummy derivatives to reduce high-index DAEs to index-1 before integration, Jacobia

Load-bearing premise

The load-bearing premise is that the gradient jax.grad returns is the exact discrete adjoint of the variable-step BDF-2 replay operator on the frozen grid, and that this gradient tracks the sensitivity of the forward adaptive BDF(1–5) solve through the conditional O(tol^{2/3}) scaling argument; if that scaling fails — for tolerances below the floor, or for grid-dependent events — the returned gradient is not the derivative of the solve the user actually ran.

What would settle it

Integrate the scalar decay y' = -ky with adaptive BDF at rtols from 1e-4 down to 1e-12, compute dL/dk for L = y(1)^2 (analytic truth -2e^{-2}) via jax.grad, and fit the relative error slope versus rtol. If the fitted slope does not approach about 2/3 before the tolerance floor (and instead plateaus or worsens), the replay adjoint is not delivering the claimed tolerance-consistent gradient. A second check: place a state-dependent event whose detection time depends on the step grid, and observe whether the gradient error fails to converge — directly testing the frozen-grid assumption (iii).

Watch this falsifier — get emailed when new claim-graph text bears on it.

If this is right

  • Gradients of assembled multi-physics DAE models become obtainable in one jax.grad call, without hand-derived adjoints, checkpointing management, or finite-difference surrogates.
  • Bayesian inference, calibration, and optimal control over plant-scale DAEs can reuse the same compiled artifact as the forward solve, since the whole pipeline is a pure differentiable function.
  • vmap(jax.grad) yields batch-invariant wall time for parameter sweeps: the measured median is effectively flat from batch 1 to 1000, with the compiler fusing the batch into a single kernel.
  • The event-handling path recovers the Galán–Feehery–Barton jump-condition adjoint by composing the implicit-function-theorem event-time sensitivity with the jump map, verified term by term in the paper.
  • Under the conditional scaling argument, the replay gradient error is predicted to decay like O(tol^{2/3}) for BDF-2, and the paper shows raw adaptive-controller AD diverges while replay converges on the tested DAEs.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • Because the replay is fixed at BDF-2 while the forward path is BDF(1–5), a natural extension is to replay with BDF-4 or BDF-5 on the same frozen grid; the paper says this is straightforward, and doing so would sharpen whether the O(tol^{2/3}) rate improves with replay order.
  • The batch-flat scaling was measured on small per-instance systems (n=2 and n=20); extrapolating to plant-scale models with hundreds of variables is plausible but untested here, and the paper itself notes per-instance work will eventually dominate.
  • The same decoupling idea — freezing discrete decisions and replaying a simpler operator for the adjoint — could plausibly extend to other adaptive implicit methods, such as Radau, which the paper currently ships as forward-only; that would be a testable generalization of the replay approach.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit.

Referee Report

4 major / 3 minor

Summary. The paper introduces jaxdae, a JAX-native DAE solver suite that combines adaptive BDF(1–5), Radau IIA and fixed-step Rosenbrock integration with Pantelides/dummy-derivative index reduction and Modelica-style acausal assembly. Its headline contribution is differentiability: adaptive BDF solves are made reverse-mode differentiable by a frozen-grid 'replay adjoint' that stops gradients on the accepted step grid and re-solves a variable-step BDF-2 on that grid with an implicit-function-theorem adjoint; fixed-step BDF and ROS2 paths differentiate by full-function AD. The authors validate forward accuracy against analytic order tests and cross-library work–precision comparisons, AD against finite differences (empirical order 1.994, floor 1.4e-8), forward-mode against reverse-mode (6.4e-14), and the GFB event adjoint term by term. They report batch-invariant GPU throughput via jit(vmap(grad)) from batch 1 to 1000, and claim to be the first JAX-native DAE suite combining these capabilities. The paper is careful to state that the replay adjoint is an engineering realization, not a new convergence theorem, and explicitly lists limitations, including the conditional nature of the convergence heuristic.

Significance. If the headline claim held unconditionally, jaxdae would fill a real gap and be broadly useful in scientific machine learning. The paper has concrete strengths: a large and reproducible test suite, a detailed reproducibility manifest, careful AD-versus-FD consistency for the non-adaptive fixed-step path, a thorough cross-library work–precision study, a clean term-by-term verification of the event-adjoint channels, and unusually candid statements of limitations. However, the central differentiability claim for the adaptive path is only conditionally supported. The gradient returned for an adaptive BDF solve is the discrete adjoint of a BDF-2 replay operator, not of the variable-order forward solve, and the paper's own convergence support is a conditional scaling heuristic whose empirical fits are partially self-referential (references generated by the same solver family on different grids) and whose single analytic-truth fit excludes the rtol floor where error rebounds. Thus the abstract's one-line 'jax.grad' claim is stronger than what is demonstrated. The contribution remains potentially significant after either restricting the claim to the fixed-order BDF-2/scan path or supplying

major comments (4)
  1. [§3.7, assumption (iv), Eq. (16); Fig. 13 and Table 7] The load-bearing claim is that the frozen-grid replay gradient converges as O(tol^{2/3}) to the gradient of the adaptive solve. The conditional argument fixes p=2 in h_max ∝ tol^{1/(p+1)}, but the forward path is variable-order BDF(1–5). For a step accepted at forward order q>2, h scales as tol^{1/(q+1)}, so the BDF-2 replay global error is h^2 = tol^{2/(q+1)}; for q=5 this is tol^{1/3}, not tol^{2/3}. The numerical evidence does not disambiguate: fitted orders 0.07–0.38 against FD references are stated to be floored by reference error, and the analytic-truth panel (fitted order 0.83) excludes the rtol=1e-8 ADAPTIVE_RTOL_FLOOR point where the error rebounds. The abstract's statement that jax.grad returns the derivative of the adaptive solve is therefore stronger than demonstrated. A concrete repair is to restrict the differentiable adaptive path to BDF-2 so that the replay order matches
  2. [§7.2 vs §3.7] The paper contains a direct inconsistency in the semantics of the returned gradient. §7.2 ('Reverse-mode AD through implicit DAE solves') states that the solver 'registers the solve as differentiable' and that the 'backward pass differentiates exactly what the forward pass computed', calling this full-function AD. §3.7 states the opposite: 'the replay gradient is the exact discrete adjoint of the BDF-2 replay operator, not of the forward adaptive BDF(1–5) operator.' Since §7.5 makes the replay adjoint the default for the adaptive solver, a user calling jax.grad on an adaptive solve receives the replay gradient, not the exact full-function AD gradient promised by §7.2. The manuscript must either distinguish clearly between the scan/BDF-2 path and the adaptive BDF path, or change the §7.2 wording so that the default semantics are unambiguous.
  3. [§7.5, Table 7 and Fig. 13] The replay-adjoint validation is partly self-referential. The 'independently-discretised FD reference' is, per the text, a high-resolution fixed-step BDF-2 scan-adjoint gradient on a different grid, and the central-FD cross-check is performed against that scan-adjoint reference rather than against an FD gradient of the original adaptive BDF solve. This does not provide an independent estimate of the mismatch between the replay gradient and the derivative of the variable-order adaptive trajectory. The only grid-independent truth is the scalar decay, which is one-dimensional and excludes the rtol=1e-8 floor point. Please add a genuinely independent reference (e.g., a high-order Radau or RK adjoint on a fine grid, or analytic sensitivities for a small DAE family) and report the floor point explicitly instead of excluding it from the fit.
  4. [Abstract, §3.7 and §9] The abstract and conclusions frame jaxdae as making 'the adaptive BDF path differentiable' and 'the DAE becomes a differentiable primitive', while §3.7 and §9 explicitly disclaim that the replay is not the forward operator's adjoint and that its convergence is only a conditional heuristic. This tension needs to be resolved in the user-facing claims. As written, a reader is entitled to believe that jax.grad on an adaptive solve returns the derivative of the solve that was actually run; the manuscript itself says this is not what the replay delivers. Either restrict the headline claim to the fixed-order BDF-2/scan path where exactness is verified, or qualify the abstract so that the adaptive-path gradient is explicitly described as the BDF-2 replay gradient with conditional convergence.
minor comments (3)
  1. [Fig. 13 and Table 7] The legend writes 'O(rtol^p)' while the text uses the same symbol for the fitted order and for the replay order p. This is confusing; use e.g. 'fit exponent' and 'replay order' separately.
  2. [Table 5 and Fig. 9] Table 5 uses SUNDIALS IDA via scikit-sundae, while Fig. 9 uses Julia Sundials.jl. The text notes this, but the two IDA rows still invite cell-by-cell comparison. Add a footnote in both captions making the harness difference impossible to miss.
  3. [Program Summary / §5.1] The test suite reports 12 known failures reproducing identically at v0.27.1. For a software paper, it would be helpful to include a short appendix table listing these failures and their causes, rather than pointing only to RELEASE_NOTES.md.

Circularity Check

0 steps flagged

No significant circularity: the replay-adjoint claim is explicitly non-self-definitional and is anchored by analytic and finite-difference checks.

full rationale

The paper's central claims are software capabilities, not derived physical predictions. The frozen-grid replay adjoint is explicitly described as the exact discrete adjoint of the BDF-2 replay operator, not of the forward adaptive BDF(1-5) operator, and the O(tol^{2/3}) rate is presented as a conditional heuristic with stated assumptions rather than as a theorem. Gradient correctness is checked against a grid-independent analytic scalar-decay truth (fitted order 0.83), central finite differences (order 1.994 to a 1.41e-8 floor on the cart-pole; 1.8e-9 to 1.1e-7 on the other testbeds), and cross-library/benchmark references. The replay-versus-scan-adjoint comparisons on different grids of the same solver family are partly self-referential, but the paper identifies the reference floor explicitly and does not use those fits to define the predicted rate. The adoption of the controller-decoupling idea from Alexe and Sandu [54] is external prior work credited as such ('We do not claim a new method'), and the only author self-citation [69] supports an illustrative PWR testbed explicitly flagged as not a standardized validation benchmark. No fitted parameter is renamed as a prediction, and no load-bearing premise reduces to the paper's own prior claims.

Axiom & Free-Parameter Ledger

3 free parameters · 4 axioms · 0 invented entities

No new physical or mathematical entities are postulated; RuntimeParam and typed connectors are software constructs, not scientific objects. The central claim rests on the standard DAE input contract, the correctness of Pantelides/dummy-derivative reduction, and the conditional assumptions behind the replay-adjoint convergence rate.

free parameters (3)
  • PI controller gains k_P, k_I = not reported
    Eq. (11): chosen per method; they shape the accepted step grid and therefore the replay-adjoint error, but are standard controller values.
  • _ADAPTIVE_RTOL_FLOOR clamp = rtol = 1e-8 (implied by error rebound at 1e-8)
    Section 3.7: below this floor the controller contract h_max ∝ tol^{1/(p+1)} is violated; the analytic-truth fit excludes that point. It is a hand-set tolerance floor that affects the convergence claim.
  • Fitted replay-error power-law slopes = 0.19, 0.38, 0.07, 0.83
    Fig. 13 and Table 7: power-law fits of gradient error versus rtol used to compare with the predicted 2/3; these are data fits, not model parameters, but the convergence-rate evidence is built on them.
axioms (4)
  • domain assumption Input DAE has semi-explicit mass-matrix form M(t,y) ydot = f, 0 = g with block-diagonal M and residual affine in ydot; fully implicit forms are rejected.
    Section 2.1: all integrators, index-reduction, and adjoint paths assume this contract.
  • domain assumption Pantelides plus dummy derivatives produces an index-1 system with nonsingular ∂g/∂z and the same solution set.
    Sections 2.2 and 3.5: invoked for high-index input; the paper itself notes mechanical systems over-expand, so this is not universal.
  • ad hoc to paper Frozen-grid replay-adjoint convergence assumptions (i)–(v): smoothness, zero-stable frozen step ratios, no events/branches, controller contract h_max ∝ tol^{1/(p+1)}, padding sanitization.
    Section 3.7: these are stated assumptions rather than theorems; (iv) is explicitly called the weakest link and is violated by the rtol floor.
  • domain assumption Event transversality: switching velocity σ = s_t + s_x xdot^- ≠ 0.
    Section 3.6: required for the IFT formula for dt*/dp; the paper limits event-adjoint claims to transversal events.

pith-pipeline@v1.3.0-alltime-deepseek · 36450 in / 12841 out tokens · 120105 ms · 2026-08-01T03:17:55.997918+00:00 · methodology

0 comments
read the original abstract

Many engineered models begin as partial differential equations. Spatial discretization converts them into ordinary differential equations coupled to algebraic constraints---conservation closures, constitutive laws, network topology---whose joint evolution is a differential-algebraic equation (DAE). Parameter inversion, uncertainty quantification, Bayesian inference, and optimal control all require gradients of this solve. The two software traditions that should supply them have not met: industrial acausal modeling tools simulate DAEs forward but stop at reverse-mode differentiation, while differentiable-physics frameworks in JAX handle explicit ODEs and PDEs but leave the algebraic-constraint layer untouched. Here we show that the forward DAE solve and its reverse-mode sensitivity can be unified in one JAX-native suite. jaxdae pairs adaptive BDF, Radau, and Rosenbrock integration with Pantelides index reduction and dummy derivatives, and makes the adaptive BDF path differentiable by freezing the accepted step grid and re-solving a variable-step BDF-2 on it for the backward pass. The full pipeline differentiates under one $\texttt{jax.grad}$ call, XLA fuses a batched parameter sweep into one program whose wall time stays nearly flat from batch~1 to~1000, and the DAE becomes a differentiable primitive for inference, control, and design.

Figures

Figures reproduced from arXiv: 2607.23202 by Chengyuan Li, Jian Deng, Shanfang Huang.

Figure 1
Figure 1. Figure 1: Main solve pipeline of jaxdae. A DAE problem is defined by acausal assembly or by a direct residual, checked [PITH_FULL_IMAGE:figures/full_fig_p007_1.png] view at source ↗
Figure 2
Figure 2. Figure 2: Term-by-term recovery of the GFB jump-condition adjoint by AD-chain composition through the IFT-differentiated [PITH_FULL_IMAGE:figures/full_fig_p011_2.png] view at source ↗
Figure 3
Figure 3. Figure 3: Code structure of jaxdae. The @connector decorator and the TypedComponent base class form an object-oriented assembly shell; System.compile() crosses into the functional core, emitting a pure residual F(t, y, y,˙ params) that the solve_dae family and the reverse-mode adjoint consume uniformly. Dashed arrows mark the tracing boundary: only the residual and the arrays to its right ever enter jax.jit or jax.g… view at source ↗
Figure 4
Figure 4. Figure 4: Capability map of jaxdae. The core solver (adaptive BDF, Radau5, fixed-step ROS2) sits at the center; around it the [PITH_FULL_IMAGE:figures/full_fig_p015_4.png] view at source ↗
Figure 5
Figure 5. Figure 5: The four core jaxdae API surfaces, shown as native Python code. (a) Acausal multi-physics assembly: [PITH_FULL_IMAGE:figures/full_fig_p016_5.png] view at source ↗
Figure 6
Figure 6. Figure 6: Order verification of the fixed-step BDF family (orders 1–5) on [PITH_FULL_IMAGE:figures/full_fig_p018_6.png] view at source ↗
Figure 7
Figure 7. Figure 7: Order verification of Rosenbrock–W (ROS2, left) and Radau IIA (Radau5, right) on [PITH_FULL_IMAGE:figures/full_fig_p018_7.png] view at source ↗
Figure 8
Figure 8. Figure 8: Three canonical stiff benchmarks solved with [PITH_FULL_IMAGE:figures/full_fig_p019_8.png] view at source ↗
Figure 9
Figure 9. Figure 9: Work–precision diagram over five benchmark problems (panels a–e), three solver libraries, and the tolerance ladder [PITH_FULL_IMAGE:figures/full_fig_p021_9.png] view at source ↗
Figure 10
Figure 10. Figure 10: Batched scalar-gradient throughput (index-1 DAE, [PITH_FULL_IMAGE:figures/full_fig_p024_10.png] view at source ↗
Figure 11
Figure 11. Figure 11: Batched Fisher-DOE sweep over candidate CSTR operating temperatures ( [PITH_FULL_IMAGE:figures/full_fig_p025_11.png] view at source ↗
Figure 12
Figure 12. Figure 12: AD-versus-FD gradient consistency, cart-pole ( [PITH_FULL_IMAGE:figures/full_fig_p027_12.png] view at source ↗
Figure 13
Figure 13. Figure 13: Replay-versus-raw gradient error versus solver tolerance, four panels. (a–c) Three DAEs benchmarked against [PITH_FULL_IMAGE:figures/full_fig_p028_13.png] view at source ↗
Figure 14
Figure 14. Figure 14: Gradient correctness on a stiff Robertson [PITH_FULL_IMAGE:figures/full_fig_p029_14.png] view at source ↗
Figure 15
Figure 15. Figure 15: Three properties of jaxdae in one figure. (a) Gradient reliability: only the jaxdae frozen-grid replay adjoint converges; [PITH_FULL_IMAGE:figures/full_fig_p030_15.png] view at source ↗
Figure 16
Figure 16. Figure 16: IEEE 14-bus transient-stability forward simulation assembled and integrated in jaxdae. Rotor-angle response of the [PITH_FULL_IMAGE:figures/full_fig_p032_16.png] view at source ↗

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

70 extracted references · 12 canonical work pages

  1. [1]

    K. E. Brenan, S. L. Campbell, L. R. Petzold, Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations, 2nd Edition, SIAM, 1996

    1996

  2. [2]

    U. M. Ascher, L. R. Petzold, Computer Methods for Ordinary Differential Equations and Differential- Algebraic Equations, SIAM, 1998

    1998

  3. [3]

    C. C. Pantelides, The consistent initialization of differential-algebraic systems, SIAM Journal on Sci- entific and Statistical Computing 9 (2) (1988) 213–231.doi:10.1137/0909014

    doi:10.1137/0909014 1988

  4. [4]

    S. E. Mattsson, G. Söderlind, Index reduction in differential-algebraic equations using dummy deriva- tives, SIAM Journal on Scientific Computing 14 (3) (1993) 677–692.doi:10.1137/0914043

    doi:10.1137/0914043 1993

  5. [5]

    Y. Cao, S. Li, L. Petzold, R. Serban, Adjoint sensitivity analysis for differential-algebraic equations: Algorithms and software, SIAM Journal on Scientific Computing 22 (5) (2003) 1816–1831.doi:10. 1137/S1064827501380630

    2003

  6. [6]

    R. J. LeVeque, Finite Volume Methods for Hyperbolic Problems, Cambridge University Press, 2002

    2002

  7. [7]

    Tarantola, Inverse Problem Theory and Methods for Model Parameter Estimation, Society for Industrial and Applied Mathematics (SIAM), 2005.doi:10.1137/1.9780898717921

    A. Tarantola, Inverse Problem Theory and Methods for Model Parameter Estimation, Society for Industrial and Applied Mathematics (SIAM), 2005.doi:10.1137/1.9780898717921

    doi:10.1137/1.9780898717921 2005

  8. [8]

    A. G. Baydin, B. A. Pearlmutter, A. A. Radul, J. M. Siskind, Automatic differentiation in machine learning: a survey, Journal of Machine Learning Research 18 (2018) 1–43

    2018

  9. [9]

    C. C. Margossian, A review of automatic differentiation and its efficient implementation, WIREs Data Mining and Knowledge Discovery 9 (4) (2019) e1305.doi:10.1002/widm.1305

    doi:10.1002/widm.1305 2019

  10. [10]

    S. S. Schoenholz, E. D. Cubuk, Jax-md: End-to-end differentiable, hardware accelerated, molecular dynamics in pure python, in: Advances in Neural Information Processing Systems (NeurIPS), 2020. arXiv:1912.04232

    Pith/arXiv arXiv 2020

  11. [11]

    Y. Hu, X. Liu, X. Zhang, Y. Fu, C. Kuang, W. Gan, N. Wan, H. Zhao, D. Xie, S.-M. Shan, et al., Difftaichi: Differentiableprogramming forphysical simulation, in: International ConferenceonLearning Representations (ICLR), 2020.arXiv:1910.00935

    Pith/arXiv arXiv 2020

  12. [12]

    B.Buhendwa, N

    D.A.Bezgin, A. B.Buhendwa, N. A.Adams, Jax-fluids: A fully-differentiable high-ordercomputational fluid dynamics solver for compressible two-phase flows, Computer Physics Communications 282 (2023) 108527.arXiv:2203.13760,doi:10.1016/j.cpc.2022.108527

    Pith/arXiv arXiv 2023

  13. [13]

    T.-K. Xue, H. Liao, P. Schmitz, S. Desai, A. Barmak, H. Gardin, W. Sun, G. Bauer, J. Zehnder, Q. Chen, et al., Jax-fem: A differentiable gpu-accelerated 3d finite element solver for automatically solving forward and inverse problems, Computer Physics Communications 291 (2023) 108802.arXiv: 2212.00964,doi:10.1016/j.cpc.2023.108802

    Pith/arXiv arXiv 2023

  14. [14]

    Ataei, H

    M. Ataei, H. Salehipour, XLB: A differentiable massively parallel lattice Boltzmann library in Python, Computer Physics Communications 300 (2024) 109187.doi:10.1016/j.cpc.2024.109187

    arXiv 2024

  15. [15]

    M. B. Giles, E. Süli, Adjoint methods for PDEs: a review, Acta Numerica 11 (2002) 145–236.doi: 10.1017/S0962492902000033

    work page doi:10.1017/s0962492902000033 2002

  16. [16]

    Modelica Association, Modelica — a unified object-oriented language for systems modeling (language specification v3.6),https://specification.modelica.org(2023)

    2023

  17. [17]

    M. M. Tiller, Introduction to Physical Modeling with Modelica, Springer, 2001. 35

    2001

  18. [18]

    Bradbury, R

    J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne, Q. Zhang, Jax: composable transformations of python+numpy programs,http://github.com/google/jax(2018)

    2018

  19. [19]

    R. T. Q. Chen, Y. Rubanova, J. Bettencourt, D. Duvenaud, Neural ordinary differential equations, in: Advances in Neural Information Processing Systems (NeurIPS), 2018.arXiv:1806.07366

    Pith/arXiv arXiv 2018

  20. [20]

    Kidger, On neural differential equations, Ph.D

    P. Kidger, On neural differential equations, Ph.D. thesis, University of Oxford (2022)

    2022

  21. [21]

    Kidger, diffrax: numerical differential equations in jax,https://docs.kidger.site/diffrax (2021)

    P. Kidger, diffrax: numerical differential equations in jax,https://docs.kidger.site/diffrax (2021)

    2021

  22. [22]

    J. M. Rader, T. Lyons, P. Kidger, Lineax: unified linear solves and linear least-squares in JAX and Equinox (2023).arXiv:2311.17283

    Pith/arXiv arXiv 2023

  23. [23]

    Rader, T

    J. Rader, T. Lyons, P. Kidger, Optimistix: modular optimisation in JAX and Equinox (2024).arXiv: 2402.09983

    Pith/arXiv arXiv 2024

  24. [24]

    Huang, C

    Y. Huang, C. Zou, Y. Li, T. Wik, MINN: Learning the dynamics of differential-algebraic equations and application to battery modeling, IEEE Transactions on Pattern Analysis and Machine Intelligence 46 (12) (2024) 11331–11344.arXiv:2304.14422,doi:10.1109/TPAMI.2024.3456475

    Pith/arXiv arXiv 2024

  25. [25]

    K. Luo, J. Tang, M. Cai, X. Zeng, M. Xie, M. Yan, DAE-KAN: A kolmogorov-arnold network model for high-index differential-algebraic equations (2025).arXiv:2504.15806

    Pith/arXiv arXiv 2025

  26. [26]

    Huang, J

    Y. Huang, J. Ding, A multiscale differential-algebraic neural network-based method for learning dy- namical systems, International Journal of Mechanical System Dynamics 4 (1) (2024) 77–87.doi: 10.1002/msd2.12102

    work page doi:10.1002/msd2.12102 2024

  27. [27]

    J. Koch, M. Shapiro, et al., Learning neural differential algebraic equations via operator splitting, in: IEEE Conference on Decision and Control (CDC), 2024.arXiv:2403.12938

    Pith/arXiv arXiv 2024

  28. [28]

    C. F. Curtiss, J. O. Hirschfelder, Integration of stiff equations, Proceedings of the National Academy of Sciences 38 (3) (1952) 235–243.doi:10.1073/pnas.38.3.235

    work page doi:10.1073/pnas.38.3.235 1952

  29. [29]

    Hairer, G

    E. Hairer, G. Wanner, Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Prob- lems, 2nd Edition, Springer Series in Computational Mathematics, Springer, 1996

    1996

  30. [30]

    J. D. Pryce, A structural analysis method for DAEs, Applied Numerical Mathematics 39 (1) (2001) 41–71.doi:10.1016/S0168-9274(01)00076-9

    work page doi:10.1016/s0168-9274(01)00076-9 2001

  31. [31]

    N. S. Nedialkov, J. D. Pryce, G. Tan, Algorithm 948: DAESA—a MATLAB tool for structural analysis of differential-algebraic equations, ACM Transactions on Mathematical Software 41 (3) (2015) Article 12.doi:10.1145/2700586

    work page doi:10.1145/2700586 2015

  32. [32]

    C. W. Gear, B. Leimkuhler, G. K. Gupta, Automatic integration of the Euler–Lagrange equations with constraints, Journal of Computational and Applied Mathematics 12–13 (1985) 77–90.doi:10.1016/ 0377-0427(85)90008-1

    1985

  33. [33]

    A. C. Hindmarsh, P. N. Brown, K. E. Grant, S. L. Lee, R. Serban, D. E. Shumaker, C. S. Wood- ward, Sundials: Suite of nonlinear and differential/algebraic equation solvers, ACM Transactions on Mathematical Software (TOMS) 31 (3) (2005) 363–396

    2005

  34. [34]

    Petzold, Y

    L. Petzold, Y. Cao, S. Li, R. Serban, S. Li, Sensitivity analysis of differential-algebraic equations and partial differential equations, Computers & Chemical Engineering 30 (10–12) (2006) 1553–1559. doi:10.1016/j.compchemeng.2006.05.015. 36

    work page doi:10.1016/j.compchemeng.2006.05.015 2006

  35. [35]

    J.E.Tolsma, P.I.Barton, Efficientsensitivityanalysisoflarge-scaledifferential-algebraicsystems, Jour- nal of Computational and Applied Mathematics 81 (1) (1997) 97–106.doi:10.1016/S0377-0427(97) 00050-0

    work page doi:10.1016/s0377-0427(97 1997

  36. [36]

    Zhang, H

    H. Zhang, H. Sundar, R. Bui, C. Gao, L. Gao, B. Smith, L. Chandra, PETSc TSAdjoint: a discrete adjoint ODE solver for first-order and second-order sensitivity analysis, ACM Transactions on Mathe- matical Software 48 (2) (2022) 1–29.arXiv:1912.07696,doi:10.1145/3547056

    Pith/arXiv arXiv 2022

  37. [37]

    Griewank, A

    A. Griewank, A. Walther, Algorithm 799: revolve: an implementation of checkpointing for the reverse or adjoint mode of computational differentiation, ACM Transactions on Mathematical Software 26 (1) (2000) 19–45.doi:10.1145/347837.347846

    arXiv 2000

  38. [38]

    J. M. Sanz-Serna, Symplectic Runge–Kutta schemes for adjoint equations, automatic differentiation, optimal control and more, SIAM Review 58 (1) (2016) 3–33.doi:10.1137/140980085

    work page doi:10.1137/140980085 2016

  39. [39]

    Blondel, Q

    M. Blondel, Q. Berthet, M. Cuturi, F. Roy, B. Vincent-Colillas, Efficient and modular implicit differen- tiation, in: Advances in Neural Information Processing Systems (NeurIPS), 2021.arXiv:2105.15183

    Pith/arXiv arXiv 2021

  40. [40]

    C.C.Margossian, M.Betancourt, Efficientautomaticdifferentiationofimplicitfunctions(2021).arXiv: 2112.14217

    Pith/arXiv arXiv 2021

  41. [41]

    diffrax contributors, Feature request: differential algebraic equations, GitHub issue #62, patrick- kidger/diffrax (open since Feb 2022),https://github.com/patrick-kidger/diffrax/issues/62 (2022)

    2022

  42. [42]

    Y. Gamaleldin, Torchdae: numerical differential-algebraic equation solvers in pytorch,https:// github.com/yousef-rafat/torchdae; PyTorch Discusshttps://discuss.pytorch.org/t/225039 (2026)

    2026

  43. [43]

    Y. Ma, S. Gowda, R. Anantharaman, C. Laughman, S. Maddhashiya, V. B. Shah, C. Rackauckas, ModelingToolkit.jl: A composable graph transformation system for equation-based modeling (2021). arXiv:2103.05244

    Pith/arXiv arXiv 2021

  44. [44]

    Rackauckas, et al., Universal differential equations for scientific machine learning (2021).arXiv: 2001.04385

    C. Rackauckas, et al., Universal differential equations for scientific machine learning (2021).arXiv: 2001.04385

    Pith/arXiv arXiv 2021

  45. [45]

    Y. Ma, V. Dixit, M. J. Innes, X. Guo, C. Rackauckas, A comparison of automatic differentiation and continuous adjoint sensitivity analysis (2021).arXiv:1812.01892

    Pith/arXiv arXiv 2021

  46. [46]

    J. A. E. Andersson, J. Gillis, G. Horn, J. B. Rawlings, M. Diehl, Casadi — a software framework for nonlinear optimization and optimal control, Mathematical Programming Computation 11 (2019) 1–36

    2019

  47. [47]

    Sulzer, S

    V. Sulzer, S. G. Marquis, R. Timms, M. Robinson, S. J. Chapman, Python battery mathematical modelling (PyBaMM), Journal of Open Source Software 6 (61) (2021) 3422.doi:10.21105/joss. 03422

    doi:10.21105/joss 2021

  48. [48]

    C. D. Freeman, E. Frey, A. Raichuk, S. Girgin, O. Bachem, Brax: A differentiable physics engine for large scale rigid body simulation, in: Advances in Neural Information Processing Systems (NeurIPS), Datasets and Benchmarks Track, 2021.arXiv:2106.13281

    Pith/arXiv arXiv 2021

  49. [49]

    Toshev, H

    A. Toshev, H. Ramachandran, J. A. Erbesdobler, G. Galletti, J. Brandstetter, N. A. Adams, JAX- SPH: A differentiable smoothed particle hydrodynamics framework, in: ICLR 2024 Workshop on AI for Differential Equations in Science, 2024

    2024

  50. [50]

    D. A. Bezgin, S. M. Schmidt, N. A. Adams, Jax-fluids 2.0: What’s new and what’s next, Computer Physics Communications 308 (2024) 109433.arXiv:2402.05193,doi:10.1016/j.cpc.2024.109433. 37

    Pith/arXiv arXiv 2024

  51. [51]

    S. Ehrlichmann, et al., JuTrack.jl: A Julia package for auto-differentiable accelerator modeling and particle tracking, Computer Physics Communications 305 (2024) 109338.doi:10.1016/j.cpc.2024. 109338

    work page doi:10.1016/j.cpc.2024 2024

  52. [52]

    Gao, et al., chemtrain: Learning deep potential models via automatic differentiation, Computer Physics Communications 310 (2025) 109512.doi:10.1016/j.cpc.2025.109512

    A. Gao, et al., chemtrain: Learning deep potential models via automatic differentiation, Computer Physics Communications 310 (2025) 109512.doi:10.1016/j.cpc.2025.109512

    arXiv 2025

  53. [53]

    Kidger, J

    P. Kidger, J. Morrill, J. Verma, T. Lyons, Neural controlled differential equations for irregular time series, in: Advances in Neural Information Processing Systems (NeurIPS), 2020.arXiv:2005.08926

    Pith/arXiv arXiv 2020

  54. [54]

    Alexe, A

    M. Alexe, A. Sandu, On the discrete adjoints of adaptive time stepping algorithms, Journal of Compu- tational and Applied Mathematics 233 (4) (2009) 1005–1020.doi:10.1016/j.cam.2009.08.109

    work page doi:10.1016/j.cam.2009.08.109 2009

  55. [55]

    C. W. Gear, The automatic integration of ordinary differential equations, Communications of the ACM 14 (3) (1971) 176–179.doi:10.1145/362566.362571

    arXiv 1971

  56. [56]

    Hairer, G

    E. Hairer, G. Wanner, Stiff differential equations solved by Radau methods, Journal of Computational and Applied Mathematics 111 (1–2) (1999) 93–111.doi:10.1016/S0377-0427(99)00134-X

    doi:10.1016/s0377-0427(99)00134-x 1999

  57. [57]

    Lang, Rosenbrock–wanner methods: construction and mission (2020).arXiv:2002.12028

    J. Lang, Rosenbrock–wanner methods: construction and mission (2020).arXiv:2002.12028

    Pith/arXiv arXiv 2020

  58. [58]

    K. Gustafsson, Control-theoretic techniques for stepsize selection in implicit Runge–Kutta methods, ACM Transactions on Mathematical Software 20 (4) (1994) 496–517.doi:10.1145/198429.198437

    arXiv 1994

  59. [59]

    Söderlind, Digital filters in adaptive stepsize control, SIAM Journal on Scientific Computing 25 (3) (2003) 888–904.doi:10.1137/S1064827502407685

    G. Söderlind, Digital filters in adaptive stepsize control, SIAM Journal on Scientific Computing 25 (3) (2003) 888–904.doi:10.1137/S1064827502407685

    work page doi:10.1137/s1064827502407685 2003

  60. [60]

    Galán, F

    S. Galán, F. Feehery, P. I. Barton, Parametric sensitivity methods for the integration of flexible multi- body dynamic systems, Applied Numerical Mathematics 31 (1–2) (1999) 17–47, jump-condition adjoint original paper

    1999

  61. [61]

    Corner, A

    J. Corner, A. Sandu, C. Sandu, Hybrid adjoint sensitivity analysis for differential-algebraic equations with implicit event detection, SIAM Journal on Scientific Computing 40 (3) (2018) B689–B713.doi: 10.1137/17M1140276

    work page doi:10.1137/17m1140276 2018

  62. [62]

    Beigel, M

    D. Beigel, M. S. Mommer, L. Wirsching, H. G. Bock, Approximation of weak adjoints by reverse automatic differentiation of BDF methods (2011).arXiv:1109.3061

    Pith/arXiv arXiv 2011

  63. [63]

    D. Phan, N. Pradhan, M. Jankowiak, Composable effects for flexible and accelerated probabilistic programming in NumPyro, in: NeurIPS 2019 Workshop on Program Transformations, 2019.arXiv: 1912.11554

    Pith/arXiv arXiv 2019

  64. [64]

    Hairer, G

    E. Hairer, G. Wanner, INdAM Bari consortium, Ode/dae test set,https://unige.ch/~hairer/ software/settest/andhttp://archimede.uniba.it/~testset(2023)

    2023

  65. [65]

    Robertson, Robertson autocatalytic reaction test problem, In Hairer-Wanner testset:http: //archimede.uniba.it/~testset/problems/rober.php(1966)

    A. Robertson, Robertson autocatalytic reaction test problem, In Hairer-Wanner testset:http: //archimede.uniba.it/~testset/problems/rober.php(1966)

    1966

  66. [66]

    Breuning, scipy_dae: a python package for solving differential-algebraic equations,https://github

    J. Breuning, scipy_dae: a python package for solving differential-algebraic equations,https://github. com/JonasBreuning/scipy_dae(2024)

    2024

  67. [67]

    National Renewable Energy Laboratory, scikit-sundae: Python bindings for sundials ida,https:// github.com/NREL/scikit-sundae(2024)

    2024

  68. [68]

    rep., Nuclear Energy Agency (2000)

    OECD/NEA, Nea/nsc/doc(2000)21: Main steam line break (mslb) benchmark, Tech. rep., Nuclear Energy Agency (2000). 38

    2000

  69. [69]

    jaxdaeauthors, Threemileislandunit1pwrloopparameters, Internalreactor-physicstestbedassembled from publicly available TMI-1 design data; not a standardized validation benchmark (2026)

    2026

  70. [70]

    IEEE Power & Energy Society, Ieee 14-bus test system, Standard power system benchmark for transient stability; available via MATPOWER and PSAT,https://matpower.org/(1962). 39

    1962