  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.
