The equations of motion may be directly integrated by a trapezoidal rule method to derive the phase plane portrait or the time series. There is an error in the size of the damping coefficient brought about by the fact that the real trajectory is everywhere convex except possibly at x=0, and therefore the calculation with tangential segments will always lie outside the real trajectory.
Other methods of numerical simulation are of course possible; these may in principle be less inaccurate, but they will not be completely precise and that just postpones the point at which the real trajectory departs from the simulation in a chaotic region; it does not eliminate the divergence.
Additionally, the real oscillator has added noise. The effect of a non-zero timestep in the calculation is very similar to the effect of adding noise in the real system; so in terms of qualitative accuracy the trapezoidal method has merit.