We introduce a method for mitigating the numerical integration errors of initial value problems. We propose a methodology for constructing an optimal state-space representation that gives minimum numerical truncation error, and in this sense, is the optimal state-space representation for modeling given phase-space dynamics. To that end, we utilize a simple transformation of the state-space equations into their variational form. This process introduces an inherent freedom, similar to the gauge freedom in electromagnetism. We then utilize the gauge function to reduce the numerical integration error. We show that by choosing an appropriate gauge function the numerical integration error dramatically decreases and one can achieve much better accuracy compared to the standard state variables for a given time-step. We illustrate the method using a few examples taken from the space systems and aeroelasticity fields. In all of our illustrating examples, the gauge-optimized integration outperforms the conventional integration.