We propose a new approach to the numerical solution of radiative transfer equations with certified
a posteriori error bounds. A key ingredient is the formulation of an iteration in a suitable (infinite dimensional) function
space that is guaranteed to converge with a fixed error reduction per step. The numerical scheme is based on approximately
realizing this outer iteration within dynamically updated accuracy tolerances that still ensure convergence to the exact solution.
On the one hand, since in the course of this iteration the global scattering operator is only applied, this avoids solving linear systems with densely populated system matrices while
only linear transport equations need to be solved. This, in turn, rests on a Discontinous Petrov—Galerkin scheme which
comes with rigorous a posteriori error bounds. These bounds are crucial for guaranteeing the convergence of the outer
iteration. Moreover, the application of the global (scattering) operator is accelerated through low-rank approximation and matrix
compression techniques. The theoretical findings are illustrated and complemented by numerical experiments with
a non-trivial scattering kernel.