We presented a multiscale nonlinear finite element simulation to analyze domain switching behaviors in ferroelectric materials. We utilized an incremental form of fundamental constitutive law to consider changes in the material properties caused by domain switching. A multiscale nonlinear problem was formulated by employing the asymptotic homogenization theory based on the perturbation method and implemented using finite element analysis. The developed simulation was applied to barium titanate with a Perovskite-type tetragonal crystal structure. The 90° and 180° domain switching behaviors of a single crystal were computed for verification. The nonlinear behaviors of a bulk polycrystal with virtual microstructure were analyzed as a case study. The variation of the crystal orientation distribution in the polycrystalline microstructure was analyzed to reveal its influence on macroscopic hysteresis and butterfly curves.