Localized Wannier functions provide an efficient and intuitive framework to compute electric polarization from first-principles. They can also be used to represent the electronic systems at fixed electric field and to determine dielectric properties of insulating materials. Here we develop a Wannier-function-based formalism to perform first-principles calculations at fixed polarization. Such an approach allows to extract the polarization-energy landscape of a crystal and thus supports the theoretical investigation of polar materials. To facilitate the calculations, we implement a quasi-Newton method that simultaneously relaxes the internal coordinates and adjusts the electric field in crystals at fixed polarization. The method is applied to study the ferroelectric behavior of $mathrm{BaTiO_3}$ and $mathrm{PbTiO_3}$ in tetragonal phases. The physical processes driving the ferroelectricity of both compounds are examined thanks to the localized orbital picture offered by Wannier functions. Hence, changes in chemical bonding under ferroelectric distortion can be accurately visualized. The difference in the ferroelectric properties of $mathrm{BaTiO_3}$ and $mathrm{PbTiO_3}$ is highlighted. It can be traced back to the peculiarities of their electronic structures.