A numerical methodology for the Painlevé equations

https://doi.org/10.1016/j.jcp.2011.04.007Get rights and content

Abstract

The six Painlevé transcendents PI  PVI have both applications and analytic properties that make them stand out from most other classes of special functions. Although they have been the subject of extensive theoretical investigations for about a century, they still have a reputation for being numerically challenging. In particular, their extensive pole fields in the complex plane have often been perceived as ‘numerical mine fields’. In the present work, we note that the Painlevé property in fact provides the opportunity for very fast and accurate numerical solutions throughout such fields. When combining a Taylor/Padé-based ODE initial value solver for the pole fields with a boundary value solver for smooth regions, numerical solutions become available across the full complex plane. We focus here on the numerical methodology, and illustrate it for the PI equation. In later studies, we will concentrate on mathematical aspects of both the PI and the higher Painlevé transcendents.

Introduction

The six Painlevé transcendents were discovered about 100 years ago [4], [13], [23], [24], and have since received extensive attention. They satisfy nonlinear second order ODEs, cannot be reduced to other types of special functions, and obey the Painlevé property of all its solutions being free from movable branch points (however, the solutions may have movable poles or movable isolated essential singularities) [18], [21]. They arise in many branches of physics, including classical and quantum spin models, scattering theory, self-similar solutions to nonlinear dispersive wave PDEs, string theory and random matrix theory. For a summary of their applications, with extensive references, we refer to [7] and to Sections 32.13–32.16 in the NIST Handbook of Mathematical Functions [21].

The attention that the Painlevé functions have received over the last century has primarily been limited to their analytic properties and applications. The numerical computation of these functions, by contrast, has received comparatively little attention. With the exceptions of [8], [20], [22] (discussed later), numerical studies have so far mostly focused on pole-free intervals along the real line [17], [21], [25] or on a pole-free region in the complex plane [11]. A possible explanation for this state of affairs is that the large pole fields in the complex plane associated with the Painlevé ODEs may be thought of as a challenge to standard numerical methods for the integration of ODEs.

We argue in the present work that these pole fields provide, in fact, especially advantageous integration opportunities. A ‘pole friendly’ algorithm based on Padé approximation is devised that is capable of integrating accurately and stably through pole fields for distances in the 105  106 range. This allows numerical studies not only of ‘near-fields’, but also of transition processes towards ‘far fields’, for which an extensive body of asymptotic work is available [16], [21]. When the pole fields do not cover the full region of interest in the complex plane, pole-free parts can be obtained by solving ODE boundary value problems (BVPs).

The first numerical initial value problem (IVP) algorithm appropriate for dealing with poles was given by Willers in 1974 [30]. It is here generalized from a single integration path along the real axis to networks of paths throughout regions of the complex plane. The present study is focused on numerical issues, and uses for illustrations the PI equationd2udz2=6u2+z.In forthcoming studies, we will change the focus towards exploring the solution properties of both the PI and the higher Painlevé transcendents.

As a background to the present work, we note first a few analytical results for the PI function. It follows from (1.1) that the ODE is left invariant by the variable changes u  ω3u, z  ωz when ω5 = 1. Related to this result, the pattern of solutions for the ODE form (at sufficiently large distances from the origin) five distinct asymptotic sectors in the complex plane, each of angular width 2π/5. Assuming that u(0) and u′(0) are both real (i.e. u(z) is real-valued along the real axis), the solutions can form pole patterns similar to the very schematic illustrations in Fig. 1.1. Further analytic motivations for these sectors, including their relations to the tronquée (truncated) cases and the regularities within the pole fields, go back to [4], and are also summarized in Section 12.3 of [16]. Within smooth areas, the left-hand side (LHS) of (1.1) becomes negligible for increasing ∣z∣, and the two terms in the RHS will balance, leading to the estimateu(z)=±-z6+o(1).More detailed asymptotic expansions were noted already in [4]. However, (1.2) suffices for the present work. In the vicinity of a singularity, local analysis yieldsu(z)=1(z-z0)2-z010(z-z0)2-16(z-z0)3+O(z-z0)4,(with a second parameter, beyond the pole location z0, entering at the order O(z  z0)4). All the singularities are therefore double poles of equal strength and with zero residue.

In the remaining sections of this work, we first introduce the numerical pole field solver. After combining it with a boundary value method, we generate a collection of illustrations displaying some different types of solutions to PI. This is followed by a discussion of some numerical implementation variations, and by some concluding remarks.

