Numerical Approach to Painlev\'e Transcendents on Unbounded Domains

A multidomain spectral approach for Painlev\'e transcendents on unbounded domains is presented. This method is designed to study solutions determined uniquely by a, possibly divergent, asymptotic series valid near infinity in a sector and approximates the solution on straight lines lying entirely within said sector without the need of evaluating truncations of the series at any finite point. The accuracy of the method is illustrated for the example of the tritronqu\'ee solution to the Painlev\'e I equation.


Introduction
Painlevé transcendents appear in many applications, see for instance [5], and have therefore been included as a chapter in the NIST digital library of mathematical functions [27], thus complementing classical transcendents as hypergeometric functions. Consequently, the numerical computation of these transcendents has seen a great deal of activity over the years. The goal of this paper is to offer an efficient numerical approach for transcendents characterized by their asymptotic behavior. The equations are studied on unbounded domains which offers a boundaryless approach to Painlevé transcendents.
Numerical approaches to Painlevé transcendents found in the literature either use the integrability of the equations and essentially solve Riemann-Hilbert problems [12], see [28] for numerical approaches, or evaluate Fredholm determinants [1,2], or they apply standard techniques for ordinary differential equations (ODEs). The latter allow also to study non-integrable deformations of the Painlevé equations with the same numerical techniques, e.g., to study characteristic features of integrable ODEs not shared by non-integrable ones. These approaches can loosely be split into techniques solving initial value problems on one hand, see, e.g., [21], and boundary value problems on the other, see, e.g., [8,10,13] and references therein. The boundaryless approach of the present paper is in particular intended to study the non-integrable deformations of Painlevé equations.
Many interesting transcendents are characterized by their asymptotic behaviour. The standard approach to numerically construct the transcendents is to provide a formal solution with the desired asymptotic in the form of an asymptotic series which is in general diverging. This series will be truncated according to standard rules, see, e.g., [15]. The truncated series will then provide boundary or initial data for the numerical solution at some finite value of the variable z in the complex z-plane. In this paper we address the obvious disadvantage of this approach, namely the use of diverging series and the arbitrariness in the choice of the value of z where to set up the initial or boundary value problem as well as in the truncation of the series. Instead, here we propose to study solutions determined by their asymptotics at infinity in some sector of the complex z-plane by numerically computing such solutions along a given straight line lying entirely within the sector with the use of the following techniques: • We use independent variable transformations to compactify the line into several finite intervals, such that the extreme endpoints of the collection of intervals correspond exactly to the points at infinity at the opposite ends of the straight line.
• We regularise the desired solution at infinity by subtracting sufficiently many terms from its large z asymptotic expansion so that the remainder tends to zero as z → ∞ along the line. Then in each interval with an endpoint corresponding to z → ∞ we solve for this normalised remainder rather than the Painlevé function itself. Of course, the remainder satisfies its own differential equation derivable explicitly from the Painlevé equation of interest.
The advantage of compactification is that we may use polynomial interpolation on finite intervals to numerically treat the boundary-value problem, and moreover working with the remainder gives the advantage that the boundary conditions are exactly known at the extreme endpoints of the collection of intervals. In this sense, our method is "boundary-less". Note that the equations studied near the infinite z-points are singular, and thus no boundary conditions will be applied there.
The example for which we illustrate the concept is the tritronquée solutions to the Painlevé I (PI) equation, which are characterized by Boutroux [3,4] by their asymptotic behaviour, Boutroux showed that there are solutions in four consecutive sectors of the complex plane of angle 2π/5 with the behaviour at infinity as in (2), which asymptotically are pole-free and which he called tritronquée. The sector can be chosen to be given by | arg z| < 4π/5 and we use the notation Ω(z) in the rest of the paper to refer to this well-defined particular solution of (1). There are four other tritronquée solutions related to this solution via Ω n (z) = e 4πin/5 Ω e 2πin/5 z , n = ±1, ±2, which corresponds to a general symmetry of the PI equation not limited to tritronquée solutions. Note that the considered tritronquée solution has the symmetryΩ(z) = Ω(z). It is shown in [20] that any solution to the PI equation is a meromorphic function on the complex plane. For a comprehensive discussion of tritronquée solutions see Joshi and Kitaev [21]. Kapaev [22,23,24] gave a complete characterization of the tritronquée solutions in terms of the Riemann-Hilbert problem associated to the PI equation.
In [10] these solutions are applied to describe asymptotically the behaviour of solutions to the nonlinear Schrödinger (NLS) equation in the semiclassical limit near a critical point. It was conjectured that the tritronquée solutions are pole-free in the whole sector | arg z| < 4π/5, not just asymptotically, which was then proven in [6]. In this paper we are interested in the tritronquée solutions in the regular sector and will construct them on unbounded lines via polynomial interpolation. If poles are to be studied, approaches based on Padé approximants as in [11,13,26] are to be used.
The paper is organized as follows: in Section 2 we present the numerical techniques we employ. In Section 3 we apply a multi-domain spectral approach to the tritronquée solution on unbounded lines. We add some concluding remarks in Section 4.

