Collisionless magnetohydrodynamic shocks are ubiquitous in the universe, occurring in locations as nearby as the Earths bow shock and as remote as the supermassive blackholes at the centers of distant galaxies. All of these shocks are believed to accelerate particles via diffusive (or first-order Fermi) shock acceleration (DSA). This thesis makes contributions to the study of collisionless magnetohydrodynamic shocks including the development of numerical solutions to the Rankine-Hugoniot shock jump conditions in relativistic shocks, independent development of a kinetic Monte Carlo simulation to model acceleration at heliospheric and astrophysical shocks, and applications of this simulation to problems in space physics. Inside the heliosphere, the simulation is used to generate simulated spectra for heliospheric shocks and these are compared to in-situ observations from the SWICS and HI-SCALE instruments on the Ulysses spacecraft. These results show that the extended power-law tail in shock 91292 can be accurately modeled by DSA using a simple scattering mechanism and only one free parameter, without the need for a pre-acceleration mechanism. To improve the fit in the energy range of slightly suprathermal particles, we explore an alternative scattering mechanism that emulates the diffusive characteristics of the field line wandering (FLW) study of Giacalone and Jokipii (1999). In supernova remnants, the effects of the cross-shock potential, an electric field created by the disparate diffusive scales of protons and electrons which create charge separation near the shock, are explored laying the groundwork for future inclusion of a self-consistent cross-shock potential in the simulation. Finally, in relativistic shocks, simulation results highlight the importance of the obliquity of the large scale magnetic field, the amount of turbulence, and even the detailed microphysics of the turbulence in determining the power-law index of accelerated particles. We also find that the super-luminality vs. sub-luminality of the shock becomes a critical characteristic in determining the behavior of the power-law index. For sub-luminal shocks, decreased turbulence and increased obliquity enhance acceleration. For super-luminal shocks, the dependence is reversed. This sudden change in behavior is understood and explained.