Geometric integration of ODEs using multiple quadratic auxiliary variables

We present a novel numerical method for solving ODEs while preserving polynomial first integrals. The method is based on introducing multiple quadratic auxiliary variables to reformulate the ODE as an equivalent but higher-dimensional ODE with only quadratic integrals to which the midpoint rule is applied. The quadratic auxiliary variables can subsequently be eliminated yielding a midpoint-like method on the original phase space. The resulting method is shown to be a novel discrete gradient method. Furthermore, the averaged vector field method can be obtained as a special case of the proposed method. The method can be extended to higher-order through composition and is illustrated through a number of numerical examples.