Since the torque characteristic of a surface permanent magnet synchronous motor (SPMSM) is strongly dependent on its magnetization distribution (MD), the prior estimation of MD is essential for its practical design. An effective method for estimating MD involves minimizing the squared function of the difference between the measured magnetic flux density and the flux density derived from the numerical analysis of the region surrounding the permanent magnet. However, as there are countless MDs for reconstructing the specified magnetic flux density, it is important to estimate the realistic MD. In this article, we proposed an estimation method to derive the realistic MD based on the Biot–Savart law incorporated with mathematical programming (steepest descent method and quasi-Newton method). In this method, MD was formulated in the spherical coordinate system. Since the constraints regarding magnetization intensity and angle were defined by the sigmoid function, both quantities were definitely distributed within the constraints. Consequently, the generation of nonphysical MD was efficiently suppressed. Furthermore, the validity of the proposed method was illustrated in comparison with other estimation methods, such as the least-square method, Tikhonov regularization, and the truncated singular value decomposition (TSVD) method.