Particle accelerators are the basic instrument for the study of fundamental interactions. By increasing the energy we obtain a larger resolving power, by increasing the intensity a higher luminosity. High intensity low energy (1-10 Gev) machines have a wide variety of applications ranging from the production of neutron beams (after the collision with a target) to be used for spectroscopic analysis, to the treatment of nuclear wastes and, in a long range perspective, thermonuclear fusion, induced by the impact of heavy ion beams on pellets containing a mixture of hydrogen and tritium. Of high interest for research in biology and condensed matter is the synchrotron light produced with electron accelerators. In protons or ions accelerators of very high energy such as LHC,one of the key problems is is the control of nonlinear effects due to field imperfections in the super-conducting dipoles. Non linear forces due to multipolar components in the magnetic field determine, in the single particle phase space, a finite stability area, which must nevertheless by sufficiently wide to prevent unacceptable reductions of luminosity. The difficulty is due to the control of beam stability on a high number of turns (about 10^9), comparable with the number of revolutions of Jupiter from the origin of the solar system to present. Since the simulations exceed the available computational power, some strategies, inspired by celestial mechanics, have been proposed for the stability analysis. At first the normal forms methods was introduced to obtain analytically the dependence of the oscillation frequencies on the amplitude for the transverse betatron motion. It is a perturbative technique which applies to the one turn map, known as Poincare' map, which allows to describe the characteristics of the main nonlinear resonances. Another method, known as frequency map analysis, has also been applied. It consists in determining the frequencies and even the invariants of motion, from the orbits computed numerically. The degradation of the beam quality is also due to low frequency components and to rapid fluctuations of the current circulating in the magnets and also to mechanical micro-vibrations. This implies the presence of a slow variation in the parameters of the one turn map and in the presence of noise. The probabilistic description of the process is obtained from the neoclassic adiabatic theory and the theory of Hamiltonian systems with noise, which, in both cases, lead to a Fokker-Planck equation for the probability density function. For intense low energy beams the relevant problem is Coulomb repulsion between the charges, which leads to an increase of the transverse section of the beam and to a decrease of the linear restoring force. The problem may be treated at two different approximation levels. The mean field theory neglects the collisions between the particles and treats the single particle motion in the self consistent field generated by the remaining particles (Liouville's equation for the single particle coupled with Poisson's equation). The microscopic theory takes into account the soft and hard collisions by solving Hamilton's equations for a system on N particles of charge q. In the numerical computations N must be chosen much smaller (10^3 - 10^4) with respect to the physical value (10^11 - 10^12): the charge q is incresed correspondingly in order to leave the product qN unchanged. A relevant aspect is the study of the transient phase which leads to the thermodynamical equilibrium and the determination of the relaxation time. The existence of scaling laws which allow the extrapolation from the values of N attainable in the simulations to the physical values is crucial. Our group has first proposed the use of normal forms for the analysis of nonlinear betatronic motion and has contributed to a better theoretical understanding of the stability boundary and of the diffusive processes induced by slow modulation and noise. More recently we have developed particle in cell (PIC) numerical procedures to integrate the 2 and 3 dimensional mean field equations in order to investigate the effects of discretization errors. We have developed a 2 dimensional microscopic model to investigate the collisional effects. The main result is a scaling law for the dependence of the relaxation time on the particles number N, which was analytically obtained from Landau's equations and confirmed by the numerical computations performed on a parallel architecture.