We construct splitting methods for the solution of the equations of magnetohydrodynamics (MHD). Due to the physical significance of the involved operators, splittings into three or even four operators with positive coefficients are appropriate for a physically correct and efficient solution of the equations. To efficiently obtain an accurate solution approximation, adaptive choice of the time-steps is important particularly in the light of the unsmooth dynamics of the system. Thus, we construct new method coefficients in conjunction with associated error estimators by optimizing the leading local error term. As a proof of concept, we demonstrate that adaptive splitting faithfully reflects the solution behavior also in the presence of a shock for the viscous Burgers equation, which serves as a simplified model problem displaying several features of the Navier-Stokes equation for incompressible flow.