Automatic differentiation for performing the Cauchy-Kovalevskaya procedure in Lax-Wendroff type discretizations

Lax-Wendroff methods combined with discontinuous Galerkin/flux reconstruction spatial discretization provide a high-order, single-stage, quadrature-free method for solving hyperbolic conservation laws. In this work, we introduce automatic differentiation (AD) for performing the Cauchy-Kowalewski procedure used in the element-local time average flux computation step (the predictor step) of Lax-Wendroff methods. The application of AD is similar for methods of any order and does not need positivity corrections during the predictor step. This contrasts with the approximate Lax-Wendroff procedure, which requires different finite difference formulas for different orders of the method and positivity corrections in the predictor step for fluxes that can only be computed on admissible states. The method is Jacobian-free and problem-independent, allowing direct application to any physical flux function. Numerical experiments demonstrate the order and positivity preservation of the method. Additionally, performance comparisons indicate that the wall-clock time of automatic differentiation is always on par with the approximate Lax-Wendroff method.

💡 Research Summary

**
The paper presents a novel application of automatic differentiation (AD) to the Cauchy‑Kowalevski (C‑K) procedure that underlies the predictor step of Lax‑Wendroff (LW) type schemes when combined with discontinuous Galerkin (DG) or flux‑reconstruction (FR) spatial discretizations. Traditional LW methods achieve high‑order accuracy by expanding the time‑averaged flux
(F = \sum_{m=0}^{N} \frac{\Delta t^{m}}{(m+1)!},\partial_{t}^{m} f(u))
and then using the conservation law to replace temporal derivatives of the flux with spatial derivatives of the solution. This replacement is usually performed analytically via the chain rule, requiring explicit Jacobians of the physical flux function and, for higher orders, increasingly complex Faà di Bruno‑type combinatorial formulas. An alternative “approximate LW” approach replaces these analytic derivatives with finite‑difference (FD) stencils, but the stencil must be changed for each order and positivity corrections are needed when the flux cannot be evaluated on non‑admissible states (e.g., relativistic hydrodynamics).

The authors propose to compute the temporal derivatives of the flux directly with AD, thereby avoiding explicit Jacobian construction and order‑dependent FD formulas. AD evaluates directional derivatives by propagating derivative information through the program that computes the flux, effectively providing Jacobian‑free, exact derivatives of any order. The key technical device is the definition of a “derivative bundle” (\mathcal{B}^{(m)}(u) = (u, u^{(1)}, \dots, u^{(m)})). The m‑th flux derivative can be expressed as a function of this bundle, (f^{(m)} = f^{(m)}(\mathcal{B}^{(m)}(u))). By applying AD recursively to the flux function, the required derivatives are obtained without manually coding the combinatorial expressions of Faà di Bruno. This recursion mirrors the nested approach used in prior works but is fully automated.

Because AD treats the flux as a black‑box, the method is problem‑independent: any physically admissible flux (Euler, MHD, relativistic hydrodynamics, etc.) can be used without deriving analytical Jacobians. Moreover, the implementation is identical for any desired order N, simplifying code maintenance and extending to arbitrarily high orders with minimal effort.

A further advantage concerns positivity preservation. The approximate LW method often requires additional corrections in the predictor step because the FD approximations may evaluate the flux at non‑admissible states (negative density or pressure). In contrast, the AD‑based predictor evaluates the flux only on the approximate solution, which is already limited to the admissible set by a subcell‑based limiter. Consequently, no extra positivity correction is needed during the predictor step. The authors combine the AD‑based predictor with the subcell‑based admissibility‑preserving limiter of