Hypergeometric First Integrals of the Duffing and van der Pol Oscillators

TL;DR

This paper demonstrates the existence of hypergeometric first integrals for the Duffing and van der Pol oscillators, revealing new non-analytic integrals expressed via hypergeometric functions and providing a general criterion for such phenomena.

Contribution

It introduces hypergeometric first integrals for these oscillators and formulates a criterion for their existence in general dynamical systems.

Findings

Existence of Liouvillian first integrals expressed through hypergeometric functions.

These integrals are neither analytic nor algebraic.

A general criterion for hypergeometric integrals in dynamical systems.

Abstract

The autonomous Duffing oscillator, and its van der Pol modification, are known to admit time-dependent first integrals for specific values of parameters. This corresponds to the existence of Darboux polynomials, and in fact more can be shown: that there exist Liouvillian first integrals which do not depend on time. They can be expressed in terms of the Gauss and Kummer hypergeometric functions, and are neither analytic, algebraic nor meromorphic. A criterion for this to happen in a general dynamical system is formulated as well.

Figures

Equations

$$\ddot{u}+\dot{u}+{\omega }_0^{2}u+u^3=0,$$

$$\ddot{u}+\alpha \dot{u}+{\omega }_0^{2}u+\gamma u^n=0.$$

$$W=\frac{3}{\sqrt{2}}{\mathrm{e}}^{t/3}u;Z=-\sqrt{2}{\mathrm{e}}^{-t/3},$$

$$W^{\prime \prime }(Z)=-W^3.$$

$$H=\frac{1}{2}{W}^{\prime 2}+\frac{1}{4}{W}^4,$$

$$I_{\mathrm{D}}={\mathrm{e}}^{\frac{4}{3}t}\left(\frac{1}{2}{\dot{u}}^2+\frac{1}{3}u\dot{u}+\frac{1}{18}{u}^2+\frac{1}{4}{u}^4\right).$$

$${\omega }_0^2=\frac{2(n+1)}{(n+3{)}^2},$$

$$I_1=V^{-\frac{1}{n+1}}\left[V^{\frac{1}{2}}+(n-1){\omega }_0\varsigma u{}_2{F}_1\left(\frac{1}{2},\frac{1}{n+1},\frac{n+2}{n+1},\frac{4u^{n+1}}{V}\right)\right],$$

$$I_2=V^{\frac{1}{2}-\frac{1}{n+1}}-\frac{n-1}{n+1}{2}^{\frac{n-1}{n+1}}{\omega }_0^2{\varsigma }^{\prime }\frac{2u+(n+3)\dot{u}}{\sqrt{V}}{}_2{F}_1\left(\frac{1}{2},\frac{n}{n+1},\frac{3}{2},1-\frac{4u^{n+1}}{V}\right),$$

$$\ddot{u}+(1+\beta u^m)\dot{u}+{\omega }_0^{2}u+u^n=0,$$

$$I_{\mathrm{P}}={\mathrm{e}}^{\frac{n}{\beta }t}\left(\dot{u}+\frac{\beta -n}{\beta }u+\frac{\beta }{n}{u}^n\right).$$

$${\omega }_0^2=\frac{n}{(n+1{)}^2},\beta =n+1\text{and}m=n-1,$$

$$I_3={\left|u^n+{\omega }_0^2(u+\beta \dot{u})\right|}^{1-\frac{1}{n}}\left[1+\frac{(n-1)u}{u+\beta \dot{u}}{}_2{F}_1\left(1,1,1+\frac{1}{n},-\frac{u^n/{\omega }_0^2}{u+\beta \dot{u}}\right)\right],$$

$$I_4={\left|u^n+{\omega }_0^2(u+\beta \dot{u})\right|}^{1-\frac{1}{n}}\left[1-\frac{u}{u+\beta \dot{u}}{}_2{F}_1\left(1,1,2-\frac{1}{n},1+\frac{u^n/{\omega }_0^2}{u+\beta \dot{u}}\right)\right],$$

$$\ddot{u}+(1+\beta u^{n-1})\dot{u}+{\omega }_0^{2}u+u^n+\phi u^{2n-1}=0,$$

$${\omega }_0^2=\frac{n}{(n+1{)}^2},\beta =n+1,\text{and}\phi =\frac{1}{4{\omega }_0^2},$$

