A simulation of the nonlinear electromechanical macroscopic behavior of ferroelectric materials by means of the finite element method is presented. A material point is depicted by a representative volume element, for which homogeneous boundary conditions are valid. The evolution of integral averages over the representative volume element is to homogenize the results. For this homogenization we favor a finite element model in which each Gauss point represents exactly one single crystal. Their number of internal variables is limited to the lattice orientation and the volume fractions of the domains. The former are randomly distributed in space. It is possible to calculate the material behavior for arbitrary coupled and nonlinear electromechanical loading cases, but the model is not effective for the solution of boundary value problems for entire bodies.