A Newton-Picard shooting method for computing periodic solutions of large-scale dynamical systems

A numerical method is presented for the efficient computation and continuation of periodic solutions of large systems of ordinary differential equations. The method is particularly useful when a shooting approach based on full Newton or quasi-Newton iteration is prohibitively expensive. Both stable and unstable solutions can be computed. The basic idea is to split the eigenspace of the monodromy matrix around the periodic orbit into two orthogonal components which contain, respectively, all eigenvalues with norm less than or greater than some threshold value less than unity. In the former ‘stable’ subspace one performs time integration, which is equivalent to a Picard iteration; in the much smaller, ‘unstable’ subspace a Newton step is carried out. A strategy for computing the stable and unstable subspaces without forming the full monodromy matrix is briefly discussed. Numerical results are presented for the one-dimensional Brusselator model. The Newton-Picard method with an appropriately chosen threshold value proves to be accurate and much more efficient than shooting based on full Newton iteration.