Phase field methods have been widely used to study phase transitions and polarization switching in ferroelectric thin films. In this paper, we develop an efficient numerical scheme for the variational phase field model based on variational forms of the electrostatic energy and the relaxation dynamics of the polarization vector. The spatial discretization combines the Fourier spectral method with the finite difference method to handle three-dimensional mixed boundary conditions. It allows for an efficient semi-implicit discretization for the time integration of the relaxation dynamics. This method avoids explicitly solving the electrostatic equilibrium equation (a Poisson equation) and eliminates the use of associated Lagrange multipliers. We present several numerical examples including phase transitions and polarization switching processes to demonstrate the effectiveness of the proposed method.