Section snippets

Padé method for solving ODE IVPs

One of the key components in the present work is an IVP solver that is able to integrate accurately throughout pole fields of any extent (including right up to poles at unknown locations). We first introduce the relevant issues in the case of the standard first order ODEy(t)=f(t,y(t)),y(t0)=y0.The most basic numerical scheme for (2.1) is known as Forward Euler. This first order accurate method is obtained by considering only the linear part of a Taylor expansion, i.e.y(t+h)=y(t)+hf(t,y(t)).

Computational aspects of the Padé method when applied to 2-D pole fields

Consider again the PI functions, defined by (1.1). Fig. 3.1 shows a typical example of a computed solution over a region surrounding the point where the initial conditions (ICs) have been provided – here u(0) = −0.1875 and u′(0) = 0.3049. Since these values are close but not exactly equal to those producing the tritronquée case (cf. Fig. 1.1), there are pole fields in all the five main sectors. However, in four of the sectors, the pole fields have moved some distance away from the origin. This is

Some further illustrations of PI solutions

In this section we test the performance and reliability of our method against some important special cases of the PI equation.

Discussion of implementation options

Several choices need to be made when converting a numerical concept into efficient code. Some of these (such as the order of accuracy in the Padé-based integrator) relate mostly to the numerical algorithm itself, whereas other choices depend more on the software system that is used for the numerical implementation.

Conclusions

In this study we have demonstrated that the pole fields of Painlevé-type equations offer excellent opportunities for accurate and computationally very fast numerical solutions. When supplemented with appropriate BVP solvers, accurate and fast computational solutions become available also across pole-free areas. Although we have here only presented results for the PI equation, the numerical approach has been found to work equally well in preliminary computations with the PII equation. We have

Acknowledgements

The work of Bengt Fornberg was supported by the NSF Grants DMS-0611681 and DMS-0914647. Part of it was carried out in the fall of 2010 while he was an Oliver Smithies Lecturer at Balliol College, Oxford, and also a Visiting Fellow at OCCAM (Oxford Centre for Collaborative Applied Mathematics). The latter is supported by Award No. KUK-C1-013-04 to the University of Oxford, UK, by King Abdullah University of Science and Technology (KAUST). André Weideman was supported by the National Research

References (31)

  • J.R. Dormand et al.

    Runge–Kutta–Nystrom triples

    Comput. Math. Appl.

    (1987)
  • P.R. Graves-Morris et al.

    The epsilon algorithm and related topics

    J. Comput. Appl. Math.

    (2000)
  • Y.L. Luke

    Computations of coefficients in the polynomials of Padé approximations by solving systems of linear equations

    J. Comput. Appl. Math.

    (1980)
  • G.A. Baker et al.

    Padé approximants

    (1996)
  • D. Barton et al.

    An implementation of the Taylor series method for ordinary differential equations

    Comput. J.

    (1971)
  • J.P. Berrut et al.

    Barycentric Lagrange interpolation

    SIAM Rev.

    (2004)
  • P. Boutroux

    Recherches sur les transcendantes de M. Painlevé et l’étude asymptotique des équations différentielles du second ordre

    Ann. École Norm.

    (1913)
  • O.P. Bruno et al.

    Approximation of analytic functions: a method of enhanced convergence

    Math. Comp.

    (1994)
  • M.P. Calvo et al.

    High-order symplectic Runge–Kutta–Nyström methods

    SIAM J. Sci. Comput.

    (1993)
  • P.A. Clarkson

    Painlevé Equations – Nonlinear Special Functions

    (2006)
  • G.F. Corliss

    Integrating ODE’s in the complex plane - Pole vaulting

    Math. Comput.

    (1980)
  • J.R. Dormand et al.

    High-order embedded Runge–Kutta–Nystrom formulae

    IMA J. Num. Anal.

    (1987)
  • B. Dubrovin et al.

    On universality of critical behaviors in the focusing nonlinear Schrödinger equation elliptic umbilic catastrophe and the tritronquée solution to the Painlevé-I equation

    J. Nonlinear Sci.

    (2009)
  • B. Fornberg

    A Practical Guide to Pseudospectral Methods

    (1996)
  • B. Gambier

    Sur les équationes différentielles du second ordre et du premier degré don’t l’intégral générale est à points critique fixés

    Acta Math.

    (1910)
  • Cited by (0)

    View full text