Numerical methods of resonant dynamics with applications to the Galaxy are considered in this thesis. We derive generating functions for first-order perturbation theory and the associated orbital frequencies by matrix calculus. For two action-angle spaces (J,θ) and (i,φ) related by a canonical map I·φ+s, we show that J can be averaged over ergodic orbits φ to provide an estimator of I to within O(|s|²). We provide examples in one and two dimensions and compare the technique to calculations of actions by numerical line integration in Poincaré sections. We then use spectral dynamics and the Laskar frequency map (Laskar, 1993) to identify the dynamically important resonances of the 'flattened' axisymmetric isochrone potential. We simulate resonant capture in a low-order resonance by populating representative tori of a spherical isochrone Hamiltonian and integrating the orbits while adiabatically introducing axisymmetry. We use the averaging technique described above to observe the fraction of orbits captured, and we compare the result to a theoretical prediction. We return to first-order perturbation theory to analyse its strengths and weaknesses, in particular near orbital pericentre, and when one action is significantly smaller than another. We also reproduce the expected pendulum dynamics in the resonant action-angle plane for orbits in our capture simulation. We develop the concept of adaptive dynamics: we vary the initial orbital energy of the particles in the capture simulations and show that resonant and non-resonant orbits can be identified as clusters in the perturbed action plane. For a given Hamiltonian, we use the perturbed frequencies and a linear regression fit in the action plane as diagnostics of a set of model Hamiltonians on a grid in a suitable parameter space. We find we are able to constrain the parameters of a model Hamiltonian by this method. Finally, we reject the null hypothesis that resonant structures in phase space can be found by traditional methods of density estimation.