Magnetic monopoles are an inevitable feature of post-inflation symmetry-breaking phase transitions in grand unified theories. Analytic estimates of their density indicate that they are compatible with standard cosmology only if their mass is less than ${10}^{11}$ GeV [1, 2]. We initiate a programme of numerical studies of monopole dynamics by simulating a gas of ’t Hooft-Polyakov monopoles formed by the Kibble mechanism after a phase transition. In this paper we simulate monopoles in a radiation background, but without interactions with the radiation, in order to resolve differences between analytical models. We find that during the radiation era, the monopoles find each other and annihilate efficiently enough to keep their density fraction constant, which supports the modelling of Zel’dovich and Khlopov [1] and Preskill [2] in the epoch when plasma interactions can be neglected. In the matter era the density fraction decreases logarithmically. Further work is needed to quantify the effect of the thermal bath, which is expected to reduce the annihilation rate at later times.