Second-Order Discrete Gradient Method

• Second-Order Discrete Gradient Method is a structure-preserving time discretization technique that exactly preserves energy while providing second order accuracy.
• It uses discrete gradients, such as the symmetrized Itoh–Abe and AVF methods, to convert continuous conservation laws into numerical schemes with O(h^3) local error.
• The method’s versatility is evident in its extensions to PDEs, constrained flows, and stochastic models, ensuring robust stability and invariant preservation across applications.

A second-order discrete gradient method is a structure-preserving time discretization whose local truncation error is $O(h^3)$ and whose update is built from a discrete gradient so that a prescribed first integral, typically a Hamiltonian or an energy functional, is preserved exactly at the discrete level. In the ODE setting emphasized in recent work, the prototype problem is a skew-gradient system $\dot x=S\nabla H(x)$ with $S^T=-S$, and the defining mechanism is the replacement of $\nabla H$ by a map $\overline{\nabla}H(x,y)$ satisfying a discrete chain rule. In the derivative-free setting, the symmetrized Itoh–Abe discrete gradient yields a canonical second-order example, while the order-theoretic analysis of average-vector-field-based schemes identifies second order as the maximal generic order attainable without introducing $h$-dependent corrections in the skew-symmetric factor (Myhr et al., 12 Jan 2026, Eidnes, 2020).

1. Discrete-gradient formulation and exact preservation

For an autonomous ODE

$\dot x=f(x), \qquad x\in\mathbb{R}^n,$

with first integral $H$, a central case is the skew-gradient form

$\dot x=S\nabla H(x),\qquad S^T=-S.$

A discrete gradient $$\overline{\nabla}H:\mathbb{R}^n\times\mathbb{R}^n\to\mathbb{R}^n$$ is defined by the two identities

$\dot x=S\nabla H(x)$0

The corresponding discrete gradient method is

$\dot x=S\nabla H(x)$1

Its conservation law follows directly: $\dot x=S\nabla H(x)$2 Thus any discrete gradient method preserves $\dot x=S\nabla H(x)$3 exactly, modulo solver tolerance (Myhr et al., 12 Jan 2026).

This exact preservation is independent of the formal order of accuracy. The second-order qualifier therefore refers not to invariant preservation, which is already built into the construction, but to the truncation error: a second-order discrete gradient method has local truncation error $\dot x=S\nabla H(x)$4 and global order $\dot x=S\nabla H(x)$5 (Myhr et al., 12 Jan 2026).

2. Second order as a symmetry condition

For the basic discrete gradient method with constant $\dot x=S\nabla H(x)$6, the order is determined by the approximation properties of the discrete gradient along the numerical trajectory. The cited derivative-free analysis states that first- and second-order discrete gradients “will thus give first- and second-order integrators if used in the scheme $\dot x=S\nabla H(x)$7.” In particular, first-order Itoh–Abe gives a first-order integrator, while symmetric second-order discrete gradients such as midpoint, AVF, and symmetrized Itoh–Abe yield second-order integrators (Myhr et al., 12 Jan 2026).

In the order-theoretic framework for skew-gradient systems, the general discrete gradient method

$\dot x=S\nabla H(x)$8

preserves $\dot x=S\nabla H(x)$9 exactly whenever $S^T=-S$0 is skew-symmetric. The same framework shows that if $S^T=-S$1 is independent of $S^T=-S$2, then one cannot, in general, exceed second order; higher order requires explicit $S^T=-S$3-dependent corrections in $S^T=-S$4 (Eidnes, 2020). In that sense, second order is the natural upper limit of the basic two-point construction.

A one-dimensional formulation makes this especially transparent. For

$S^T=-S$5

the standard discrete gradient scheme is

$S^T=-S$6

This method is explicitly identified with the modified midpoint rule; it is denoted GR, GR-1, or GR-2, and “GR-1 is of 2nd order” (Cieśliński et al., 2010).

3. Canonical second-order constructions

A standard second-order example for canonical Hamiltonian systems is the average vector field method

$S^T=-S$7

obtained by using the AVF discrete gradient

$S^T=-S$8

For constant $S^T=-S$9, this is the classical AVF integrator; it is a B-series method, preserves $\nabla H$0 exactly, and has order $\nabla H$1 (Eidnes, 2020).

A derivative-free second-order construction is the symmetrized Itoh–Abe discrete gradient. If $\nabla H$2 denotes the Itoh–Abe discrete gradient, then the symmetrized form is

$\nabla H$3

The associated basic second-order scheme is

$\nabla H$4

Its key reported properties are: exact preservation of $\nabla H$5, second-order accuracy, derivative-free dependence on $\nabla H$6 alone, and implicitness of the update (Myhr et al., 12 Jan 2026).

In one spatial degree of freedom, the modified midpoint rule quoted above provides the canonical second-order discrete gradient method. The same paper states that it preserves the energy exactly and, in the one-dimensional Hamiltonian setting, preserves the trajectories in phase space in the sense that the discrete solution stays on the exact energy level curve (Cieśliński et al., 2010).

4. Order theory and the transition to higher order

The modern order theory of discrete gradient methods is formulated in terms of B-series for constant $\nabla H$7 and P-series for general skew-gradient systems. For canonical Hamiltonian systems, the AVF-based discrete gradient method admits a B-series expansion, and order $\nabla H$8 is characterized by matching the B-series coefficients with those of the exact solution for rooted trees of order up to $\nabla H$9. For non-constant $\overline{\nabla}H(x,y)$0, the corresponding expansion is a P-series built from bi-colored rooted trees that encode derivatives of both $\overline{\nabla}H(x,y)$1 and $\overline{\nabla}H(x,y)$2 (Eidnes, 2020).

