Isogeometric analysis (IGA) is a recently developed technology in computational mechanics. Its main idea is to use the same smooth and higher-order basis functions, e.g. non-uniform rational B-splines (NURBS), for the representation of both the geometry in CAD and the approximation of solution fields in analysis. Thus, the complex geometry of materials can be modeled exactly. As a consequence, geometric errors are eliminated. In addition, the high order continuities of the basis functions are ideally suitable for solving high order partial differential equations, like the equation of flexoelectric problems and Cahn-Hilliard equation, which require at least C¹ continuous approximations. These attributes permitted us to derive accurate, efficient, and geometrically flexible methods for the problems with higher-order derivatives. In this paper, a NURBS-based variational formulation for the phase field equations of ferroelectric materials is established. Several numerical examples are provided in this paper to confirm the accuracy of this method. After that, the polarization distributions inside the PbTiO₃ material are simulated with the help of IGA. Moreover, the polarization distribution inside a three-dimensional ferroelectric material with the consideration of flexoelectricity is obtained in this paper.