In the formulation of synchrotron radiation (SR), the Wigner function represents the photon flux density in the 4-dimensional phase space, and thus defines the brilliance and coherence. Although we have to rely on a numerical method to evaluate the Wigner function for a realistic electron beam, a convolution approach is available under the assumption that each electron travels along a straight line apart from a sinusoidal motion driven by an undulator magnetic field. As a result, the computation of Wigner functions for typical SR sources is usually completed within an acceptable time. It should be noted, however, that the above assumption is violated in a multiple-beam-waist (MBW) configuration, in which the lattice function is designed to have more than one beam waist position where the vertical betatron function is optimized to accommodate a narrow-gap undulator. In such a case, we have to carry out a Monte Carlo simulation, which requires a huge amount of time to get a reliable result. To overcome this difficulty, we introduce a numerical method based on a theoretical approach, to aim at reducing the numerical cost to compute the Wigner functions in the MBW configuration. We show that the computation time can be saved by 2-4 orders of magnitude by means of the proposed method.