The notion of error in practical molecular and Langevin dynamics simulations of large biomolecules is far from understood because of the relatively large value of the timestep used, the short simulation length, and the low-order methods employed. We begin to examine this issue with respect to equilibrium and dynamic time-correlation functions by analyzing the behavior of selected implicit and explicit finite-difference algorithms for the Langevin equation. We derive : local stability criteria for these integrators; analytical expressions for the averages df the potential, kinetic, and total energy; and various limiting cases (e.g., timestep and damping constant approaching zero), for a system of coupled harmonic oscillators. These results are then compared to the corresponding exact solutions for the continuous problem, and their implications to molecular dynamics simulations are discussed. New concepts of practical and theoretical importance are introduced : scheme-dependent perturbative damping and perturbative frequency functions. Interesting differences in the asymptotic behavior among the algorithms become apparent through this analysis, and two symplectic algorithms,"LIM2" (implicit) and "BBK" (explicit), appear most promising on theoretical grounds. One result of theoretical interest is that for the Langevin/implicit-Euler algorithm ("LI") there exist timesteps for which there is neither numerical damping nor shift in frequency for a harmonic oscillator. However, this idea is not practical for more complex systems because these special timesteps can account only for one frequency of the system,and a large damping constant is required. We therefore devise a more practical, delay-function approach to remove the artificial damping and frequency perturbation from LI. Indeed, a simple MD implementation for a system of coupled harmonic oscillators demonstrates very satisfactory results in comparison with the velocity-Verlet scheme. We also define a probability measure to estimate individual trajectory error. This framework might be useful in practice for estimating rare events, such as barrier crossing. To illustrate, this concept is applied to a transition-rate calculation, and transmission coefficients for the five schemes are derived.