A numerical method that solves the KZK equation and satisfies the Rankine-Hugoniot relation for shock-wave propagation is described and characterized. By comparison with a known planar solution, it is demonstrated that current numerical methods in both the time and frequency domains predict a shock front that is stationary relative to the propagation phase and that the proposed method predicts the correct speed. These methods are then compared in the context of shock-wave lithotripsy and high intensity focused ultrasound. At the focus, axisymmetric shock-wave lithotripter simulations show that the Rankine-Hugoniot method predicts peak positive pressures that are 20% smaller, peak negative pressures that are 10% larger, and a full width at half maximum that is 45% larger. High-intensity focused ultrasound simulations for an axisymmetric transducer in water have even larger variations in peak positive pressures and intensity but smaller variations in peak negative pressure. Simulations for a rectangular transducer in tissue, where absorption is more prominent, exhibit smaller but significant variations in peak pressures and intensity. Satisfying the Rankine-Hugoniot condition plays an important role in numerical simulations of the KZK equation. [This work was supported by NIH Grants R01-HL075485 and 1R01-CA114093-02.]