The article considers the construction of a parallel algorithm on the GPU computing architecture forfinding the Wigner function of a quantum system with a polynomial potential. A numerical-analyticalmethod for constructing the Wigner function, which is based on calculating the trace of the product ofthe density matrix and the matrix of the Weyl operator, is described. The operators are represented inthe basis of a quantum harmonic oscillator, for which the Moyal equation transforms into the Liouvilleequation. This approach enables to visually analyze the degree of anharmonicity of the system in termsof the off-diagonal elements of the density matrix. The parallel implementation on the GPU massivelyparallel architecture has reduced the calculation time by two orders of magnitude compared to thesingle-threaded version on the x86 architecture.