{'fraction_non_alphanumeric': 0.05013943760167325, 'fraction_numerical': 0.013769463165233559, 'mean_word_length': 4.424834541443429, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Volterra-series approach to stochastic nonlinear dynamics: linear response of the Van der Pol oscillator driven by white noise

The Van der Pol equation is a paradigmatic model of relaxation oscillations. This remarkable nonlinear phenomenon of self-sustained oscillatory motion underlies important rhythmic processes in nature and electrical engineering. Relaxation oscillations in a real system are usually coupled to environmental noise, which further complicates their dynamics. Determination of the equation parameter values becomes then a difficult task. In a companion paper we have proposed an analytical approach to a similar problem for another classical nonlinear model--the bistable Duffing oscillator. Here we extend our techniques to the case of the Van der Pol equation driven by white noise. We analyze the statistics of solutions and propose a method to estimate parameter values from the oscillator's time series.


I. INTRODUCTION
Balthazar van der Pol introduced the concept of relaxation oscillations together with his eponymous equation for a simplified dynamics of a triode electric circuit [1,2]. Regarded as a power-series approximation for a more general class of Lienard systems [3,Sec. 7.4 and 7.5], this model became a paradigm of self-sustained oscillatory motion [2]. Besides its applications in engineering, the Van der Pol equation and its generalizations are used to describe various rhythmic processes in biology [4][5][6][7][8][9][10][11][12].
A general form of the Van der Pol equation that we consider in this paper extends the harmonic oscillator by introducing a nonlinear dissipative term in the equation for an unknown function of time x(t) and an external force f (t); the constants a, b, and c are, respectively, the friction coefficient, the stiffness, and the Van der Pol damping parameter. Because the above equation is of second order in time, the phase of this system is specified by two degrees of freedom (x,ẋ). The self-sustained oscillations of the Van der Pol equation correspond to its limit-cycle solution. In absence of the external force f (t), all trajectories of this system relax to a periodic orbit in phase space. Self-sustained oscillations exist if the friction constant in Eq. (1) is negative (a < 0). The Van der Pol system is stable when the parameters b and c are both positive. The amplitude of the limit cycle, which encircles an unstable equilibrium point in phase space, shrinks to zero when a = 0 and * belousov.roman@gmail.com disappears for a > 0. Therefore the Van der Pol oscillator with a ≥ 0 behaves as a monostable system. This dynamical regime is not studied in the present paper and should be treated by a different approach [13,Appendix A].
Environmental noise, which intertwines with relaxation oscillations in real systems, is often modeled by a stochastic force f (t) = Aẇ(t), with a constant amplitude A > 0 and a Gaussian white noiseẇ(t) of zero mean and unit intensity. One must usually resort to complex measures to determine the model's parameter values for this class of stochastic nonlinear problems [9][10][11][12].
In the companion paper [13] we demonstrated that time series of a second-order dynamical system-the stochastic Duffing oscillator-contain enough information to infer the parameter values of the underlying nonlinear model. Here we extend our analysis to the case of Van der Pol relaxation oscillations driven by white noise. We derive analytical expressions for approximate solutions and time-series statistics of Eq. (1). These formulas are then used to devise a parametric method of inference.
Our approach is based on the functional series of Volterra [14,15], which we expand up to the linearresponse term. The analytical results and the inference method that we propose are therefore applicable to relatively small noise amplitudes A; more details on the system's physical scales are given in Sec. III. Because we must approximate the limit-cycle solution of Eq. (1), for which no closed-form expression is known, our development is restricted to moderately nonlinear regimes.
Even in the absence of external driving, the statistical properties of the relaxation oscillations are far from trivial (Appendix A). This feature of the Van der Pol equation renders the time-series analysis more difficult than in the case of the Duffing oscillator [13]. Our inference method can therefore estimate the noise amplitude A only to the order of its magnitude (Sec IV).
Volterra series provide a representation of solutions for nonlinear stochastic equations. Other theoretical approaches, such as the Fokker-Planck equation and path integrals, focus instead on statistical properties of an ensemble of systems' realizations and offer less information about their dynamics. The theoretical tools may complement each other, for instance, the Volterra series may be used to an advantage in ergodic problems when time averaging is more convenient than ensemble averaging for the evaluation of statistical properties.

II. THEORY
A. Linear response of the Van der Pol oscillator A Volterra series is a polynomial functional expansion of the form in which g 1 (t) and g 2 (t) are the Volterra kernels of the linear and quadratic terms in the argument function f (t). Provided that the above series exists, a truncated expansion (2) approximates solutions of Eq. (1) driven by a small external force: (3) in which we neglect terms of the second and higher orders in f (t). The functions x 0 (t) and γ 1 (t) can be found by using the variational approach [13,15,Sec. 3.4], which yields a set of equations · · · Equation (4), which uniquely defines x 0 (t) for a given initial condition x(0),ẋ(0) , is equivalent to the homogeneous Van der Pol problem-Eq. (1) with f ≡ 0. The linear Eq. (5), which determines the first-order Volterra term γ 1 (t), in general contains time-dependent coefficients.
Because the Volterra series generalizes the Taylor-Maclaurin expansion of functions in calculus [15, Sec. 1.5], Eq. (2) may be restricted by a radius of convergence or may even fail to exist for some choices of x 0 (t). The equilibrium point x 0 (t) ≡ 0, which is a convenient choice for the monostable case of Eq. (1), is unstable in the regime of relaxation oscillations and yields a divergent kernel g 1 (t). With x 0 (t) ≡ 0 we therefore cannot construct an approximate representation (3) that is valid for long time scales [13].
In the following the homogeneous term x 0 (t) represents the stable limit cycle of Eq. (4). Because a closed-form expression for this solution is unknown, as its approximation one may adopt a single-mode Fourier expansion obtained by the method of harmonic balance [16,Sec. 4.4]: in which α = 2 −a/c. Equation (6), substituted into (5) for x 0 (t), yields in which are the time-dependent periodic coefficients. To continue the chain of simplifications, we replace the periodic coefficients in Eq. (7) by their mean values in which we relate the ensemble average of a periodic function to its time average over one period 2π/ √ b. In this approximation Eq. (7) describes a harmonic oscillator γ ξ ≈ γ 1 :γ with the linear response function The linear-response term γ ξ (t) in the above equation has a Gaussian probability density, with a zero mean γ ξ = 0 and an autocovariance function [17] Even without external driving f (t), the relaxation oscillations of Eq. (1) have nontrivial statistics that are accurately described by the approximate solution ξ(t) (Appendix A). Because the covariance ξ(t)γ ξ (t ) = 0 of the two independent contributions ξ(t) and γ ξ (t) vanishes, one easily obtains from Eqs. (6), (14), and (15) the autocorrelation function: The approximate autocorrelation function χ ξ (t) is one of our main results.

B. Parametric inference
The approximate solution x ξ (t) is limited to small noise amplitudes not only due to the truncation error in Eq. (2). When the external force is strong enough to drive the system's trajectory through or very close to the special point (x,ẋ) = (0, 0), the phase of the periodic term ξ(t) can be shifted by an angle φ: ξ(t) → ξ(t + φ).
As in the case of the stochastic Duffing oscillator [13], we circumvent the above issue by applying Eq. (14) locally: the trajectory of x(t) can be split into pieces x + (t) and x − (t) for respectively x(t) > 0 and x(t) < 0. The duration of each component corresponds approximately to a half period π/ √ b of ξ(t). The above phase shifts which are not taken into account by Eq. (14), also decorrelate the time series of x(t). The persistent periodic term ∝ cos 2 ( √ bt) in Eq. (16) is therefore accurate only at small time scales. Unlike χ ξ (t), the actual autocorrelations of x(t) decay to zero for t → ∞ unless f (t) ≡ 0. For small amplitudes A the function χ ξ (t) approaches its first zero as t → τ ≈ π/(2 √ b). We may then apply the criterion of Lagarkov and Sergeev [17,19] to select the interval 0 ≤ t ≤ τ over which Eqs. (14) and (16) are expected to be accurate.
Although Eq. (16) predicts accurately the initial decay of χ(t) (Sec. IV), fitting this expression to autocorrelations of the time series x(t) overestimates the absolute magnitude of the parameter a. As noted in the companion paper [13], the value of a is sensitive to errors of the approximations used. In Appendix C we propose therefore a correction formula Equations (14) and (17) are sufficient to implement a simple technique of parametric inference. First one can fit empirical autocorrelations of x(t) byχ(t) over the time interval 0 ≤ t ≤ τ and thus determine the parameter values a and b. Considering γ ξ (t) in Eq. (14) as a Gaussian error term, one may then fit the trajectory pieces x ± (t) by in which the unknown parameter values account for the relative phase shift √ bt 0 of each trajectory piece. One then estimates The numerical error of fitting the approximate expression ξ(t) to trajectory pieces x ± (t) eventually exceeds the uncertainty of the driving noise f (t) of a small ampltitude A → 0. One may even find a situation, in which the parameter A cannot be evaluated reliably through Eq. (20). For this reason we propose to use as a proxy the following formula in which is the standard deviation of α over the sample of trajectories x ± (t). This expression is motivated by the fact that, when the homogeneous term ξ(t) ∝ α dominates the statistical variability of the data, the magnitude of A → 0 cannot be inferred with a precision superior to √ā b . We also adopt ∆A = √ā b as the error estimate for A.
in which the parameter µ = −a/ √ b controls the nonlinear character of the dynamics. The greater its value, the larger is the amplitude of the relaxation oscillations.
Two control parameters of Eq. (22) that are not fixed in the system of reduced units are µ and A. Without external driving the Van der Pol oscillator, which orbits around the origin of the phase space with the amplitude α = 2 −a/c in the harmonic potential U (x) = bx 2 /2 [Eqs. (1) and (6)], has an energy scale U (α) = bα 2 /2 = 2µ b 3/2 /c. For increasingly large values of µ Eq. (6) becomes progressively less accurate. As a typical value we choose therefore µ = 1. The techniques that we propose in the present paper should work also for µ > 1, but their precision deteriorates for larger values of this parameter. The energy scale of the external force f (t) = Aẇ(t) is A 2 / √ b. One might therefore expect the small-force expansion Eq.

IV. NUMERICAL RESULTS
To test the theoretical expressions of Sec. II A, we simulated Eq. (1) for selected values of the parameters µ and A (Sec. III). The computational details are summarized in Appendix B. Below we report the results of our simulations in the reduced units.
A comparison of the noisy oscillations of x(t) with the theoretical expression for ξ(t) is shown in Fig. 1. As a representative example we use the parameter values µ = 1 and A = 0.6 and the initial conditions (x,ẋ) = (α, 0). The trajectory ξ(t) approximates well the period of oscillations and the overall trend of x(t).
For the same parameters as above, we compare the autocorrelation function χ(t) of the simulated time series x(t) and the theoretical curve given by Eq. (16) in Fig. 2. The analytical expression of χ ξ (t) predicts well the initial decay of the autocorrelation. However, a moderate phase difference accumulates at longer times. Moreover, as explained in Sec. II B, χ(t) decays to zero with time to zero unlike the approximate Eq. (16). An amplitude difference between χ ξ (t) and the actual autocorrelations therefore becomes evident in Fig. 2 at longer times.
The efficiency of the parametric-inference method we proposed in Sec. II B is demonstrated in Table I. The best results are obtained for a moderate level of noise 0.5 A 1.0. The parameters a and b are inferred most reliably. Because relaxation oscillations have nontrivial statistics even in the absence of external forces, it is difficult to discriminate between the numerical errors of fitting approximate expressions and the uncertainty due to the driving noise. Our method estimates the ampli-    tude A of the stochastic force f (t) only within its order of magnitude.

V. CONCLUSION
Using the Volterra series we have analyzed statistical features of a noisy Van der Pol equation. Perhaps surprisingly, its solution can be decomposed within the linear order of the driving force into two independent contributions: a deterministic part that describes relaxation oscillations, and a stochastic linear-response term. With the help of simple approximation schemes we showed that the deterministic contribution has a singular probability density, whereas the stochastic part can be described by a Gaussion process with a second-order autocorrelation function.
The inference method based on our analytical results allows us to estimate parameter values of the stochastic Van der Pol model from observed time series of oscillations for moderate levels of the driving noise. For small noise amplitudes only an order of magnitude can be determined due to insufficient accuracy of our simple approximations. More accurate approximation approaches exist, for example, the two-timing method [3, Sec. 7.6], and could be pursued to improve the precision of the the analytical formulas.
The mean and the autocovariance function of ξ(t) can be evaluated by time averaging The probability density of ξ has two singularities at the ends of its support interval ξ = ±α. Histograms of the time series ξ(t) and x(t) have two distribution modes near these points. In the companion paper [13] we have succeeded in fitting a bimodal probability density of the noisy Duffing oscillator to an approximate expression that was derived from a power series for the exponential family of random variables. This approach unfortunately fails however in the case of the noisy Van der Pol oscillator: such an expansion may not exist near the two singularities at which p(x 0 ) tends to infinity.
For the probability density of x ξ (t) [Eq. (14)], regarded as the sum of two independent variables ξ(t) and Gaussian variable γ ξ (t), there is no simple analytical expression. However the Fourier image η(X ξ )-the characteristic function of x ξ for the reciprocal variable X ξ -can be obtained in a closed form. Because γ ξ (t) is Gaussian, we have in which is the characteristic function of ξ and J 0 (·) is the zero thorder Bessel function of the first kind. The autocovariance of x ξ is simply the sum of Eqs. (A3) and (16). In Fig. 3 we compare the empirical characteristic function of x(t) (Appendix B) with Eq. (A4) for two representative examples. Our analytical expression for η(X ξ ) is accurate at least for X ξ x −1/2 even for large noise amplitudes A. The theory is in excellent agreement with the simulations for A ≡ 0.

Appendix B: Simulation algorithm
In our computational experiments we use a companion system of Eq. (1) with X = (x, y) = (x,ẋ): We adopt a second-order operator-splitting approach [21] for stochastic systems [13,22,Appendix C], by decomposing the time-evolution operator T aṡ in which The can be approximated by The action of individual operators of the form exp(L∆t) is then determined by the differential equatioṅ The operators T x , T y , and T f produce linear equations of the above type. The composite operator (B3) leads then to the following integration algorithm for Eq. (B2) y(t + ∆t) = y(t) exp{−[a + cx(t + ∆t/2) 2 ]∆t} − exp{−[a + cx(t + ∆t/2) 2 ]∆t/2}bx(t + ∆t/2)∆t x(t + ∆t) = x(t + ∆t/2) + y(t + ∆t)∆t/2.
For each value of the parameter A, the simulations reported in Sec. IV involved 5 × 10 5 time steps of a size ∆t = 0.01 in the reduced units. Statistics were calculated from single trajectories sampled at time intervals of 0.05, which included 10 5 observations.
Appendix C: Correction of the linear-response formula Equation (17) improves on Eq. (16) by replacing the effective friction parameter a ξ = −a withā [cf. Eqs. (7)- (13)]. To obtain this correction, we first refine the approximation of the homogeneous term x 0 (t) by use of a variational Green-function approach.

Simulation
Eq. This iterative method resembles that of Refs. [23,24], but is based on a different Piccard operator.