Abstract The solution of the equations of motion defined by the Hamiltonian function H = H o (q k ,p k ) + ∑ n=1 λ n H n (q k ,p k ) = E, k = 1,…,N is found through the use of power series expansions in the perturbation parameter λ. The solution is in the form of 2N independent integrals of motion j k = P k + ∑ n=1 λ n J kn (Q l ,P l ) = const W k = Q k + ∑ n=1 λ n W kn (Q l ,P l ) = δ kN ·t + const where the perturbation terms of these integrals are given by recursion formulas J kn = ∫ C Q N δH n δQ k − ∑ s=1 n−1 {J ks , H n−s } dQ′ N for k≠N, J N n = H n (Q l ,P l ) W kn = ∫ C Q N δH n δP k + ∑ s=1 n−1 {W ks , H n−s } dQ′ N and Pk(ql, pl), Qk(ql, pl) is the canonical transformation for which PN = H0(ql, pl), (δkN = Kronecker delta, {} = Poisson bracket). Quantities Jk, wk are the canonical variables giving the perturbed Hamiltonian the form H = JN. There are 2N − 1 time-independent integrals: J1 ,…, JN, W1 ,…, WN−1 describing the trajectory in 2N-dimensional phase space. The last integral WN − t = const, determines the position on the trajectory at any given time.