The evaluation of the free energy of hydration depends essentially upon a calculation of the difference in polarization energy in a vacuum and in water. This is feasible by electrostatic continuum theory which is lucid and computationally effortless. It is, however, insufficient to apply the frequently used linear response approximation and nonlinear response effects must be integrated with the continuum model. We start from the time-honored Langevin-Debye theory and modify it to even describe polar liquids such as aqueous solvents. The modified Langevin-Debye model is then applied to study the hydration of cations. It is to be noted that reaction-field models are unapplicable to a nonlinear dielectric in an inhomogeneous field. The nonlinear effects are less than 10 kJ/mol or 2% for monovalent ions. However, for multiply charged ions the effect can be as large as 1000 kJ/mol or 15% and the nonlinear response model reduces the mean error of the calculated hydration free energy by at least 60%, which results in a mean error of only 4% for such ions. The mean error is of the same magnitude as the electrostatic free energy uncertainty.