Abstract
We propose and implement a third-order accurate numerical scheme for the Landau-Lifshitz-Gilbert equation, which describes magnetization dynamics in ferromagnetic materials under large damping parameters. This method offers two key advantages: (1) It solves only constant-coefficient linear systems, enabling fast solvers and thus achieving much higher numerical efficiency than existing second-order methods. (2) It attains third-order temporal accuracy and fourth-order spatial accuracy, and is unconditionally stable for large damping parameters. Numerical examples in 1D and 3D simulations verify both its third-order accuracy and efficiency gains. However, when large damping parameters and pre-projection solutions are involved, both this proposed method and a second-order method of the same style fail to capture reasonable physical structures, despite extensive theoretical analyses. Additionally, comparisons of domain wall dynamics among BDF2, BDF3, and BDF1 show that BDF2 and BDF3 yield failed simulations, while BDF1 performs marginally better.