Based on a fast and accurate method for the stationary random seismic response analysis for discretized structures(Lin 1992, Lin et al. 1992), a Ritz method for dealing with such responses of continuous systems in developed. This method is studied quantitatively, using cantilever shear beams for simplicity and clarity. The process can be naturally extended to deal with various boundary conditions as well as non-uniform Bernoulli-Euler beams, or even Timoshenko beams. Algorithms for both proportionally and non-proportionally damped responses are described. For all of such damping cases, it is not necessary to solve for the natural vibrations of the beams. The solution procedure is very simple, and equally efficient for a white or a non-white ground excitation spectrum. Two examples are given where various power spectral density functions, variances, covariances and second spectral moments of displacement, internal force response, and their derivatives are calculated and analyses. Some Ritz solutions are compared with "exact" CQC solutions.