The 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 CASSCF derivatives and nonadiabatic couplings computed at a single equilibrium geometry and electronic structure. We apply our analytic approach to the computation of Orbach and Raman relaxation rates in a solvated bis-cyclobutadienyl Dy(III) sandwich complex, where the solution environment is modelled by electrostatic multipole expansions, and benchmark our findings against results obtained using conventional numerical derivatives and a fully electronic description of the whole system. Among the ghost, charge and charge+dipole representations, we find that the point charge model shows the best agreement with the reference calculation, 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. We demonstrate that our novel LVC approach exhibits high accuracy over a wide range of coupling strengths and enables 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.