Within that framework, the second-order AVF discrete gradient method for general skew-gradient systems can be written as

$\overline{\nabla}H(x,y)$3

and the paper shows that this scheme is second order (Eidnes, 2020).

The derivative-free fourth-order construction developed later makes the dependence of higher order on a second-order base method explicit. Its starting point is a second-order symmetric discrete gradient, chosen there as $\overline{\nabla}H(x,y)$4. One first constructs a third-order method with a specially designed skew matrix $\overline{\nabla}H(x,y)$5, then symmetrizes it to obtain a fourth-order method with

$\overline{\nabla}H(x,y)$6

If one sets $\overline{\nabla}H(x,y)$7, thereby removing the higher-order corrections involving $\overline{\nabla}H(x,y)$8 and $\overline{\nabla}H(x,y)$9, the construction reverts to the second-order SIA-based discrete gradient method (Myhr et al., 12 Jan 2026).

A different higher-order generalization proceeds through discontinuous Galerkin time stepping. In that setting, the classical discrete gradient method appears as the $h$0 member of a DG hierarchy, while the $h$1 scheme is globally second order and third order at the time nodes. The DG formulation preserves or dissipates the energy exactly through an integral-in-time identity and extends the classical discrete gradient paradigm to arbitrary polynomial degree $h$2 (Kemmochi, 2023).

5. Implicit realization, derivative-free Jacobians, and computational behavior

Second-order discrete gradient methods are typically implicit. For the SIA scheme, one solves

$h$3

The exact Jacobian is

$h$4

but the derivative-free implementation replaces $h$5 by symmetric finite differences. The discrete gradient itself remains exact within machine precision because it uses only function values of $h$6; the approximation enters only through the Newton Jacobian (Myhr et al., 12 Jan 2026).

The convergence analysis is formulated in the language of inexact Newton methods. For the derivative-free implementation, the paper states that if $h$7, where $h$8 is the function evaluation precision, then the finite-difference error in derivatives is negligible compared to the basic discretization error $h$9 for reasonable step sizes. The reported conclusion is that overall second-order accuracy is preserved despite the derivative-free Jacobian, and numerical experiments verify the order and show that the derivative-free method is significantly faster than obtaining derivatives by automatic differentiation (Myhr et al., 12 Jan 2026).

In the one-dimensional Hamiltonian setting, the implicit nonlinear system is solved by fixed-point or Newton iterations, iterated until the accuracy $\dot x=f(x), \qquad x\in\mathbb{R}^n,$0 is obtained. There the exact energy conservation property is retained, and the numerical experiments emphasize very good long-time behavior, including stable periods of oscillation and excellent behavior near the separatrix (Cieśliński et al., 2010).

6. PDE, constrained-flow, fractional, and stochastic extensions

Second-order discrete gradient ideas extend well beyond finite-dimensional Hamiltonian ODEs. For two-phase flow in porous media, a refactorized Cauchy one-leg method with $\dot x=f(x), \qquad x\in\mathbb{R}^n,$1 is equivalent to the implicit midpoint method, and a specially constructed discrete chemical potential $\dot x=f(x), \qquad x\in\mathbb{R}^n,$2 satisfies

$\dot x=f(x), \qquad x\in\mathbb{R}^n,$3

The resulting midpoint scheme is second order in time and obeys a discrete Gibbs free energy balance without numerical dissipation (Jones et al., 2023).

For Oseen–Frank gradient flows, the rotational discrete gradient scheme rewrites the unit-length constrained evolution in the form

$\dot x=f(x), \qquad x\in\mathbb{R}^n,$4

Its second-order time discretization is

$\dot x=f(x), \qquad x\in\mathbb{R}^n,$5

and the paper proves that the scheme is strictly length-preserving and unconditionally energy-stable. Alongside mean-value and Gonzalez discrete gradients, it proposes a problem-specific Oseen–Frank discrete gradient satisfying the exact energy difference relation (Xu et al., 2023).

For nonlinear integro-differential models with fractional memory, second-order variable-step approximations of the Riemann–Liouville integral and Caputo derivative are combined with discrete complementary kernels to build a discrete gradient structure for the convolution operators. This yields discrete variational energy dissipation laws for time-fractional Allen–Cahn and Klein–Gordon type models, and the schemes are shown to be asymptotically compatible with the corresponding classical energy laws in the relevant fractional-order limits (Liao et al., 2023).

The terminology also appears in discrete sampling. In Preconditioned Discrete-HAMS, the algorithm is described as a second-order discrete gradient method because it uses a global quadratic approximation of the log target, a preconditioning matrix $\dot x=f(x), \qquad x\in\mathbb{R}^n,$6, and a momentum correction term

$\dot x=f(x), \qquad x\in\mathbb{R}^n,$7

that vanishes when the target is exactly quadratic. In that setting, the method is irreversible, satisfies generalized detailed balance, and is rejection-free for quadratic potentials with curvature $\dot x=f(x), \qquad x\in\mathbb{R}^n,$8 (Zhou et al., 29 Jul 2025).

Second-order discrete gradient methods therefore occupy a central position in structure-preserving computation: they are the basic exact-invariant or exact-energy-balance schemes from which higher-order discrete gradient integrators, DG-in-time variants, constrained rotational flows, and derivative-free constructions are built. The cited literature consistently treats second order not as an endpoint of accuracy, but as the foundational level at which exact structure preservation and practical numerical efficiency first coexist.