This work is concerned with the proof of optimal error bounds for the discretization of $H^1$-elliptic minimization problem with solutions taking values in a Riemannian manifold. The discretization is done using Geodesic Finite Elements, a method of arbitrary order that is invariant under isometries. The discretization error is considered both intrinsically in a specially introduced Sobolev-distance as well as extrinsically. Optimal estimates of $H^1$- and $L^2$-type are shown, that have been observed experimentally in previous works of other authors. Using the Rothe method consisting of an implicit Euler method for the time discretization and Geodesic Finite Elements for the spatial discretization, error estimates for $L^2$-gradient flows of $H^1$-elliptic energies are derived as well. The core of the work is formed by the discretization error estimates for minimization problems in instrinsic $H^1$- and $L^2$-distances. To derive these, inverse estimates and interpolation errors for Geodesic Finite Elements and their discrete variations are shown. Using a nonlinear Cea's Lemma, this leads to the $H^1$-error estimate for minimizers of $H^1$-elliptic energies. A generalization of the Aubin-Nitsche-Lemma shows optimal $L^2$-error estimates for (essentially) semilinear energies, as long as the dimension of the domain of the minimizer is limited to $d<4$ for technical reasons. All results are illustrated using harmonic maps into a smooth Riemannian manifold satisfying certain curvature bounds as an example.