Large particle accelerators and storage rings are represented by exceedingly complicated Hamiltonians. On the one hand the relativistic motion in electromagnetic fields leads to nonlinear connections between positions and momenta due to the presence of a square root. On the other hand, different from other examples including typical many body problems, an accelerator has thousands of separate sections with different localized fields. Finally it is often necessary to predict orbit stability for around billions of revolutions. All these make even attempts at symplectic tracking very difficult. One method to describe the nonlinear motion is via perturbation theory describing the motion on a Poincare section for one revolution. Recent advances based on differential algebraic methods have allowed the generalization of such perturbation analysis from typically orders two and three to in principle arbitrary order, and in practice to orders around ten. Such representations of the Poincare map allow symplectic tracking in various ways. More importantly, the high order Poincare maps can be used to perform rigorous stability analysis. For this purpose, the differential algebraic approach allows performing automated high order normal form analysis of the motion, which upon adjusting system parameters of the accelerator leads to dynamics that is nearly integrable. Utilizing modern methods of rigorous global optimization based on high-order Taylor model methods, it is possible to perform rigorous Nekhoroshev-type stability estimates. Results are shown about such stability estimates, which can establish guaranteed stability within the regions typically occupied by the particle beam for operationally required times.