We compare three methods for computing invariant Lyapunov exponents (LEs) in general relativity. These methods involve the geodesic deviation vector technique (M1), the two-nearby-orbits method with projection operations and with coordinate time as the independant variable (M2), and the two-nearby-orbits method without projection operations and with proper time as the independent variable (M3). An analysis indicates that M1 and M3 do not need any projection operation. In general, the values of LEs from the three methods are almost the same. However, M2 fails for some specific cases. As a result, M3 is the most preferable to calculate LEs in most cases. In addition, we propose to construct the invariant fast Lyapunov indictor (FLI) with two-nearby-trajectories and give its algorithm in order to quickly distinguish chaos from order. Taking a static axisymmetric spacetime as a physical model, we apply different algorithms of the FLI to explore the global dynamics of phase space of the system where regions of chaos and order are clearly identified.