Accurate modeling of vibronically driven magnetic relaxation from ab initio calculations is of paramount importance to the design of next-generation single-molecule magnets (SMMs). Previous theoretical studies have been relying on numerical differentiation to obtain spin-phonon couplings in the form of derivatives of spin Hamiltonian parameters. In this work, we introduce a novel approach to obtain these derivatives fully analytically by combining the linear vibronic coupling (LVC) approach with analytic complete active space self-consistent field derivatives and nonadiabatic couplings computed at the equilibrium geometry with a single electronic structure calculation. We apply our analytic approach to the computation of Orbach and Raman relaxation rates for a bis-cyclobutadienyl Dy(III) sandwich complex in the frozen-solution phase, where the solution environment is modeled by electrostatic multipole expansions, and benchmark our findings against results obtained using conventional numerical derivatives and a fully electronic description of the whole system. We demonstrate that our LVC approach exhibits high accuracy over a wide range of coupling strengths and enables significant computational savings due to its analytic, "single-shot" nature. Evidently, this offers great potential for advancing the simulation of a wide range of vibronic coupling phenomena in magnetism and spectroscopy, ultimately aiding the design of high-performance SMMs. Considering different environmental representations, we find that a point charge model shows the best agreement with the reference calculation, including the full electronic environment, but note that the main source of discrepancies observed in the magnetic relaxation rates originates from the approximate equilibrium electronic structure computed using the electrostatic environment models rather than from the couplings.