Multidomain spectral method
In this section we briefly summarize the numerical methods we use. We are interested in solutions on lines z = ax + b in the complex plane, where a, b ∈ C are constant with respect to x ∈ R. The line will be divided into several intervals for x each of which will be mapped to the interval [−1, 1]. On the latter, polynomial interpolation is used to approximate the solution.
As discussed for instance in [13], it is numerically problematic to set up an initial value problem for the tritronquée solution for large |z|, say for x 0 1 and to integrate the solution up to the value −x 0 since cancellation errors will lead to a loss of digits which can destabilize such a shooting approach. In practice it is numerically more stable to set up a boundary value problem with boundary data at ±x 0 .

Domains and equations
The line z = ax + b, x ∈ R will be divided into three domains, x < x l , x l < x < x r and x > x r numbered I, II, III respectively; here x l < x r are real constants.
The tritronquée solution we are interested in is characterized by the asymptotic behaviour (2). This implies, see [21], that near infinity the solution can be represented by a formal series, The series in (3) is divergent and therefore of little use from an analytical point of view. However, its form suggests a numerical approach to solve for the function This function is regular at infinity in the sector under consideration: v → 0 as z → ∞ for | arg z| < 4π/5. The PI equation (1) implies the following equation for v, Note that the same approach can be used for similar Painlevé transcendents though exponentially small terms in the asymptotics might require the use of more general transformations than considered here. Furthermore, starting with a function that is singular at infinity there are other ways to generate a function that vanishes there, e.g., in the considered case one could take 1/Ω(z) that would have a prescribed algebraic decay in s at infinity. This approach will not be explored here. This equation will be solved in domains I and III, equation (1) in domain II. Each of these domains will be mapped to the interval [−1, 1] in the following way (l ∈ [−1, 1]), in domain II and in domain III. Note that the choice of x l and x r can be optimized according to [9] (we mainly check this via the Chebyshev coefficients as shown in the next section), and that it is straight forward to generalize the approach to more than 3 domains. This was just chosen here to keep the presentation simple. At the domain boundaries, the function Ω has to be C 1 in x. Together with the PI equation, this guarantees that the solution will be smooth on the whole considered line x ∈ R. Since equation (5) is singular for s = 0, i.e., for z → ∞, no condition needs to be given there (the numerical approach will automatically produce the regular solution at this point).
Remark 1. All roots appearing in this paper are to be understood as being defined on a twosheeted Riemann surface. The sign of the roots will be fixed as usual at some point which is not a branch point. By analytic continuation as for instance in [14], sheets on the Riemann surface can be numerically established independently of the definition of the root in Matlab.

