We simulate spin relaxation processes, which may be measured by either continuous wave or pulsed magnetic resonance techniques, using trajectory-based simulation methodologies. The spin-lattice relaxation rates are extracted numerically from the relaxation simulations. The rates obtained from the numerical fitting of the relaxation curves are compared to those obtained by direct simulation from the relaxation Bloch-Wangsness-Abragam-Redfield theory (BWART). We have restricted our study to anisotropic rigid-body rotational processes, and to the chemical shift anisotropy (CSA) and a single spin-spin dipolar (END) coupling mechanisms. Examples using electron paramagnetic resonance (EPR) nitroxide and nuclear magnetic resonance (NMR) deuterium quadrupolar systems are provided. The objective is to compare those rates obtained by numerical simulations with the rates obtained by BWART. There is excellent agreement between the simulated and BWART rates for a Hamiltonian describing a single spin (an electron) interacting with the bath through the chemical shift anisotropy (CSA) mechanism undergoing anisotropic rotational diffusion. In contrast, when the Hamiltonian contains both the chemical shift anisotropy (CSA) and the spin-spin dipolar (END) mechanisms, the decay rate of a single exponential fit of the simulated spin-lattice relaxation rate is up to a factor of 0.2 smaller than that predicted by BWART. When the relaxation curves are fit to a double exponential, the slow and fast rates extracted from the decay curves bound the BWART prediction. An extended BWART theory, in the literature, includes the need for multiple relaxation rates and indicates that the multiexponential decay is due to the combined effects of direct and cross-relaxation mechanisms.