A new method for description of viscoelastic functions with Prony series of exponentials is presented and compared with existing methods. This method places constraints on the coefficients of the series which guarantees smoothness of the functions and a discrete spectral representation. An optimization technique is used to obtain the solution. The resulting series permit the use of Whittaker's method of solution of integral equations which is specialized to the case of viscoelastic analysis, and an exact relationship is obtained between creep and relaxation functions.