An efficient Newton-Picard algorithm for computing periodic solutions
of autonomous partial differential equations
D. Roose and K. Lust
Dept. of Computer Science, K.U.Leuven,
Celestijnenlaan 200A, B-3001 Leuven
Belgium
We present an algorithm to compute and analyse periodic solutions of certain large-scale dynamical systems represented by partial differential equations (PDEs). Nowadays it is well known that many dissipative systems only have low-dimensional long-term dynamics. Our method tries to exploit the low-dimensional behaviour by preserving traditional space- and time discretisation techniques, but adapting the numerical procedure.
We developed a hybrid Newton-Picard scheme based on the condensed Newton - supported Picard method of Jarausch and Mackens and on the recursive projection method of Shroff and Keller [3]. We combine straightforward time integration and single Newton-based shooting. The shooting method is only used in the low-dimensional subspace where the limit cycle is weakly attracting or repelling, while time integration is used in the (high-dimensional) orthogonal complement of that subspace. The low-dimensional subspace can be computed by iterative methods, which do not need the explicit building of the monodromy matrix. The method is particularly suited to be used in conjunction with continuation. We then have good starting values for both the limit cycle and the low-dimensional subspace basis. Since the dominant Floquet multipliers are a natural by-product of the algorithm, the method is also very well suited to perform a bifurcation analysis. We developed variants of the method for the fixed parameter case and for pseudo-arclength continuation. The different variants have different trade-offs with respect to convergence speed, convergence region, cost per iteration step and inherent parallelism. For the fixed parameter variants, it is easy to show that the asymptotic convergence rate is fully determined by the largest eigenvalue of the monodromy matrix restricted to the high-dimensional subspace. The size of the low-dimensional subspace is independent of the number of degrees of freedom of the space discretisation. The number of time integrations needed to decrease the residual with a given amount is thus independent of the space discretisation. More information can be found in [1] and [2].
During the presentation, we will shortly introduce the different variants and then elaborate on the convergence properties, both from the theoretical viewpoint and as observed in implementations. The results will be illustrated with some testcases.