$$I_5=|V{|}^{-\frac{1}{n}}\left[\exp \left(\frac{2(n-1)\phi u^n}{nV}\right)V+(n-1)u{}_1{F}_1\left(\frac{1}{n},1+\frac{1}{n},\frac{2(n-1)\phi u^n}{nV}\right)\right],$$

$$\dot{x}$$

$$\dot{y}$$

$$DF=KF,\text{for some }K\in \mathbb{C}[x_1,\dots ,x_d],$$

$$F_1$$

$$F_2$$

$$F_{3,4}=2\pm 3\sqrt{2}\mathrm{i}x+6y,K_{3,4}=-\frac{2}{3}\pm \sqrt{2}\mathrm{i}x-y,$$

$$\frac{D\left(x^2(9x^2+2(1+3y{)}^2)\right)}{x^2(9x^2+2(1+3y{)}^2)}=-\frac{4}{3},$$

$$J:=x^2(9x^2+2(1+3y{)}^2)={\mathrm{e}}^{-4t/3}{J}_0,$$

$$D(F_1/F_2)=-({\partial }_{x}P+{\partial }_{y}Q)(F_1/F_2).$$

$${\alpha }_1{k}_1+{\alpha }_2{k}_2+{\alpha }_3{k}_3={\alpha }_0\in \mathbb{C},{\alpha }_i\in \mathbb{C}\text{not all zero},$$

$$k_1-k_2=f_3,$$

$$\int J\mathrm{d}t$$

$$\frac{J}{{\alpha }_0}+C_1$$

$$I=\frac{J}{{\alpha }_0}-\frac{{\zeta }^{{\alpha }_1}}{{\alpha }_1}{}_2{F}_1(1-{\alpha }_3,{\alpha }_1,{\alpha }_1+1,\zeta ),$$

Full text

1

The Duffing and the van der Pol oscillators are among the simplest (at least in form) dynamical systems, present in biology [1, 2], nano-electronics and mechanics [3] as well as many subfields of physics [4, 5] (see the last one for a more comprehensive list of references). The systems’ integrability, or solvability, is still an active topic with only partial results available [6, 7, 8] and the aim of the present article is to use them as examples of how non-analytic Liouvillian first integrals can arise in polynomial differential equations. This is the content of Lemma 1 in Section 2 and its generalization in Section 5.

The force-free form of the basic Duffing oscillator is

[TABLE]

or, more generally,

[TABLE]

With a suitable rescaling of $u$ and the time, $\alpha$ and $\gamma$ can be made equal to 1, so the only essential parameters will be $\omega_{0}$ and $n$. Although the case of immediate interest is that of integer $n>1$, the extension to real $n$ is straightforward if $u^{n}$ is replaced by $|u|^{n-1}u$.

For the specific value of $\omega_{0}^{2}=\tfrac{2}{9}$, explicit solutions of (1) in terms of Jacobian elliptic functions were found in [6]. The authors note that the equation admits a transformation

[TABLE]

which turns the equation into a solvable one

[TABLE]

It can be regarded as Hamiltonian with

[TABLE]

which is then conserved, but upon going back to $u$ and $\dot{u}$, the first integral defined by $H$ becomes time dependent:

[TABLE]

A generalization of the above approach, which utilizes transformations of the form (3), was formulated in [7], and gives a more constructive method of finding first integrals.

As the solutions are explicitly available, this time dependence of $I_{\mathrm{D}}$ does not seem to be a practical problem – one could argue that a first integral can be obtained (at least locally) by elimination of time between the solutions, but that is seldom practically computable. Thus, the question emerges of whether a “proper” time-independent first integral exists. The goal of this investigation is to show that this is indeed the case for the general unforced Duffing oscillator, and partly so for the Duffing-van der Pol oscillator. What is more, the existence of such an integral can be ascertained without the a priori knowledge of the solutions.

The first main result of the present letter is the following

Theorem 1**.**

If the frequency $\omega_{0}$ of the Duffing oscillator (2) with general nonlinear term of degree $n$ satisfies

[TABLE]

then the equation has a Liouvillian first integral. When $n$ is odd, in the half-planes separated by $2u+(3+n)\dot{u}=0$ the integral is given by

[TABLE]

where $V=4u^{n+1}+(2u+\omega_{0}^{2}(n+3)\dot{u})^{2}$ and $\varsigma=\mathrm{sign}(2u+(n+3)\dot{u})$; alternatively, in the half-planes separated by $u=0$ it is

