  We analyze the cosmic-ray magnetohydrodynamic (CR MHD) equations to improve
the numerical simulations. We propose to solve them in the fully conservation
form, which is equivalent to the conventional CR MHD equations. In the fully
conservation form, the CR energy equation is replaced with the CR "number"
conservation, where the CR number density is defined as the three fourths power
of the CR energy density. The former contains an extra source term, while
latter does not. An approximate Riemann solver is derived from the CR MHD
equations in the fully conservation form. Based on the analysis, we propose a
numerical scheme of which solutions satisfy the Rankine-Hugoniot relation at
any shock. We demonstrate that it reproduces the Riemann solution derived by
Pfrommer et al. (2006) for a 1D CR hydrodynamic shock tube problem. We compare
the solution with those obtained by solving the CR energy equation. The latter
solutions deviate from the Riemann solution seriously, when the CR pressure
dominates over the gas pressure in the post-shocked gas. The former solutions
converge to the Riemann solution and are of the second order accuracy in space
and time. Our numerical examples include an expansion of high pressure sphere
in an magnetized medium. Fast and slow shocks are sharply resolved in the
example. We also discuss possible extension of the CR MHD equations to evaluate
the average CR energy.
