The formation of domain configuration in ferroelectrics is a consequence of energy minimization, and critically depends on their transformation strain and spontaneous polarization. In this article, we develop an energetic analysis on ferroelectric domain patterns using equivalent inclusion method, treating ferroelectric domain as an ellipsoidal inhomogeneous inclusion in a ferroelectric matrix. The potential energy of the domain is calculated in terms of its orientation and shape, and the energy minimizing configurations have been identified. Both tetragonal and rhombohedral crystals have been analyzed, and the lamellar domain configurations as predicted by the compatibility analysis have been recovered. Additional energy minimizing states have also been revealed, including needle type of domains and charged domains. Different contributions of strain compatibility and polarization compatibility have also been analyzed.