In this paper we present results on the stability of perturbation methods for the evaluation of Dirichlet-Neumann operators (DNO) defined on domains that are viewed as complex deformations of a basic, simple geometry. In such cases, geometric perturbation methods, based on variations of the spatial domains of definition, have long been recognized to constitute efficient and accurate means for the approximation of DNO and, in fact, several numerical implementations have been previously proposed. Inspired by our recent analytical work, here we demonstrate that the convergence of these algorithms is, quite generally, limited by numerical instability. Indeed, we show that these standard perturbative methods for the calculation of DNO suffer from significant ill-conditioning which is manifest even for very smooth boundaries, and whose severity increases with boundary roughness. Moreover, and again motivated by our previous work, we introduce an alternative perturbative approach that we show to be numerically stable. This approach can be interpreted as a reformulation of classical perturbative algorithms (in suitable independent variables), and thus it allows for both direct comparison and the possibility of analytic continuation of the perturbation series. It can also be related to classical (preconditioned) spectral approaches and, as such, it retains, in finite arithmetic, the spectral convergence properties of classical perturbative methods, albeit at a higher computational cost (as it does not take advantage of possible dimensional reductions). Still, as we show, an alternative approach such as the one we propose may be mandated in cases where substantial information is contained in high-order harmonics and /or perturbation coefficients of the solution.