[TABLE]

with $\varsigma^{\prime}=\mathrm{sign}(u)$ and $V$ as above. Each ${}_{2}F_{1}$ is real when its argument is less than 1.

When $n$ is even, so that $V=0$ defines a real invariant curve $\mathcal{C}$, the above expressions remain valid if $I_{1}$ is taken to the right of $\mathcal{C}$, and $I_{2}$ to its left with $\varsigma^{\prime}=(-1)^{n/(n+1)}$.

The second equation of interest is the Duffing-van der Pol oscillator, which exhibits the same type of time-dependent first integrals. The equation is

[TABLE]

and it was shown in [8] to admit a whole parametric family of time-dependent integrals. The required condition is that $m=n-1$ and $n^{2}+\beta^{2}\omega_{0}^{2}=n\beta$, and the conserved quantity is

[TABLE]

Similarly to the previous system, this quantity can be used to obtain a time independent first integral, leading to the second result:

Theorem 2**.**

If the parameters of the Duffing-van der Pol oscillator (10) with general nonlinear term of degree $n$ satisfy

[TABLE]

then the equation has a Liouvillian first integral given by

[TABLE]

which is real when $-\omega_{0}^{-2}u^{n}/(u+\beta\dot{u})<1$; or alternatively by

[TABLE]

which is real when $\omega_{0}^{-2}u^{n}/(u+\beta\dot{u})<0$. Additionally, the singularities of the hypergeometric function correspond to the invariant sets $\omega_{0}^{2}(u+\beta\dot{u})+u^{n}=0$ and $u+\beta\dot{u}=0$.

The above system can be further generalized to

[TABLE]

which was analyzed with the Prelle-Singer procedure in [9], and with the Lie symmetry method in [10]. It differs from the first two in that the specific structure of its invariant sets makes the hypergeometric first integral confluent:

Theorem 3**.**

If the parameters of the generalized Duffing-van der Pol oscillator (15) satisfy

[TABLE]

then the equation has a Liouvillian first integral given by

[TABLE]

where $V=u+\beta\dot{u}+2\varphi u^{n}$. The above integral is real, and its only singularity corresponds to the invariant curve $V=0$.

The appearance of hypergeometric functions is quite remarkable – it shows how restrictive the notion of analytic, meromorphic or algebraic integrability can be. Although various parametric forms of solutions of polynomial dynamical systems can be given in terms of hypergeometric functions [11], this seems to be the first example of a natural (i.e. not artificially contrived) system with first integrals of such form. This is the price we have to pay for the transition from relatively simple, but time-dependent, integrals to time-independent ones, which permit us to call an autonomous system integrable in the standard sense.

A starting point in integrability analysis might be looking in the class of analytic functions, but because trajectories converge to the node at the origin, this approach will not work here. As the systems are of dimension two, it is not possible to directly apply the basic methods of differential Galois theory [12] either, because the first normal variational equation will be of order 1 and thus soluble. Higher variational equations could be used, but in addition to complicating the analysis, it will still be limited to meromorphic first integrals. To make progress, we will have to turn to the Darboux polynomials, which are most often used in the context of rational first integrals, but can lead to algebraic or even Liouvillian expressions.

The remainder of the article is devoted to introducing the basic machinery (Section 2.) and providing a general criterion for the appearance of hypergeometric first integrals in polynomial systems (Lemma 1). The proofs of the theorems are given in Sections 2 through 4, and the central criterion is further extended in Lemma 2. The possibility of reduction to more elementary functions is discussed in the final section.

2

The first step in the construction of first integrals is to consider the dynamical system corresponding to (1)

[TABLE]

where $x=u$ and $y=\dot{u}/u$. The vector field associated to the flow defines a derivation over the polynomial ring $\mathbb{C}[x,y]$, and one can look for the Darboux polynomials of the derivation $D:=P\partial_{x}+Q\partial_{y}$. If enough of them are found, the second step is to combine them into a first integral.

To briefly review the general setting, consider a $d$-dimensional system or, equivalently, a derivation over complex polynomials of $d$ variables. A Darboux polynomial $F$ is defined as an elements of $\mathbb{C}[x_{1},\ldots,x_{d}]$, such that

[TABLE]

