Existing element types available in finite element codes typically utilize polynomial shape functions to define the displacement field in the problem of interest. The polynomial shape functions serve the purpose adequately in static analysis where the displacements and the stresses in a structure are of primary interest. These shape functions give rise to increasing inaccuracy for higher natural frequencies. It is shown that harmonic shape functions yield better results for the higher natural frequencies with same element count. Axial vibrations of bars and transverse vibrations of beams have been analyzed using harmonic shape functions. A comparative analysis has been made between results predicted by the harmonic interpolation functions and polynomial interpolation functions using same number of nodes.