Graphene plasmons are important collective excitations in graphene, which play a key role in determining the optical properties of graphene. They have quite lots of unique features in comparison with classical plasmons in noble metals. Of them, the active tunability is the most attractive, which is realized by external gating (equivalently electric field). As is well known, graphene also has strong magnetic response (e.g. room temperature quantum Hall effect), so magnetic field can act as another degree of freedom for actively tuning graphene plasmons, with the new quasi particles being so-called graphene magneto-plasmons. Because of the two-dimensional nature of graphene, the numerical studies (or full wave simulations) of graphene magneto-plasmons are usually carried out through a three-dimensional approximation, e.g. treating two-dimensional graphene as a very thin three-dimensional film. Actually, this treatment takes quite some time and requires high memory consumption. Herein, starting from Coulomb law and charge conservation law, we propose an alternative numerical method, namely, two-dimensional finite element method, to solve this problem. All the calculations are now performed in two-dimensional graphene plane, and the usual three-dimensional approximation is not required. To characterize the excitations of graphene magneto-plasmons, the eigenvalue loss spectrum is introduced. Based on this method, graphene magneto-plasmons in graphene rings of four kinds are investigated. The strongest magneto-optic effect is observed in circular ring, which is consistent with its highest rotational symmetry. In all the rings, the lowest dipolar graphene magneto-plasmon always supports symmetric mode splitting, which can be further modified by the interaction between inner edge and outer edge of ring. As the hole size is very small, the edge current confined to the outer edge dominates, and that confined to the inner edge can be ignored; while increasing the hole size, the interaction between these two edges increases, which results in the reduction of the symmetric mode splitting; when the hole size is larger than a critical value, the symmetric mode splitting will disappear.