and $K$ is called the cofactor. Since along a solution one has $DF=\dot{F}$, it follows from the definition that $F=0$ is an invariant set, and if $K\equiv 0$ then $F$ is just a polynomial first integral. Another basic property holds for the product: $D(F_{1}F_{2})=(K_{1}+K_{2})F_{1}F_{2}$, so that $F_{1}F_{2}$ is itself a Darboux polynomial. In particular, any integral power $F^{m}$ is a Darboux polynomial with the cofactor $mK$. A converse result can also be proven: if $F$ is reducible, its factors must be Darboux polynomials too.

Importantly, the existence of such polynomials allows one to analyze integrability through rational functions, because if $f/g\in\mathbb{C}(x_{1},\ldots,x_{d})$ is a first integral, then necessarily $f$ and $g$ must be Darboux polynomials with the same cofactor; conversely, if there are enough Darboux polynomials, so that $K$’s are linearly dependent over $\mathbb{Z}$, a rational first integral exists. For the full Darboux theorem see [13].

Going beyond polynomials, $F^{q}$, with rational $q$, is an algebraic function which still has a polynomial cofactor. Indeed for any function $f$, Liouvillian over $\mathbb{C}(x_{1},\ldots,x_{d})$, there might exist a polynomial $K$ such that $Df=Kf$; for brevity the term Darboux pair will be used to refer to $\{f,K\}$ then.

The last general remark to be made is that the degree of $K$ has to be lower than that of the derivation, and finding variables for which $D$ is quadratic will simplify things considerably. However, neither this, nor $d=2$ is assumed for the general lemmas given below. For an exposition of the subject in the context of polynomial derivations see also [14].

Coming back to the present case of (18), the derivation is of degree 2, so the cofactors can be at most linear, and assuming some degree of $F$, the equation $DF=KF$ can be solved term by term. It is a straightforward calculation to verify that when $\omega_{0}^{2}=\tfrac{2}{9}$ the following quadratic Darboux polynomials exist

[TABLE]

As mentioned, the irreducible factors of $F_{2}$ are themselves Darboux polynomials

[TABLE]

but with a view to finding a real first integral, $F_{2}$ will be used.

It should also be noted, that there is no effective general tool to give the bounds on the degree of $F$, so when $\omega_{0}^{2}\neq\tfrac{2}{9}$ the question of existence of higher degree $F$ remains open. In general it is a difficult task to exclude all possible degrees, see for example [15], but fortunately the goal here is not to find all the polynomials or rational first integrals.

Thanks to the product property, Darboux polynomials (20) can be combined, and the new cofactor can be made constant

[TABLE]

so that we immediately have

[TABLE]

where $J_{0}=J(x(0),y(0))$ is equivalent to $I_{\mathrm{D}}$ of (6).

If $K$’s are dependent over $\mathbb{Z}$, the resulting cofactor can even be made zero, so a rational first integral is found (or at least algebraic, were they dependent over $\mathbb{R}$). This is not the case here, and there are no new Darboux polynomials of degree 3 either, that is, they are all products of $F_{1}$ and $F_{2}$. In principle, one could try looking at higher degrees, but as it turns out this will not be necessary.

An important thing to notice is also that the following combination of (20) is related to the divergence of the flow (18)

[TABLE]

This means that a Liouvillian first integral can be obtained, as described by Theorem 1 of [13] with $F_{1}/F_{2}$ being the integrating factor of the form $P\mathrm{d}y-Q\mathrm{d}x$. Instead of applying that theorem, however, a slightly different derivation will be given here, which leads directly to a more concise expression for the first integral. It relies on a result valid in a very general setting:

Lemma 1**.**

Let $D=\sum_{i=1}^{d}P_{i}\partial_{i}$ be a derivation which admits two Darboux polynomials $f_{1}$ and $f_{2}$ with cofactors $k_{1}$ and $k_{2}$, respectively. If $f_{3}:=f_{2}-f_{1}$ is also a Darboux polynomial with cofactor $k_{3}$, such that

[TABLE]

and $f_{3}$ is the cofactor of $f_{1}/f_{2}$, i.e.,

[TABLE]

then the dynamical system $\dot{x}=P_{i}(x)$ associated with the derivation $D$ has a Liouvillian first integral, expressible locally in terms of the hypergeometric function.

