  Magnetic skyrmions widely exist in a diverse range of magnetic systems,
including chiral magnets with a non-centrosymmetric structure characterized by
Dzyaloshinkii-Moriya interaction~(DMI). In this study, we propose a generalized
semi-implicit backward differentiation formula projection method, enabling the
simulations of the Landau-Lifshitz~(LL) equation in chiral magnets in a typical
time step-size of $1$ ps, markedly exceeding the limit subjected by existing
numerical methods of typically $0.1$ ps. Using micromagnetics simulations, we
show that the LL equation with DMI reveals an intriguing dynamic instability in
magnetization configurations as the damping varies. Both the isolated
skyrmionium and skyrmionium clusters can be consequently produced using a
simple initialization strategy and a specific damping parameter. Assisted by
the string method, the transition path between skyrmion and skyrmionium, along
with the escape of a skyrmion from the skyrmion clusters, are then thoroughly
examined. The numerical methods developed in this work not only provide a
reliable paradigm to investigate the skyrmion-based textures and their
transition paths, but also facilitate the understandings for magnetization
dynamics in complex magnetic systems.
