The paper studies numerical methods that preserve a Lyapunov function of a dynamical system, i.e., numerical approximations whose energy decreases, just like in the original differential equation. With this aim, a discrete gradient method is implemented for the numerical integration of a system of ordinary differential equations. In principle, this procedure yields first-order methods, but the analysis paves the way for the design of higher-order methods. As a case in point, the proposed method is applied to the Duffing equation without external forcing, considering that, in this case, preserving the Lyapunov function is more important than the accuracy of particular trajectories. Results are validated by means of numerical experiments, where the discrete gradient method is compared to standard Runge–Kutta methods. As predicted by the theory, discrete gradient methods preserve the Lyapunov function, whereas conventional methods fail to do so, since either periodic solutions appear or the energy does not decrease. Moreover, the discrete gradient method outperforms conventional schemes when these do preserve the Lyapunov function, in terms of computational cost; thus, the proposed method is promising.