Proof. The first condition means that we have a function $J=\prod f_{i}^{\alpha_{i}}$ such that $DJ=\alpha_{0}J$ or $J(t)=\mathrm{e}^{\alpha_{0}t}J(0)$. Let us thus take the time integral of $J$ and re-express it with $f_{i}$ and a new variable $\zeta:=f_{1}/f_{2}$ to get

[TABLE]

where the determination of the complex argument of $\zeta$ or $1-\zeta$ was ignored, with the understanding that the form of the integral remains the same save for a multiplicative constant. This constant depends on the region in which $\zeta$ (and hence $x$) lies, which is why the word “locally” is necessary in the conclusion.

Because $\alpha_{i}$ are defined up to rescaling, we can take $\alpha_{1}+\alpha_{2}+\alpha_{3}=1$ and the last integral becomes the Euler representation of the hypergeometric function, so that

[TABLE]

is the sought first integral. Because the Euler representation is obtained from an indefinite integral, the result is a Liouvillian function of $\zeta$, which in turn is a rational function of the original variables. ∎

Thus, in the standard case of equation (1), the above lemma can be applied with

[TABLE]

but the resulting first integral will then explicitly contain the imaginary unit. A slight modification is required if one insists on real expressions, and it can be effected for the more general equation (2).

Proof of Theorem 1. Let us introduce a new set of variables $x=u^{n-1}$ and $y=\dot{u}/u$, in which system (2) reads

[TABLE]

As before, $F_{1}=x$ is a Darboux polynomial, and direct computation shows that

[TABLE]

is another one if the condition $\omega_{0}^{2}=2(n+1)/(n+3)^{2}$ is satisfied. As above, a function $J$ with constant logarithmic derivative can be found to be

[TABLE]

The new variable to use is

[TABLE]

for which we have

[TABLE]

and additionally

[TABLE]

This makes it possible to calculate the time integral of $J^{\alpha}$, with $\alpha$ to be determined,

[TABLE]

where $\gamma=(n+1)/(n-1)$, and $|\varsigma|=1$ is a sign factor required because in this formal manipulation of exponents we have ignored the usual problems arising from complex arguments. It will be more convenient to determine it later, in the final variables $u$ and $\dot{u}$.

Taking $\alpha=1/(2\gamma)$ in the above leads again to the Euler integral and the hypergeometric functions

[TABLE]

These are two forms of the same first integral, but the former is more convenient around $\zeta=1$, which corresponds to the line $2+(n+3)y=0$, and has a singularity at $\zeta=0$; conversely, the latter works around $\zeta=0$ ($x=0$) and has a singularity at $\zeta=1$.

When re-expressing $z$ and $J$ in terms of $x$ and $y$ and then $u$ and $\dot{u}$, we again ignore the complex arguments and expand the powers freely. For $\widetilde{C}_{1}$ ($C_{1}$ is analogous) this gives

[TABLE]

with $V=4u^{n+1}+\omega_{0}^{2}(2u+(n+3)\dot{u})^{2}$. Instead of dealing with all the regions of $\zeta$ or $x$ and $y$ to determine $\varsigma$ it seems easier to differentiate $I_{1}$ and ensure that the pseudo-polynomial expression in $u$ and $\dot{u}$ vanishes identically. This happens for

[TABLE]

When $n$ is odd, $V$ is always positive, so this reduces to $\varsigma=\text{sign}(2u+(n+3)\dot{u})$; when $n$ is even, $V=0$ defines a real invariant curve $\mathcal{C}$, which lies in the left half-plane and has a cusp at the origin. $V$ is positive to the right of this curve, so (38) applies there with the same sign factor. For $C_{1}$, and consequently $I_{2}$ in (9), the factor becomes

[TABLE]

which is just $\text{sign}(u)$ for odd $n$; when $n$ is even, $I_{2}$ is to be used to the left of $\mathcal{C}$, where both $u$ and $V$ are negative, so $\varsigma^{\prime}$ simplifies to $(-1)^{n/(n+1)}$ as in the theorem. Both hypergeometric functions are real when their argument is less than 1, so only in the last case does the integral acquire a (constant) complex phase via the roots of $V$.∎

An interesting difference between between this and Lemma 1 is that we have only used a relation between $K_{1}-K_{2}$ and $F_{2}-F_{1}$, without assuming that the latter is linked with $F_{3}$. Because $J$ involved only $F_{1}$ and $F_{2}$ that was enough to express the integrand as a function of $\zeta$ only. This observation will be crucial in proving a more general lemma in section 4.