Polynomial interpolation and τ -method
Since each of the above three domains have been mapped to the interval [−1, 1], the task is now to approximate a smooth function on this interval. A standard approach is to choose a suitable discretisation of the independent variable l, see for instance [8,29,31], and to use polynomial interpolation on the collocation points l j , j = 0, . . . , N c . This leads to an approximation of the derivatives in terms of so-called differentiation matrices obtained by differentiating the interpolation polynomials. In [17], a collocation method with cubic splines distributed as bvp4 with Matlab was applied. In [18], a Chebyshev collocation method on Chebyshev collocation points l j = cos(jπ/N c ), j = 0, . . . , N c was used. The latter is related to an expansion of the solution in terms of Chebyshev polynomials T n (l) = cos(n arccos(l)), n = 0, 1, . . . of some function u(l), Putting y = arccos(l), a 0 = c 0 , and a n = a −n = c n for n = 1, 2, . . . , N c , one gets for relation (6) u(l) ≈ Nc n=−Nc a n e iny , i.e., a discrete Fourier transform in the variable y with collocation points y j = πj/N c , j = −N c , . . . , N c . Thus one advantage of the use of Chebyshev polynomials is that the expansion in terms of the latter can be efficiently computed with a fast cosine transform (FCT) which is related via (7) to a fast Fourier transform (FFT) for which efficient algorithms exist (note that special care has to be taken of the terms j = ±N c in order to establish the relation to an FFT, see the discussion in [29]).
A further advantage of the relation between Chebychev and Fourier series is that results for the decrease of the Fourier coefficients of a given function f for large n in dependence of the regularity of f can be applied directly to Chebychev series, see [29,30]. Since Fourier (6) and Chebychev sums (7) can be seen as truncated series, the numerical error in approximating a given function by a sum is of the order of the coefficients c Nc and a Nc respectively. We recall from [29,30] that for a function f ∈ C p−1 (T) with a pth derivative of bounded variation, the Fourier coefficient a n = O n −(p+1) for n → ∞. For a function f ∈ C ∞ (T), the Fourier coefficients decrease as a n = O(n −m ), ∀ m ∈ N. In numerical applications where one works with finite precision, this means that C ∞ functions are approximated with an essentially exponentially small error. Note that because of the relation between Chebychev sums and discrete Fourier transforms, functions u ∈ C ∞ ([−1, 1]) are approximated with essentially exponential accuracy by the sum (6). This allows to control the numerical resolution via the spectral coefficients.
Remark 2. The above description of properties of Chebychev sums motivates our choice of the independent variable s = 1/ √ z for the remainder function v defined in intervals I and III when we compute the Painlevé-I tritronquée solution Ω(z). Indeed, the asymptotic expansion (3) implies that v(s) has an asymptotic expansion as s → 0 with | arg(s)| < 2π/5, and this expansion is a (divergent) power series in s. Since Ω(z) is known to be analytic in the sector | arg(z)| < 4π/5, it follows that v(s) is analytic for |s| > 0 in the sector | arg(s)| < 2π/5, and then the Cauchy integral formula shows that the asymptotic power series for v(s) is differentiable term-by-term, i.e., arbitrary derivatives of v(s) have asymptotic power series expansions as s → 0 in the indicated sector given by formal derivatives of the power series for v. In particular, this means that if the line z = ax + b is taken with | arg(a)| < 4π/5, all derivatives of v(s) will be continuous up to the endpoints when intervals I and III are mapped to [−1, 1], which guarantees that the Chebychev coefficients will satisfy c n = O(n −m ), ∀ m ∈ N as n → ∞. However if | arg(a)| equals or exceeds 4π/5 we can draw no such conclusion from the asymptotic series (3).
The PI equation (1) is thus replaced by N I + 1, N II + 1 and N III + 1 algebraic equations where N I , N II , N III give the number of collocation points on the respective domain. The resulting system of algebraic equations is solved using Newton's method. As the initial iterate we use in domains I and III v = 0; for domain II, we compute the first 6 terms of the asymptotic series (3) for x l and x r and use the linear interpolate between these values as an initial iterate.
The normal Newton iteration for the solution of an equation F (u) = 0 takes the form where Jac F is the Jacobian of F and u n the nth iterate. The accessible precision is mainly limited by the conditioning of the Chebyshev differentiation matrices which is of the order of N 2 c , see the discussion in [29] for second order differentiation matrices. Note that this problem can be addressed as mentioned by introducing more than 3 domains. The iteration is stopped typically at a residual ||F (u n )|| ∞ of 10 −10 .
The junction conditions at the domain boundaries are implemented via Lanczos' τ -method [25]. This means that the C 1 conditions are simply replacing the lines corresponding to x = x l and x = x r in the equation Jac F ((u n ))(u n+1 − u n ) = −F (u n ), before numerically solving for u n+1 . The derivatives in the differentiability conditions on u at the boundary are again computed via Chebyshev differentiation matrices. It is known, see [29], that the boundary conditions implemented in this way will be satisfied with the same spectral accuracy as the solution.

Tritronquée solutions on lines in the complex plane
In this section we consider two examples for the above numerical scheme, the tritronquée solution to PI on the imaginary axis and a line close to the Stokes' lines arg z = ±4π/5 where the solution is known to show an oscillatory behavior, see Kapaev [24]. On the Stokes' lines, the tritronquée solutions have an oscillatory singularity at infinity, and the spectral approach will obviously have problems to approximate such a behavior.

