To simulate the polarization dynamics of ferroelectrics, the model of local fields is improved by considering thermal vibrations. The local electric field at a dipole is the superposition of the external field, the fields of all other charges and dipoles and the fields of the image charges at the electrodes. The model is based on a sequence of single dipole flips which are thermally activated. The time to flip a single dipole depends on its deterministic transition rate depending on the local field and on a probabilistic factor. In each step, the dipole with the shortest flip time is switched. Thermal vibrations change the distances between the dipoles. The variation of distances modifies the local fields at the dipoles. The extended model considers these variations by multiplying the local fields in each step with Gaussian distributed random numbers. The model is applied to compute the polarization of two and of three dimensional systems based on the barium titanate structure. Polarization switching, hysteresis loops and the decay of the polarization during heating are simulated.The simulations yield intrinsic dead layers close to the electrodes and around defects which cannot be switched even in very strong fields. These nonswitchable layers are nuclei for domains and thus nuclei for polarization switching. The switching time of the system vastly decreases with the amplitude of the thermal vibrations. Moreover, the thermal vibrations enable the polarization switching in low external fields, decrease the coercive fields and cause a sharp phase transition from ferroelectric to paraelectric.