  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 numerical integration of a system of ordinary
differential equations. In principle, this procedure yields first order
methods, but the analysis paves the way to 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 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. Besides, 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.
