This paper deals with contact analysis of rolling bearings. A new mathematical model and numeric method are presented in this paper for the purpose of contact analysis of rolling bearings based on the principle of mathematical programming method. Three-dimensional (3D), finite element method (FEM) is introduced to calculate deformation influence coefficients and gaps of assumed pairs of contact points between the contact surfaces in the mathematical model. Special software is developed to realize numeric calculation procedures of the contact analysis. With the help of the developed software, loaded bearing contact analyses are conducted for a deep groove ball bearing and a cylindrical roller bearing. In the case of the ball bearing, it is found that calculated contact pressure on the ball surfaces is correct elliptical distribution (usually called contact ellipse). This result is more reasonable than the results obtained by commercial software SolidWorks and other researchers. In the case of the roller bearing, it is found that edge-loads (non-Hertz contact that cannot be analyzed with Hertz theory) on the roller surface are analyzed when the roller is not crowned. It is also found that the edge-loads disappear and contact pressure on the roller surface becomes uniform distribution when the roller is crowned with Johson-Gohar curve. Since the most important features (contact ellipse and edge-loads) of the bearing contacts are analyzed successfully by the developed software and these results cannot be obtained simply by general methods, the mathematical model and numeric method presented in this paper are of practical meaning in engineering design and the developed software can be used as a tool for contact strength and rigidity calculations of the rolling bearings.