Computation of phase response curves via a direct method adapted to infinitesimal perturbations

A new numerical algorithm for computation of phase response curves of stable limit cycle oscillators is proposed. The idea of the algorithm originates from a direct method that is based on computation of the oscillator response to short finite pulses delivered at different phases of oscillations. Here we adapt the direct method to the case of infinitesimal perturbations and compare our algorithm with the standard algorithm based on the backward integration of the adjoint equations. In contrast to the standard algorithm, our algorithm does not require any backward integration and it is easier to program since a necessity of numerical interpolation for the Jacobian matrix is avoided. In addition, we demonstrate by examples that our algorithm is faster than the standard algorithm and this advantage is especially notable for weakly stable limit cycle oscillators.