Abstract: In natural element method,the trial and test functions were constructed with the natural neighbor interpolation method.The interpolation was based on the Voronoi tessellation of the scattered nodes in the problem domain.The integration of the weak form was performed in the Delaunay triangles which were the dual diagram of the Voronoi tessellation when the Galerkin method was used to form the discrete system equation.But there was obvious error in the numerical integration of the natural element method due to the characteristics of natural neighbor interpolation function.Every factor that produced the error possibly was analyzed,and a method using Monte Carlo integration was proposed to solve this problem.The weight coefficient was directly related to the precision,and its determination was also simple and effective.Integral point was cast in the Delaunay triangle,so the result of probability integral was close to mathematical expectation.The definite method of the least integral points,was given so as to enhance computational efficiency of the Monte Carlo integral as far as possible.Finally,the numerical example of the patch experiment and the cantilever beam confirmed the validity and feasibility of using Monte Carlo integration to solve these errors.