It should also be noticed that the Gaussian integral happens to be expressible via the incomplete beta function, namely

[TABLE]

Together with the special case $n=3$ discussed in the introduction, this suggests a link with inverse elliptic functions or even reduction to elementary functions, which is discussed for all three systems in Section 5.

An example of the phase portrait as expressed via levels of the first integrals for odd $n$ is shown in Figure 1. As the origin is an attracting node, there cannot exist a global analytic first integral. Upon passing between the half-planes, the values of the two integrals $I_{1}$ and $I_{2}$ need to be adjusted if one wants to keep the formulae intact, alternatively they can be regarded as multivalued functions. The slightly more complicated situation with even $n$ is shown in Figure 2, in addition to the line of discontinuity there appears an invariant curve with another critical point on it.

Finally, one has to remember that to generalize the $u^{3}$ term in (1) to arbitrary exponents, while keeping the “harmonic” character of the force, one should take $|u|^{n-1}u$ instead of just $u^{n}$. The argument of the hypergeometric function, which corresponds to $\zeta$, then satisfies $0\leq 4|u|^{n+1}/V\leq 1$ everywhere except the origin, and the situation is again like that in Figure 1.

3

All of the above can be immediately applied also to the Duffing-van der Pol equation (10), although just the existence of a time-dependent first integral will not be enough. As in [8], the exponents will have to satisfy $m=n-1$, so that excludes the classical van der Pol ($m=2$, $n=1$), and Duffing ($m=0$, $n=3$) oscillators.

Proof of Theorem 2. Adopting the same variables as before, i.e, $x=u^{n-1}$ and $y=\dot{u}/u$, the system is now

[TABLE]

The Darboux polynomial $F_{1}=x$ is self-evident, and if the condition

[TABLE]

holds, a second linear one can be found by direct computation:

[TABLE]

These two are enough for the construction of (11), because we have

[TABLE]

However, this alone is not enough to use Lemma 1 because the cofactor of $F_{1}/F_{2}$, or a similar combination, is a function of $y$ alone so not necessarily a Darboux polynomial. At the same time, $F_{2}-\beta^{2}F_{1}$, which could be the candidate for $f_{2}-f_{1}$, is linear in $y$ but in order for it to be a Darboux polynomial, an additional condition is necessary:

[TABLE]

It is independent of the one in (43), and it guarantees the existence of the Darboux pair

[TABLE]

Notice that just $F_{1}$ and $F_{3}$ are insufficient for the time-dependent integral under consideration, because no linear combination of their cofactors can be made constant. It should also be said that unless a full characterization of the system’s Darboux polynomials is given, it remains an open question if the condition (43) alone is not enough to proceed with the proof.

If both (43) and (46) are to be satisfied, it follows that $\beta=n+1$ and the conditions (12) are obtained. The situation thus resembles that of the Duffing case because the value of $\omega_{0}$ is strictly determined by the degree $n$.

Is is now straightforward to take

[TABLE]

and check that Lemma 1 can be applied with $\zeta=f_{1}/f_{2}$ and $J=f_{1}^{1/n}f_{3}^{1-1/n}$, which leads to the following Euler integrals

[TABLE]

A difference with the previous situation is that $|\zeta|$ might be greater than 1, even when we take $x=|u|^{n-1}$. Still, the above hypergeometric functions can be continued along the real axis past $\zeta=-1$ and the only potential problematic points are $\zeta=1$ and $\zeta=\infty$. The former corresponds to $f_{1}=f_{2}\iff F_{2}=0$, and the latter to $f_{2}=0\iff F_{3}=0$, so they are both invariant sets. In the original variables, they become $\omega_{0}^{2}(u+\beta\dot{u})+u^{n}=0$ and $\beta\dot{u}+u=0$, respectively. These curves separate the phase space into regions in which one of the above hypergeometric function can be chosen. Finally, substituting for $\zeta$, $x$ and $y$, the sign factors can be determined as in the previous proof, and turn out not to depend on the region giving the stated formulae. Lastly, the integrals might acquire a complex phase due to fractional powers of $V=u^{n}+\omega_{0}^{2}(u+\beta\dot{u})$, but since this is just a multiplicative constant, its absolute value can be taken on each side of $V=0$.∎

