This paper proposes a new version (fundamentally different from the ex
isting ones) of finite element method for the mean and covariance func
tions of the displacement for bending beams with spatially random stif
fness. Apart from the conventional finite element method for stochasti
c problems, which utilizes either perturbation or series expansion tec
hnique or the Monte Carlo simulation, the present method is based on t
he newly established variational principles. The finite element scheme
is formulated directly with respect to the mean function and covarian
ce function, rather than perturbed components of the displacement. It
takes into account an information on joint probability distribution fu
nction of the random stiffness to obtain the covariance function of th
e displacement. Therefore, the accurate solution can be obtained even
if the coefficient of variation of the random stiffness is large, in c
ontrast to conventional technique. Several examples are given to illus
trate the advantage of the proposed method, compared with the conventi
onal ones.