A formulation based on the finite element method and the periodic boundary condition for plane wave scattering by rough surfaces is derived and implemented. The method is applied to Monte Carlo simulations of one-dimensional random rough surface scattering with a Gaussian roughness spectrum satisfying Dirichlet boundary conditions. Convergence of the method is demonstrated by varying the number of sampling points. The percent error in conservation of energy is shown to be less than 0.007% for all the examples presented. The numerical solutions are compared with integral equation methods. The results from the finite element method are in perfect numerical agreement with the integral equation results for each surface realization if the periodic boundary condition is also used in the integral equation method. They are in excellent numerical agreement with a tapered wave solution of the integral equation for large surface lengths and upon averaging over many realizations. For the examples discussed, it is shown that the CPU time and memory storage requirements for the finite element method are much less than those of the integral equation method for cases when the number of horizontal sampling points is much larger than the number of vertical points in the region of discretization. Numerical results are presented for various rms surface heights, correlation lengths and surface lengths. The results are also shown to be consistent with the corresponding Kirchhoff approximation and the small perturbation theory within the domains of validity of the two approximate methods.