Quantitative interpretations of the large Coulombic effects of changes in cation concentrations on processes involving polyanionic DNA require accurate theoretical descriptions of the thermodynamic consequences of cation-DNA interactions. In the present study, the thermodynamic consequences of accumulation of divalent and/or univalent cations in the vicinity of polyionic double-stranded DNA and of exclusion of univalent co-ions are characterized by ion-polyion preferential interaction coefficients Gamma(i) (i = 2+, +, or -). These are calculated using integrals of ion distributions generated from either canonical Monte Carlo (CMC) simulations or numerical solutions of the cylindrical Poisson-Boltzmann (PB) equation, for the same minimally parameterized cylindrical cell model over experimentally relevant ranges and ratios of the uni-and/or divalent cations. For solutions containing both types of cations, trends in Gamma(i)(MC) and in Gamma(i)(PB) With changes in the absolute and relative values of the divalent and univalent cation concentrations are examined and compared with trends calculated for solutions containing only one type of cation. Differences between Gamma(i)(MC) and Gamma(i)(PB) are quantified and related to differences between MC and PB predictions of the extents of local cation accumulation within 3 Angstrom of the polyion surface. Discrepancies between Gamma(2+)(MC) and Gamma(2+)(PB), and between Gamma-(MC) and Gamma-(PB), are significant whether or not univalent cations are present, but the difference Delta Gamma(2+) = Gamma(2+)(MC) - Gamma(2+)(PB) is relatively insensitive to changes in the concentration of salt (excess 1:1). Therefore, PB calculations may provide a satisfactory alternative to more computationally demanding MC simulations as a basis for analyzing the salt-concentration dependences of Gamma(2+) and of the closely related measurable thermodynamic properties that reflect the importance of Coulombic interactions.