Here again, the integral itself can be written as the incomplete beta function

[TABLE]

although it is more convenient to use the Gauss function, so that the fractional power of $V$ can be written as a common factor. In fact, this beta function turns out to be elementary, albeit with no simple general formula, as is shown in Section 5.

Like before, for non-integer $n$, or to keep the sign as in the harmonic case, $u^{n}$ can be replaced by $|u|^{n-1}u$. A plot of the level sets of $I_{3}$ and $I_{4}$ is presented in Figure 3. Both first integrals are valid everywhere between the red curves, and the two formulae merely reflect the fact that they are not always real-valued. As opposed to the previous system, there is no need for the regions of definition to overlap, because the red curves, being invariant sets, are impassable barriers.

4

The third system presents a new challenge, because of non-polynomial Darboux elements. To wit, in the coordinates $x=u^{n-1}$, and $y=\dot{u}/u$ it reads

[TABLE]

If, in addition to both the previous conditions $\beta=n+1$, and $\omega_{0}^{2}=n(n+1)^{-2}$, the new parameter satisfies

[TABLE]

then, there exist three Darboux pairs:

[TABLE]

Moreover, they generate a time-dependent first integral through

[TABLE]

where the exponents were chosen to be

[TABLE]

The appearance of an exponential element along with polynomials is not surprising: these are the only two types that can arise for a two-dimensional system with a Liouvillian integral, as proved by Christopher [16]. It is thus natural to include exponentials, but at the same time one would wish for a criterion or procedure which, like Lemma 1, works in any dimension.

The first step in getting there is to recall that the cofactor $k_{1}-k_{2}$ had to be reexpressed as a function of $f_{i}$; likewise, the time-dependent integral contained $f_{3}$ which had to be eliminated. In other words, we are really working with Darboux elements like $f_{1}/f_{2}$, $J$ or $(f_{1}-f_{2})^{\alpha_{3}}$, and they can be grouped as we please, with the aim of integration in the variable $\zeta$. This suggests the following generalization:

Lemma 2**.**

*Let $D=\sum_{i=1}^{d}P_{i}\partial_{i}$ be a derivation which admits at least two Darboux elements which can be combined to yield Darboux pairs $\{f_{1},k_{1}\}$ and $\{f_{2},k_{2}\}$ such that:

1) the element $J:=f_{1}f_{2}$ has a constant cofactor $\alpha_{0}$;*

2)* the cofactor of $\zeta:=f_{1}/f_{2}$ satisfies*

[TABLE]

for some Liouvillian $L$. Then, the dynamical system $\dot{x}=P_{i}(x)$ associated with the derivation $D$ has a Liouvillian first integral.

In particular, a binomial $L(\zeta)=\zeta^{a}(1-\zeta)^{b}$ leads to the Gauss function, while $L(\zeta)=\zeta^{a}\mathrm{e}^{b\zeta}$ to the Kummer function in the first integral.

Proof. The element $J$ considered as a function of time satisfies $J(t)=\mathrm{e}^{\alpha_{0}t}J(0)$, so the integral of $J^{\gamma}$ can be formally transformed, using $\zeta=f_{1}/f_{2}$ as follows

[TABLE]

The same considerations of complex arguments as in Lemma 1 apply, so a suitable constant $\varsigma$ of modulus 1 will have to be added once a particular region of $\zeta$ is fixed. The two special cases then give the time-independent first integral $I$ through

[TABLE]

or

[TABLE]

where $\Gamma(s,z)$ is the incomplete gamma function. ∎

Both the previous proofs are special cases, and they are important preliminary results showing how to combine the Darboux polynomials. In Lemma 1, denoting the “old” polynomials by $g_{i}$, we can take $f_{1}=g_{1}^{(\alpha_{1}+1)/2}g_{2}^{(\alpha_{2}-1)/2}g_{3}^{\alpha_{3}/2}$, and $f_{2}=g_{1}^{(\alpha_{1}-1)/2}g_{2}^{(\alpha_{2}+1)/2}g_{3}^{\alpha_{3}/2}$, so that $\zeta=f_{1}/f_{2}=g_{1}/g_{2}$, and $J=f_{1}f_{2}=g_{1}^{\alpha_{1}}g_{2}^{\alpha_{2}}g_{3}^{\alpha_{3}}$, while the relation between the cofactors yields

