We apply a combination of N-body modeling techniques and automated data fitting with Monte Carlo Markov Chain uncertainty analysis of Keplerian orbital models to radial velocity data to determine long term stability of the planetary system GJ 581. We find that while there are stability concerns with the 4-planet model as published by Forveille et al. (2011), when uncertainties in the system are accounted for, particularly stellar jitter, the hypothesis that the 4-planet model is gravitationally unstable is not statistically significant. Additionally, the system including proposed planet g by Vogt et al. (2012) also shows some stability concerns when eccentricities are allowed to float in the orbital fit, yet when uncertainties are included in the analysis the system including planet g also can not be proven to be unstable. We present revised reduced chi-squared values for Keplerian astrocentric orbital fits assuming 4-planet and 5-planet models for GJ~581 under the condition that best fits must be stable, and find no distinguishable difference by including planet g in the model. Additionally we present revised orbital element estimates for each assuming uncertainties due to stellar jitter under the constraint of the system being gravitationally stable.