Rail transportation plays an important role in our everyday life, and there is fast development and modernization in the railway industry to meet the growing demand for swifter, safer and more comfortable trains. At the same time, the security of train operation and the maintenance of rails have to be considered. A lot of research on these issues has been carried out, among which the study of the contact between a train's wheel and the rail is particularly significant. The contact problem considers two elastic bodies. When they are pressed together, a contact area is formed where the two body surfaces coincide with each other. Moreover, an elastic field of stress, strain and displacement arises in each body. These stresses consist of normal stress (pressure) and frictional stress (traction) acting in the tangential direction. When solving the so-called {\it normal contact problem}, we search for the contact area and the pressure on it. The {\it tangential contact problem} is studied when the two bodies are brought into relative motion. If the relative velocity of the two surfaces is small, a creeping motion may be observed which is largely caused by the elastic deformation at the contact region. In those parts of the contact area where the tangential stress is small, the surfaces of the two bodies stick to each other. Otherwise, local relative sliding may occur. The research question is to find the adhesion and slip areas, and the tangential tractions. The solution methods for contact problems have been studied from the late nineteenth century, resulting in a variety of analytic and numerical approaches, w.r.t. their own specific applications. Motivated by the requirement of fast computation for involved applications such as the simulation of railway wheel-rail dynamics, we aim at developing fast numerical solvers for concentrated elastic contact problems in this thesis. Our work focuses on the contact between bodies of linear homogeneous elastic material. Moreover, it is a concentrated contact, i.e. the contact area is small compared to the dimensions of the contacting bodies. The models in use are provided by a variational formulation, which is based on a boundary element method (BEM). It gives rise to a convex optimization problem with linear or nonlinear constraints. The corresponding Karush-Kuhn-Tucker conditions provide the governing equations and contact conditions, that are numerically solved. The most time-consuming part attributes to solving a Fredholm integral of the first kind, resulting from the BEM. The corresponding Green's function expresses the relation between tractions and deformation, using a half-space approach. This integral yields linear systems with coefficient matrices that are dense, symmetric and positive definite. Moreover, they are Toeplitz in two-dimensional (2D) problems and block Toeplitz with Toeplitz blocks in three-dimensional (3D) problems. Fast computing techniques such as the fast Fourier transform (FFT) are explored. We start our work by solving the normal contact problem in Chapter 2. It is modeled by a linear complementarity problem, for which a full multigrid method (FMG) is presented. This method combines a multigrid (MG) method, an active set strategy and a nested iteration technique. It is applied to a Hertzian smooth contact and a rough surface contact. The results show the efficiency and robustness of the FMG method. Tangential contact is considered in Chapter 3 and Chapter 4. A 2D no-slip tangential problem is first studied in Chapter 3, where we mainly solve the surface integral. A fast MG method is proposed with an FFT smoother, where a Toeplitz preconditioner is constructed to approximate the inverse of the coefficient matrix. This smoother reduces many error components but enlarges some smooth error modes. Techniques such as subdomain deflation and row sum modification (RSM) are incorporated. Numerical experiments indicate rapid convergence and mesh-independence of MG with the FFT+RSM smoother. Moreover, FFT+RSM as a stand-alone solver also shows its efficiency. The complexity of these two methods is $\mathcal{O}(n\log (n))$, with $n$ the number of unknowns. We work on the 3D tangential contact in Chapter 4, where a nonlinear constrained optimization problem arises. A fast solver, called TangCG, is proposed. It combines an active set strategy and a nonlinear conjugate gradient method. The most pronounced component of this method is that it employs two types of variables in the adhesion and slip areas. Techniques including the FFT and diagonal preconditioning are also incorporated. The TangCG method is tested for Cattaneo shift problems, with different amounts of slip. It dramatically reduces the computational time, compared to the state-of-art ConvexGS method. The numerical methods presented above are based on the influence coefficients (ICs) that give the relation between tractions and deformation. In Chapter 5, we investigate ICs by computing them numerically. Based on a concentrated contact setting, an elastic model is built for this purpose and a finite element method (FEM) is employed. Suggestions about the FEM meshing and element types are given, considering the accuracy and computational cost. The effects of employing the numerical ICs on contact solutions are examined. The work in this chapter provides a guidance for developing fast solvers for conformal contact problems, which typically are governed by a larger and curved contact region. With the research presented in the present PhD thesis on efficient numerical solution techniques, the numerical solution of full-scale train-rail contact problems may have come one step closer. With the research presented in present PhD thesis, and with the resulting improved numerical solution techniques, it becomes one step closer to incorporate detailed contact models in the numerical simulation of rail vehicle dynamics and in the simulation of rail and wheel wear and track deterioration.