[TABLE]

when $\alpha_{1}+\alpha_{2}+\alpha_{3}=1$ as in the proof of Lemma 1. Similarly in the proof of Theorem 1, we could take $f_{1}=F_{1}F_{2}^{(n-3)/4}$, $f_{2}=F_{2}^{(n+1)/4}$, and then

[TABLE]

The crucial quantity is always $k_{1}-k_{2}$, which is polynomial in the dynamical variables, but not necessarily a polynomial, or even algebraic, in $f_{i}$. This is ostensibly so in the present system (51), to which we now turn.

Proof of Theorem 3. The exponential factor suggests the choice of $\zeta=F_{1}/F_{2}$, and Lemma 2 can be used with

[TABLE]

The relation between the cofactors is

[TABLE]

so, by Lemma 2, the first integral is

[TABLE]

or, going back to the original variables,

[TABLE]

where $V=u+\beta\dot{u}+2\varphi u^{n}$. The curve given by $V=0$ is invariant, so the dynamics is separated into two regions, in which the Kummer function is analytic and real. To get a single real formula it is thus enough to take the absolute value in the $V^{-1/n}$ factor. Differentiating then gives $\varsigma=1$ as in the theorem.∎

Phase portraits for this system with $n=3$ and $n=4$ are shown in Figure 4. As before, if one wishes to consider arbitrary $n$, or just to keep a single critical point at the origin, a change of $u^{n}$ to $|u|^{n-1}u$ is necessary. The invariant curve corresponds to the irregular singularity of the Kummer function at infinity, so a series expansion around it is not readily available due to the Stokes phenomenon. Still, as with the exponential function, the series around zero has infinite radius of convergence.

5

Having obtained the hypergeometric expressions, the natural question to ask is whether they can be reduced to elementary functions. In the case of ${}_{2}F_{1}$, the answer follows from a result of Chebyshev’s [17]:

Theorem 4**.**

If $p$, $q$ and $r$ are rational numbers and $a$ and $b$ are nonzero real numbers, the indefinite integral $\int z^{p}(a+bz^{r})^{q}\mathrm{d}z$ is elementary if and only if at least one of $(p+1)/r$, $q$, or $q+(p+1)/r$ is an integer.

In the first system, $a=-b=r=1$ while $p=-1+1/(n+1)$ and $q=-1/2$, and the three numbers to check are $1/(n+1)$, $-1/2$, and $-1/2+1/(n+1)$ – none of them are integers when $n>1$, so the hypergeometric function is truly transcendental in this case.

In the second system, $b=-a=r=1$ while $p=-1+1/n$ and $q=-1/n$, which means that $q+(p+1)/r=1$ regardless of $n$. However, the theorem still requires that the exponents be rational, and the specific change of variables that make the integrand rational depends on those exponents. In the central case of integer $n>1$ the change is $\zeta=1/(1+z^{n})$, and the relevant expression is

[TABLE]

The evaluation of such an integral is straightforward for a specific value of $n$, e.g. $n=2$ gives $I_{z}=-2\text{arcsin}\left(\sqrt{1-\zeta}\right)$, but it cannot be given as a simple expression for symbolic $n$ – other than by using the hypergeometric or beta functions.

Additionally, the transcendental functions of the first theorem can become inverses of elliptic functions. The three cases when this happens for $B_{\zeta}\left(\frac{1}{n+1},\frac{1}{2}\right)$ are given in [11], and correspond to $n=2,3,5$. This is best seen on the initial case, which admits a change to the Hamiltonian system (5), and can be integrated for energy $E$ by

[TABLE]

where, in the transformed variables, $\zeta=W^{4}/(4E)$. This is an inversion of the direct solution of the Hamiltonian system in terms of the Jacobian elliptic function $W=\sqrt{2}E^{1/4}\mathrm{sn}(E^{1/4}Z)$.

Finally, the Kummer function of Theorem 3 is transcendental for integer $n>1$, because it can be rewritten as the classical exponential integral

[TABLE]

where $\zeta=s^{n}$. The Gauss error function is a familiar example of the above for $n=2$.

Acknowledgements

This work was supported partially by the grant No. DEC-2013/09/B/ST1/04130 of the National Science Centre of Poland, and partially by the Japan Society for the Promotion of Science, Grant-in-Aid for Scientific Research (B) (Subject No. 17H02859).