This paper addresses the approximation of fractional harmonic maps. Besides a unit-length constraint, one has to tackle the difficulty of nonlocality. We establish weak compactness results for critical points of the fractional Dirichlet energy on unit-length vector fields. We devise and analyze numerical methods for the approximation of various partial differential equations related to fractional harmonic maps. The compactness results imply the convergence of numerical approximations. Numerical examples on spin chain dynamics and point defects are presented to demonstrate the effectiveness of the proposed methods.