The goal is to develop a rigorous theory to study the effects of stochastic fluctuations of self amplified spontaneous emission free-electron laser (SASE FEL) on the near-resonant ionization of atomic/ionic systems. To this end, density matrix equations of motion are utilised in their raw form for the sake of Monte Carlo simulations and in their particular averaged forms. First, a couple of simpler averaging methods were used to quickly understand the general behavior of the interaction of fluctuating fields with atoms. To this end, the particular case of neon (Ne) and helium (He) were used, where the effects of the interplay of the pulse duration and the coherence time on the oscillations of the yield profile, which manifest due to Rabi oscillations, were studied in comparison to a coherent pulse. Second, a rigorous method of perturbative approach is developed, which is based on the expansion in terms of multitime cumulants which are a particular combination of the field’s coherence functions. The range of validity of the model is tested in terms of the field’s coherence temporal length and peak intensity. Particularly, the photoionization of helium (He) and lithium ion (Li+ ) via their doubly-excited state 2s2p 1 P has been studied with the interacting FEL’s 1st-order coherence function modelled as square-exponentially dependent. The traditional asymmetric resonant Fano-profile is broadened and is shown to acquire a Voigt profile. The effects of pulse duration, coherence time and peak intensity on the lineshape are clearly shown. The two approaches of averaging were compared and the supremacy of the latter is shown. When compared to Monte Carlo, it is seen that the rigorous averaging approach results in competing values up to higher peak intensities. This suggests that in the scenarios of low coherence time and peak intensities, the averaging method developed in this thesis is highly preferable over the traditional Monte Carlo given its high computational demands.