Some years ago, the harmonic polynomial was introduced to study the harmonic topological index. Here, using this polynomial, we obtain several properties of the harmonic index of many classical symmetric operations of graphs: Cartesian product, corona product, join, Cartesian sum and lexicographic product. Some upper and lower bounds for the harmonic indices of these operations of graphs, in terms of related indices, are derived from known bounds on the integral of a product on nonnegative convex functions. Besides, we provide an algorithm that computes the harmonic polynomial with complexity O(n 2 ).This work was supported in part by two grants from Ministerio de Economía y Competititvidad, Agencia Estatal de Investigación (AEI) and Fondo Europeo de Desarrollo Regional (FEDER) (MTM2016-78227-C2-1-P and MTM2017-90584-REDT), Spain