{'fraction_non_alphanumeric': 0.0367504835589942, 'fraction_numerical': 0.024259908861423467, 'mean_word_length': 4.317991631799163, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Computing the solutions of the van der Pol equation to arbitrary precision

We describe an extension of the Taylor method for the numerical solution of ODEs that uses Pad\'e approximants to obtain extremely precise numerical results. The accuracy of the results is essentially limited only by the computer time and memory, provided that one works in arbitrary precision. In this method the stepsize is adjusted to achieve the desired accuracy (variable stepsize), while the order of the Taylor expansion can be either fixed or changed at each iteration (variable order). As an application, we have calculated the periodic solutions (limit cycle) of the van der Pol equation with an unprecedented accuracy for a large set of couplings (well beyond the values currently found in the literature) and we have used these numerical results to validate the asymptotic behavior of the period, of the amplitude and of the Lyapunov exponent reported in the literature. We have also used the numerical results to infer the formulas for the asymptotic behavior of the fast component of the period and of the maximum velocity, which have never been calculated before.


Introduction
The Taylor method is one of the oldest methods for the numerical solution of ODEs: in fact the well-known Euler method is essentially a Taylor method of order one. Higher order Taylor methods have been explored in the literature, but their use is still not widespread. From our (partial) exploration of the literature, we have found refs. [1,2,3,4,5,6,7,8] to be the most relevant to our discussion.
The paper by Deprit and Zahar [1] is possibly the first to show the advantages of using the Taylor method compared to a more standard Runge-Kutta-Nystrom algorithm: these authors found that the latter needs to use step size about 4000 smaller to reach the same accuracy of the former. Corliss and Chang [2,3] have applied the Taylor method to a series of problems (using fortran). The remaining papers, ref. [4] by Jorba and Zou and refs. [5,6,7,8] by Barrio and collaborators, are particularly important because of their emphasis on achieving high precision. In general it is found that the high-order Taylor method is very competitive with other numerical methods when high accuracy is needed; moreover this method can be used to calculate the numerical solutions to certain class of ODEs with unprecedented accuracy (from the abstract of ref. [7]: "this methodology is the only one capable of reaching precision up to thousands of digits for ODEs").
At the same time, Jorba and Zou clearly state a limitation to the Taylor method: "The main drawback is that the Taylor method is an explicit method, so it has all the limitations of these kind of schemes. For instance, it is not suitable for stiff systems.". A discussion of the applicability of Taylor method for stiff ODEs is also contained in [5] 1 .
The main goal of our paper is to describe an extension of the Taylor method that allows to obtain arbitrarily accurate numerical solutions even for stiff differential equations, such as the van der Pol equation that will be used in this paper as the ideal testing ground. The method that we propose is based on the use of Padé approximants rather than Taylor polynomials; we will show that the performance of the standard Taylor method is drastically improved both by enlarging the class of problems that can be attacked (e.g. stiff ODEs) and also by increasing the accuracy of the solutions with respect to the Taylor method with the same number of coefficients.
We will illustrate our method by applying it to the van der Pol oscillator, the most famous example of ODE supporting a limit cycle. For this problem there is an extensive literature, but surprisingly enough numerical solutions have never been calculated to a high level of accuracy, particularly for the stiff regime. Moreover the behavior of the period and amplitude of the limit cycle has been studied in depth both for small coupling (using perturbation theory, i.e. the Lindstedt-Poincaré method) and for large coupling (using asymptotic methods): by producing extremely accurate numerical solutions for a wide range of couplings we plan to show that it is possible to reproduce the known analytical results, but also to predict the asymptotic behavior of quantities not yet studied.
In our opinion, the quest for extremely precise numerical solutions is not a mere exercise of virtuosity but it has real applications: one cannot overstate the importance of predicting the behavior of a dynamical system while keeping under control the precision of the results. We believe that the technique described in this paper may be applied to a large number of problems for which accuracy is important.
The paper is organized as follows: in Section 2 we introduce our Padé -Taylor method (PTM) and discuss its main features; in Section 3 we consider the van der Pol oscillator and discuss the application of the Padé-Taylor method to it; in the subsections 3.1 and 3.2 we present the numerical results and apply appropriate extrapolation method to extract the leading asymptotic behavior for different quantities (particularly the period and the amplitude of the oscillator); finally, in Section 4 we draw our conclusions and discuss the possible directions of future work.

The Padé -Taylor method
Consider the ODE with initial conditions At t > t 0 the solution of eq. (1) can be approximated by a Taylor polynomial of order N The Taylor coefficients in eq. (3) can be obtained by repeated derivations of eq. (1) and by subsequently setting t = t 0 . Let T 0 be the radius of convergence of the power series obtained by taking the limit N → ∞: clearly we are justified to use the polynomial (3) only for t 0 < t ≪ t 0 + T 0 . Without going to the extreme case where the power series is divergent (T 0 = 0), the reader may recognize that a finite value of T 0 sets an upper limit to the stepsize in the numerical integration of eq. (1). As a time step τ 0 ≪ T 0 is chosen (in some way -we will discuss this point soon), a new approximation to the function at t 1 ≡ t 0 + τ 0 : from this approximate solution one can estimate the new Taylor coefficients corresponding to eq. (3) expanded arout t 1 and repeat the same arguments. In general we will have τ 1 = τ 0 (variable stepsize): additionally, at each step one can use polynomials of different order if needed (variable order).
Of course we run into trouble if the power series becomes divergent (vanishing radius of convergence) at any finite time t i : in this case τ i = 0 is the only allowed stepsize.
The natural choice to extend the range of applicability of the power expansion is to substitute the Taylor expansion with a Padé approximant of appropriate order: where P m (t − t 0 ) and Q n (t − t 0 ) are polynomials of order m and n and the expansion of eq. (4) reproduces the first N + 1 Taylor coefficients (m + n = N ). Let τ (0) i , with i = 1, . . . , n, be the (complex) roots of the polynomial Q n (τ ): the smallest positive real root (if any) will limit the region of applicability of (4) to smaller values of the stepsize.
Let us now how come to the important question: how do we select a suitable stepsize τ 0 ? The only limitations on τ 0 that we have discussed so far are related to the radius of convergence of the power series at t = t 0 and to the possible presence of spurious singularities of the Padé approximant on the positive real axis. As a matter of fact one needs to select time-steps that are far enough from this spurious pole of the Padé approximant to avoid catastrophic loss of precision.
To determine the stepsize in a suitable way, we consider the residual obtained by substituting the approximate solutionX(t) ≡ P [X(t)] [m,n] inside eq. (1). Since the residual vanishes at t = t 0 , R[X] t=t0 = 0, we expect it to grow in size at later times t > t 0 , unlessX is a solution of eq. (1) (in this case it vanishes identically). In other words, the residual grows monotonically in a sufficiently small region to the right of t 0 and we can constrain the loss of precision of the numerical solution for t 0 < t < t 0 + τ 0 , by choosing τ 0 in a suitable way. This is done by introducing a real parameter δ, with 0 < δ ≪ 1, such that The value of the time step dt 0 found in this way ensures that the error made in approximating the exact solution withX is bounded at all times t 0 ≤ t ≤ t 0 + dt 0 .
The process can now be iterated: moving at the time t 1 = t 0 + dt 0 we can repeat the previous steps by considering eq. (1) with initial conditions X(t 1 ) = X(t 1 ) and find a new time step and consequently new initial conditions at a time t 2 = t 1 + dt 1 .
At each stepsize there is the freedom of changing the order of the Padé approximant used in the calculation, so that we can regard our method as a variable stepsize variable order Padé-Taylor (VSVO-PT) method.
In the next section we illustrate the method by applying it to the van der Pol equation.

The van der Pol oscillator
The van der Pol equation [9] x is an example of nonlinear ordinary differential equation that supports a limit cycle: for µ > 0, regardless of the initial conditions (with the exclusion of x(0) =ẋ(0) = 0), the system evolves in time towards a limit cycle with fixed period and amplitude (both functions of the parameter µ). The trajectory that the solution describes in phase space spirals in or out (depending on the initial conditions) towards a close orbit that for 0 < µ ≪ 1 closely resembles that of a harmonic oscillator while for µ ≫ 1 results in a highly deformed loop (reflecting the presence of two very different -slow and fast -regimes). The van der Pol equation constitutes the perfect testing ground for our method, because it represents a challenging problem, which has attracted (and continues to attract) much interest in the literature. This interest is reflected in Fig. 1 where we report the number of publications per year found in google scholar under the keyword "van der Pol equation". The recent paper by Ginoux and Letellier provides an historical account on the origins of the van der Pol equation [10].
Following [11] we classify the strategies for studying the limit cycle of the vdP equation into three categories: • numerical studies [12,13,14,15,16,17], where the eq. (6) is solved numerically for some finite value of µ (apparently the largest value of µ for which the period and amplitude have been calculated is µ = 200 in ref. [14]); the accuracy of the results progressively worsens as µ gets larger 3 ; • asymptotic studies [23,24,25,26,27,28,29,14,30,31,32], where the behavior of the solution (particularly the period and the amplitude of the limit cycle) is established for µ → ∞; it is remarkable that as early as in 1944 Haag [25] and then in 1947 Dorodnitsyn [26] were able to approximate the asymptotic behavior of the period and the amplitude of the limit cycle of the van der Pol oscillator. For a discussion of the work of Dorodnitsyn, see ref. [33], written to celebrate the centenary of the birth of A.A. Dorodnitsyn.
Later works on the van der Pol oscillator are mainly interested in the study of the perturbation series for the solution, which is most efficiently obtained by applying the Lindstedt-Poincaré method. The numerical results of refs. [12,13,15,16,17] are used in those works to validate the accuracy of the resummation schemes adopted (see for instance Table 4 of [38], Table 2 of [39] and Table  2 of [40]). Strictly speaking, the use of the perturbation series for the period (amplitude) of the oscillator is limited to the region µ ≤ µ T (µ ≤ µ A ), where µ T is the radius of convergence of the series 4 . For µ > µ T,A , different techniques can be applied to improve the convergence of the series, with some success (see [38,40,11]). In this regime however the numerical results are far more precise than the results based on perturbation theory, which can explain only in part the lack of interest in deriving far more precise numerical results. In our view another reason for the scarcity of numerical results is the difficulty in obtaining accurate numerical solutions of eq. (6) in the stiff regime (µ ≫ 1).
We will now apply to the van der Pol equation the procedures explained in the previous section for the general case.
We start at t = t 0 , where the initial conditions x(t 0 ) = x 0 andẋ(t 0 ) = x ′ 0 are provided and consider the power series solution with c 0 = x(t 0 ) and c 1 =ẋ(t 0 ). The coefficients c k with k ≥ 2 can be obtained in terms of the nonlinear recurrence relation This equation allows one to calculate exactly a large number of coefficients of the power series and thus estimate its radius of convergence. In Fig. 2 we estimate the radius of convergence of the power series (7) for µ = 100, x 0 = 2 and x ′ 0 = 0, which turns out to be ≈ 1/26. In Fig. 3 we plot the real poles of the Padé approximant of order [n, n] for µ = 100, x 0 = 2 and x ′ 0 = 0: we see that Padé of order 3,5,7 and 9 have a positive real root. For example for n = 3 the root is approximately τ ≈ 0.03091366759, which constitutes an upper limit to the stepsize. Similar behaviors are also observed for smaller values of µ, but in that case the real positive roots are much larger (thus allowing to span a much greater interval of time).
In Fig. 4 we plot the residual for µ = 100, x 0 = 2 and x ′ 0 = 0, again using only the first 13 coefficients: the residual obtained using a diagonal Padé approximant is about 10 6 smaller than the corresponding residual built on the Taylor polynomial with the same number of coefficients. To get an idea, in this case the region for which R < 10 −30 is about 3 times larger when the Padé approximants are used.
There is an important advantage in using our method (or the Taylor method) for the calculation of the period and amplitude of the limit cycle (or more in general of any periodic solution): although these methods work by performing discrete time steps, t 0 → t 1 → t 2 → . . . , the solution is also known at any arbitrary time between any two consecutive times steps, in terms of a Padé approximant (in our case) or of a Taylor polynomial (for the Taylor method). In our method, when moving from time t i to time t i+1 , we expect that the maximum error occurs at t = t i+1 , where the residual takes the maximum value δ. In this way we can estimate within the same accuracy both the location of a maximum or a minimum and the value of the solution at this point.
The situation is be completely different using the Runge-Kutta method: for with a strong numerical evidence that µ 2 T = µ 2 A ≈ 3.420187909357 . . . (see ref. [11] for the value of µ accurate to 100 digits).  example, Odani [43] has estimated the amplitude of the van der Pol oscillator for a limited set of values of the coupling using the Runge-Kutta method with a stepsize of 2 −20 ! Even so the results reported by Odani contain only 5 decimals. Notice that refs. [43,44,45,46,47] deal with the problem of determining the maximum amplitude of the limit cycle.
In our case the procedure to estimate the period and the amplitude of the oscillator is straightforward: we pick the initial conditions (typically we choose x 0 = 2 and x ′ 0 = 0) and iterate our method until a change in the sign of the velocity is detected; at this point we know that a minimum of the solution can be found between the two last time steps and we can look for the zero of the derivative of the Padé approximant contained in this region. The value of the solution at this point (changed of sign) is then used as the new x 0 in the initial conditions (we set the time to zero at this point) and the whole procedure is repeated. The time elapsed from t = 0 to the instant where the derivative changes sign is essentially the half-period. This process is iterated until the amplitude and the period have converged to a sufficiently large number of digits (alternatively the same process can be carried out by looking for two changes of sign in the derivative, which in turn would provide directly the period).
Importantly, we have the nice bonus that the number of iterations that are needed to reach the desired level of accuracy is much smaller in the stiff regime, where the transient behavior of the solution decays more rapidly with time (the number of time steps however increases, both because of the larger period, which grows with µ, and because of the need for smaller steps in the regions where the derivative is large).
The accuracy that we can hope to reach for the period and the amplitude is essentially determined by δ (provided that the number of digits used in the calculation is sufficiently large): in this way, one can expect to improve the numerical estimates by using a sufficiently small δ and working with large enough digits. In the next section we provide an independent check of this claim can be done by comparing the purely numerical results with the analytical results both for µ ≪ 1 (perturbative) and for µ ≫ 1 (asymptotic -stiff).

Numerical results
We have obtained numerical results for the period and the amplitude of the limit cycle of the vdP oscillator for a series of values of µ, starting at µ = 1 and going up to µ = 500 (well beyond the range of values studied in the literature). In the tables 1 and 2 we report the period and the amplitude of the limit cycle, obtained by working with the Padé-Taylor method, for µ = 100 (stiff), with different values of δ. From these tables we see that the results appear to converge quite fast (we expect that the best results in the table, corresponding to δ = 10 −100 , have more than 100 digits correct).
Although Tables 1 and 2 provide an indication of the accuracy of our results, it is possible to perform an independent check by comparing these numerical results with the perturbative results obtained with the Lindstedt-Poincaré method (in ref. [11] we have worked out analytically the LP to order 308 and numerically to order 859 -in this case the coefficients are calculated with 1000 digits).
This comparison can be done directly only for µ < µ T,A ≈ 1.85 (radius of convergence of the perturbative series): for µ = 1 the sum of the first 308 terms of the LP series for the period and the amplitude are in excellent agreement with the numerical results obtained with the PTm: and The perturbative results, however, can be extended beyond the radius of convergence of the series, by using a Padé approximant (we opted for a diagonal Padé [150, 150]). For µ = 1, the agreement between the Padé approximant and the numerical results is now extended to more than 100 decimals: and Even for µ > µ T,A , where the LP series does not converge, there is a remarkable agreement between the numerical results and those obtained using the Padé approximants, for µ ≤ 10.
This comparison is performed in Fig. 5, where we appreciate that, for µ ≤ 4, the results agree to more than 20 decimals. This observation justifies using the Padé approximant for the amplitude to estimate quite precisely the maximum amplitude of limit cycle of the van der Pol equation, that has been studied in a series of papers, ref. [43,44,45,46,47], but never calculated with high precision. In Fig. 6 we plot the Padé approximant for the amplitude (solid line) and compare it with the upper bound conjectured by Odani [43]. Our estimated value for the maximum amplitude is A (max) = 2.0234222556061133094 for µ = 3.2940126635728034197 (which complies with Odani's conjectured upper bound 2.0235).
In Fig. 7 we compare the stepsize of the PTm (normalized to the period) for different values of µ; observe that the number of steps required to cover the whole period appears to scale as ∝ µ 5/3 (see Fig. 8). The corresponding solutions are displayed in Fig. 9.   The reader may be skeptical on the necessity of calculating the period, the amplitude and more in general the limit cycle with such large precision as we claim having reached, but we hope that our next application will convince the most critical reader.

Extrapolation
In this section we adapt the discussion of Ref. [48], where the Richardson extrapolation was applied to the finite difference calculation of the Laplacian eigenvalues on complicated non-tensor domains in two dimensions.
In the previous section we have obtained very precise numerical results for different quantities (period, amplitude, maximum velocity, Lyapunov exponent, etc) related to the limit cycle of the van der Pol equation for different values of µ. Let O be such quantity and consider a sequence of values O 1 , O 2 , . . . , corresponding to different values of the parameter, µ 1 , µ 2 , . . . .
For µ ≫ 1, the behavior of O(µ) may be approximated as where f 0 (µ) ≫ f 1 (µ) ≫ f 2 (µ) ≫ . . . , for µ → ∞. For the case O(µ) = T (µ), for example, several terms of the asymptotic formula has been derived in [23,24,25,26,27,28,29] and correspond to choosing If we dispose of N values of O, for different values of µ, and restrict our calculation to the first N functions f i (µ), we can write the system of equations which can be cast in matrix form as The solution to eq. (15) can be obtained either by using the inverse matrix or by using Cramer's rule.
A further improvement can be obtained by performing a Richardson-Padé extrapolation, as done in ref. [48]; in this case we assume where lim µ→∞ g k (µ) = 0 (a convenient choice is g k (µ) = 1 µ k , k = 1, 2, . . . ). By using a rational extrapolation, we may be able to describe more precisely also the behavior of the observable for smaller values of µ.
Using N + M + 1 numerical results one can write the system of equations which can be cast in matrix form as Once again we can solve these equations either by using the inverse ofR or by applying Cramer's rule.
Of course the implementation of the Richardson and Richardson-Padé extrapolation requires to identify as many f i (µ) as possible, which may be a difficult task in general. Some of these functions may be provided by previous asymptotic calculations, as for the case of the period and the amplitude of the limit cycle, but the remaining functions need to be "guessed". In an empirical approach one may select functions that lead to a faster convergence of the extrapolated results and discard functions whose coefficients are abnormally small.
Determining these functions rigorously is most likely a very difficult task: for the case of the Dirichlet eigenvalues of the Laplacian on a L-shaped membrane, for which the first correction to the lowest eigenvalue has been proved to behave as h 4/3 (h is the grid spacing), the next few non-integer exponents have not been derived rigorously and, quoting Kuttler and Sigillito [51], "the exact form of the first several terms in the asymptotic formula for specific regions where no boundary interpolation is required is a nice problem at about the level of a doctoral thesis". The fact that this problem is still open can be taken as an indication of its difficulty.
The extrapolation carried out in [48], where we tried to guess these exponents, resulted in obtaining very precisely the lowest eigenvalues of the L-shaped membrane (later confimed by the results obtained in ref. [52] using a different method) and the most precise estimate to date of the lowest eigenvalues of the pair of isospectral domains of ref. [53] (see also Refs. [49,50] for further example of applications of the Richardson and Richardson-Padé extrapolations).
In Tables 12, 13 and 14 we report the precise numerical values of the times T fast for different values of µ.
We perform a Richardson extrapolation using the functions and use the last 50 numerical values to estimate the leading behavior of T fast (µ) for µ ≫ 1.
Our analysis shows that where one can observe that 1.33333333333333333334 − 4/3 ≈ 6.6 × 10 −21 and 3.5305579106695527878 − 4 3 + log 9 ≈ 7.2 × 10 −17 . We then have The slow component of the period, T slow is obtained from the condition T slow = T − T fast . In this case the Richardson extrapolation provides T slow (µ) ≈ 1.6137056388801093811656 µ + 7.0143222313792851710947 The accuracy of the first 3 numerical coefficients allows to guess their exact values where α is the first zero of the Airy function Ai(x). We can improve the estimate of the last coefficient by performing the Richardson extrapolation on the data where the contributions of the three first terms is taken out: in this case we obtain the slightly better value −4.85382301173.
By combining the formulas for the fast and slow periods, we obtain the asymptotic formula for the period of the vdP oscillator in the form which coincides with eq. (4) of Ref. [14], apart from last term that reads − 1.3246 µ . The coefficient −1.3246 used by Ponzo and Wax had been estimated earlier by Urabe [13].
It is interesting to observe that the formulas for the period given by Dorodnitcyn [26,27] and Urabe [13] contain a different coefficient for the term log µ/µ, −22/9 and −1/3 respectively. Ponzo and Wax were able to establish what is now accepted as the correct coefficient −2/3 by performing a numerical calculation up to µ = 200 (see Table 1 of ref. [14]).
By performing a Richardson-Padé extrapolation we are able to describe the numerical values for the amplitude in an excellent way, up to µ ≈ 2 (see Fig. 10). In fact the Richardson-Padé extrapolation is so accurate that it reproduces the first 10 decimals of the value of the maximum amplitude obtained using the Padé approximant of order [150, 150] for the Lindstedt-Poincaré expansion! Similarly the value of µ corresponding to this maximum agrees within 8 decimals with the value previously obtained.
While the expression for the amplitude obtained with Richardson-Padé is too lengthy to be reported here, it is useful to approximate it for µ ≫ 1, where it reads: It is easy to guess the exact form for the first three coefficients and obtain which agrees with the asymptotic formula for the amplitude, while adding few extra terms.
In the case of the maximum ofẋ (corresponding to the numerical results in Tables 9, 10 and 11) using Richardson extrapolation we have found, on purely numerical grounds, the behavior This formula has not been obtained before. Finally we consider the Lyapunov exponent, for which the asymptotic behavior has been established in ref. [34]. By performing the Richardson-Padé extrapolation of the numerical results, we have found that, for µ ≫ 1, the  quantity −λ 2 /µ behaves approximately as The coefficient of the term µ 0 agrees to 24 digits with the predicted coefficient of [34], 3+log (16) 6−log (16) ; the coefficient of the term µ −4/3 agrees to 15 digits with the predicted coefficient of [34], − 3(log(256)−3)α 2(log(4)−3) 2 . Alternatively one can obtain a perturbative expression for −λ 2 /µ using the Lindstedt-Poincaré method; using the results previous obtained in [11] we can write the first few terms read 6 +O µ 12 (33) Notice that our coefficient of order µ 4 does not agree with the corresponding coefficient of eq.(14) of [34]. In Fig. 11 we plot the numerical results for −λ 2 /µ and compare them with the Padé approximant of order [50,50] built on the perturbative LP series (dashed line) and the Richardson-Padé extrapolation of the numerical results (solid line). In Fig. 12 we plot the quantity (−λ 2 /µ) P ade LP − (−λ 2 /µ) numerical . We observe a remarkable agreement between the perturbative Padé approximant and the numerical results for µ ≤ 5. The numerical results presented in [34] appear to have limited accuracy (even for µ = 1, the result reported in table I of [34], − λ2 µ ≈ 1.0648, appears to have only two digits correct).

Conclusions
The sentence that stiff ODE cannot be solved by means of explicit methods, such as the Taylor method, is sometimes used as a definition of stiffness (pag.2 of ref. [18]). In particular, refs. [4] and [5], that describe explicit Taylor methods, also state this limitation and consider only non-stiff applications.
In this paper we have devised a different explicit method, that uses Padé approximants rather than Taylor series; our method does not suffer from this limitation and in fact it can be used to obtain arbitrarily accurate numerical solutions both for stiff and non-stiff problems. The accuracy of the solutions is limited only by the computing time and by the computer memory, provided that one works with a sufficiently large number of digits.
To illustrate the virtues and potentialities of our method we have applied it to the van der Pol equation, that is often used as an example of stiff ODE. The advantage of considering this equation lies on the fact that there it has been widely studied, both numerically and analytically: in particular, the perturbative Lindstedt-Poincaré expansion has been calculated to very high order in our previous paper [11], while the asymptotic behavior for large coupling has been established with the pioneer work of Dorodnitcyn [26,27] and subsequent improvements by Urabe [13] and by Ponzo and Wax [14]. All the numerical calculations of both period and amplitude of the vdP oscillator are quite old and with a limited precision. Despite the fact that the vdP equation is a standard example in many textbooks, we are unaware of recent attempts in improving these numerical results. With the exception of ref. [14], where the numerical experiment has lead to the determination of the correct coefficient of the term 1/µ 2 in the expression for the period, the numerical approximations have only been used to provide a qualitative validation of the asymptotic behaviors of the period and amplitude.
By disposing with our method of extremely precise numerical results (roughly 100 digits) we have been able to infer the asymptotic behaviors of the period and of the amplitude of the vdp oscillator by means of both Richardson and Richardson-Padé extrapolations. The expressions for the period and the amplitude reproduce the asymptotic formulas of [26,27,13,14]. Similarly our extrapolation of the Lyapunov exponent agrees with the formula of ref. [34]. Additionally, we have been able to infer the asymptotic behavior of the portion of the period corresponding to fast (slow) motion, and of the maximum amplitude of the velocity, which have never been derived before. It would be nice that our numerical work could stimulate the interest in deriving such behaviors from an asymptotic analysis.
We plan to apply our method to discuss further examples in future work.