Imaginary axis
As a first example, we will study the case z = ix. We use x r = −x l = 10 and N I = N III = 20 and N II = 256 collocations points in the respective domains. The solution Ω is shown for |x| < 15 in Fig. 1. It can be recognized that the solution is smooth on the whole interval and in particular at the domain boundaries x l , x r . The symmetryΩ(z) = Ω(z) of the solution is obvious. The function Ω in domains I and III follows from the computed functions v via (4). The computed v in dependence of s, i.e., on the unbounded domains, is shown on the left of Fig. 2. It can be seen that the functions are of the order of 10 −4 and smooth on the considered intervals. As expected they vanish at infinity. This shows that the tritronquée solution is in the considered example already very close to the asymptotic solution.
As already stated, an expansion in terms of Chebyshev polynomials offers the possibility to check the numerical resolution via the decrease of the spectral coefficients c n (6) as computed by an FCT. The coefficients in domain II are shown on the right of Fig. 2. They decrease to the order of the rounding error for n ∼ 128. Therefore it would have been possible to use just half of the value of N II without loss of accuracy. The corresponding coefficients for domain I and III can be seen on the left and right respectively of Fig. 3. Maximal accuracy is reached in this case with roughly 15 polynomials. The spectral coefficients indicate that the solution is as expected smooth on the whole line.

Line close to the Stokes' lines
The situation becomes numerically more demanding close to the Stokes' lines arg z = ±4π/5 because of the oscillatory behavior of the solutions on the latter. We consider the case z = exp(i(4π/5 − 0.05))x with x ∈ R. Again we use x r = −x l = 10 and N I = 20, but this time N II = N III = 256 to provide more resolution near the oscillatory singularity. Though the solution is still obeying the asymptotics (2) on this line, the vicinity to the Stokes' line produces oscillations with decreasing amplitude towards infinity as can be seen in Fig. 4. Note that there are no oscillations for x < 0. Again the solutions are smooth on the whole shown interval and thus also at the domain boundaries. The actually computed functions v in domains I and III can be seen in Fig. 5 on the left. In domain I the solution is once more very close to the asymptotic behaviour, the solution is of the order of 10 −4 . However, this is not the case in domain III where the deviation near x r is of the order 0.04.
The spectral coefficients c n (6) in domain II can be seen on the right of Fig. 5. Due to the oscillatory nature of the solution, 256 Chebyshev polynomials are needed to reach the level of the rounding error. The spectral coefficients for domains I and III are shown in Fig. 6. As expected the saturation level is reached in domain I with roughly 15 polynomials. But the vicinity of a singularity at infinity on the Stokes' line leads to a much smaller decrease of the coefficients in domain III. Though the solution in Fig. 5 appears rather innocent, with 256 polynomials just a resolution of the order of 10 −6 is reached. Note that this does not change much if a larger x r is chosen for the same N III which increases the numerical resolution. This is a clear hint on

Outlook
In this note we have shown for the example of the tritronquée solutions of PI that Painlevé transcendents can be computed on unbounded domains even if they have an algebraic increase in a local parameter near infinity which can be the branch point of a Riemann surface. A related example would be the tritronquée solutions of the second equation in the Painlevé hierarchy studied in [16]. If a Painlevé transcendent decreases exponentially as the Hastings-McLeod solution [19] u(x) of the Painlevé II equation on the real line for x → +∞, it is of course sufficient to impose the Dirichlet condition u(x r ) = 0 for a sufficiently large x r such that u vanishes there with numerical accuracy. However, the asymptotic behaviour of the same solution for x → −∞ is given by u ∼ √ −x. Thus the approach illustrated for the tritronquée solutions could be also applied to the region x < 0 of the Hastings-McLeod solution.
Since Painlevé transcendents can be extended to meromorphic functions in the complex plane, it would be interesting to extend the approach to whole domains of C. As was pointed out in [10], it is in this case better to set up a boundary value problem for the Laplace equation where the solution is known to be holomorphic since it is known to be also harmonic there. Since the Laplace equation is linear, no computationally expensive iteration will be needed in this case. In the case of the tritronquée solutions one could think of computing the solution on lines close to the Stokes' lines as in the previous section and use polar coordinates for the Laplace equation to obtain the solution on part of the regular sector. As was mentioned already in [10], the coordinate singularity at the origin leads to a loss of accuracy there. Such problems can be avoided if as in [7] for the hypergeometric equation elliptic or rectangular domains are used. It would be also interesting to explore the sector with poles in this way by combining the treatment of the asymptotics of this paper with the use of Padé approximants as in [11,13]. This will be the subject of further research.