My research is concerned with the analysis and numerical simulations of elliptic partial differential equations that model steady state flow (electric, thermal, fluid) in high contrast composite materials consisting of conducting and insulating inclusions that are close to touching. The coefficients in these equations vary rapidly, thus modeling the micro scale of the composites, and have large (even infinite) ratios of their maximum and minimum values. It is difficult to simulate numerically flow in high contrast composites because of singularities of the gradient of the solution between the high contrast inclusions. Solvers need fine meshes to resolve these singularities and lead to very large linear systems that are poorly conditioned. In this thesis, we first approximate the Dirichlet to Neumann (DtN) map for high contrast two phase composites. The mathematical formulation of the problem is to approximate the energy for an elliptic equation with arbitrary boundary conditions. The boundary conditions may have high oscillations, which makes our problem very interesting and challenging. Our main result is more than general homogenization of problems in high contrast composites because we consider the problem with arbitrary boundary conditions. In order to approximate the energy of the problem with arbitrary boundary conditions, we propose a method to divide the original problem into two subproblems in two separated subdomains. One subdomain is close to the boundary, i.e. the boundary layer, and the other subdomain is far from the boundary. We approximate the energy in these two subdomains separately and then combine them together to obtain the approximation in the whole domain. In the subdomain far from the boundary, the energy is not influenced that much by boundary conditions and methods are studied before. In the boundary layer, the energy strongly depends on the boundary conditions. We use a new method to approximate the energy there such that it works for arbitrary boundary conditions. We then directly apply the approximation of DtN map into numerical methods for solving problems in high contrast media. We use this approximation to construct preconditioners in nonoverlapping domain decomposition methods. Preconditioners constructed from the approximation of DtN map almost work as well as preconditioners from solving problems numerically, however it is much cheaper to construct preconditioners from theoretical approximation results. This leads to the idea of coupling theoretical results and numerical methods in order to save computational time for solving problems numerically in high contrast media.