The article describes the role of polarisation in calculation of multiple reflections. A mathematical model of multiple reflections based on the Stokes vector for beam description and Mueller matrices for description of surface properties is presented. On the basis of this model, the global illumination equation is generalised for the polarisation case and is resolved into volume integration. This allows us to obtain an expression for the Monte Carlo method local estimates and to use them for evaluation of light distribution in the scene with consideration of polarisation. The obtained mathematical model was implemented in the software environment using the example of a scene with its surfaces having both diffuse and regular components of reflection. The results presented in the article show that the calculation difference may reach 30 % when polarisation is taken into consideration as compared to standard modelling.