Qualitative properties of systems of 2 complex homogeneous ODE’s: a connection to polygonal billiards

In the following, we study a special case of the system of 4 real ordinary differential equations with homogeneous right-hand sides:

$$\dot{y}_{i}=P_{i}^{(r)}(y_{1},y_{2},y_{3},y_{4})\qquad(1\leq i\leq 4),$$

(1)

where the $P_{i}^{(r)}$ are homogeneous polynomials of degree $r$ in the 4 variables $y_{i}$. Such equations are a bit anomalous due to the absence of linear terms, nevertheless they are studied to a considerable extent in various fields, such as reaction kinetics and ecology.

The case we shall study is special in 2 respects

1. 1.

   First, and most importantly, the system (1) is assumed to arise from a complex system of 2 ODE’s:

   	$\displaystyle\dot{x}_{1}$	$\displaystyle=$	$\displaystyle p_{r}(x_{1},x_{2})$		(2a)
   	$\displaystyle\dot{x}_{2}$	$\displaystyle=$	$\displaystyle q_{r}(x_{1},x_{2}).$		(2b)

   Here $p_{r}(x,y)$ and $q_{r}(x,y)$ are both homogeneous polynomials of degree $r$.

2. 2.

   Second, the coefficients of the polynomials $p_{r}(x_{1},x_{2})$ and $q_{r}(x_{1},_{2})$ are real.

Clearly, the most important restriction is the first one. The second can certainly be significantly generalised. However, the first one is essential, and is responsible for the highly atypical behaviour we identify for such systems.

As we shall see, these systems can be completely understood, once a corresponding polygonal billiard problem is solved: each orbit of (2) can be brought uniquely into correspondence with the orbit of a given polygonal billiard, the characteristics of which do not depend on the initial conditions, but only on the system (2) itself.

In particular, we find the following properties, which differ rather strongly from the generic properties of such ODE’s: first and foremost, the dynamics is regular: that is, with probability one, the orbits are defined and finite for all times. As follows, for example, from escape , within the sets of quadratic systems, there is an open set with the property that the time for a solution to diverge is finite for an open set of initial conditions. On the other hand, for systems of the type (2), solutions are non-singular for all times with probability one. Note that this does not mean that the orbit remains bounded: in a rather general case, as we shall see, the orbit never diverges, but with probability one assumes arbitrarily large values over small time intervals.

In a polygonal billiard, the singular orbits (that hit a corner) have an important characteristic: in their vicinity, the regular orbits vary discontinuously: Figure 1 shows how this occurs for an orbit hitting a $2\pi/3$ corner, but the effect occurs generally, except for angles of the form $\pi/n$. In the ODE system, this is reflected in the fact that the orbits hitting a corner diverge, so that the theorems on continuous dependence of the solution of an ODE on initial conditions, fail in their vicinity.

Another striking feature of these systems is the fact that their Lyapunov exponents chaos1 ; chaos are always zero. This means that the dynamics can actually be predicted efficiently over large times, though, as we shall see, the dynamics can indeed be quite complicated.

We thus show that a correspondence exists between any given system of ODE’s of the type (2) and the dynamics of a free particle bouncing elastically within a polygon, the shape of which is uniquely and elementarily determined by the system (2). Since many results on the existence and nature of periodic orbits, on ergodicity, Lyapunov exponents and several other properties, are known for polygonal billiards, it turns out that they trivially translate into corresponding properties for the system of ODE’s (2). This is the essence of this paper.

In Section II we show how the system (2) can be solved by quadratures. Since the integrals involved are complicated, inverting them is not a trivial task, and an understanding of the orbit requires an additional remark. In Section III, we establish in detail the correspondence, in Section IV we establish the general consequences of this correspondence. The results are different depending on whether the polygon is bounded or unbounded, so this Section is divided in 3 subsections, which treat the properties valid in either case (Subsection IV.1) in the unbounded case, in which the orbits are scattering orbits (Subsection IV.2) and finally the most important csase in which the polygon is bounded and the motion is finite(Subsection IV.3). In Section V, we illustrate numerically some of the predictions made here, and in Section VI we present conclusions , Finally, to keep this article self-contained, we quote without proof the various properties of polygonal billiards used in this paper in Appendix A.

We consider the system (2) of ordinary differential equations. These are viewed as complex equations, that is, viewed as a system involving real quantities, they correspond to a system of 4 equations in 4 unknowns, corresponding to the real and imaginary parts of $x_{1}$ and $x_{2}$. We rewrite these equations as follows:

	$\displaystyle\dot{x}_{1}$	$\displaystyle=$	$\displaystyle x_{2}^{r}P_{r}(x_{1}/x_{2})$		(3a)
	$\displaystyle\dot{x}_{2}$	$\displaystyle=$	$\displaystyle x_{2}^{r}Q_{r}(x_{1}/x_{2}).$		(3b)

Here $P_{r}(u)$ and $Q_{r}(u)$ are both polynomials of degree $r$. There is a minor loss of generality in this description, as we assume that there exists in both equations a term $x_{2}^{r}$. This can always be reached by a linear transformation of the dependent variables. By an appropriate scaling of $x_{2}$ we may further choose $Q_{r}(u)$ to be a monic polynomial.

In the following, we present a solution by quadratures of this system. This is not new, but has been explicitly formulated by nicklason , and similar calculations have been performed by Garnier garnier1 ; garnier2 . This approach has also been used in CCL to identify particularly simple special cases of (2) for $r=2$.

We define the following quantities

	$\displaystyle u$	$\displaystyle=$	$\displaystyle x_{1}/x_{2},$		(4a)
	$\displaystyle x_{1}$	$\displaystyle=$	$\displaystyle uR(u),$		(4b)
	$\displaystyle x_{2}$	$\displaystyle=$	$\displaystyle R(u).$		(4c)

From (3) follows the following equation for $u$

$$\dot{u}=-x_{2}^{r-1}\left[uQ_{r}(u)-P_{r}(u)\right]=:-x_{2}^{r-1}S_{r+1}(u)$$

(5)

where the final equation defines $S_{r+1}(u)$, which is a monic polynomial of degree $r+1$. Note the use we have made of the normalisation of $Q_{r}(u)$ as a monic polynomial.

From (5) and (4c), one obtains

$$\dot{u}=-R(u)^{r-1}S_{r+1}(u).$$

(6)

Now from (3b) and (4c) one finds

	$\displaystyle\dot{x}_{2}$	$\displaystyle=$	$\displaystyle R^{\prime}(u)\dot{u}$		(7)
		$\displaystyle=$	$\displaystyle-R(u)^{r-1}R^{\prime}(u)S_{r+1}(u)$
		$\displaystyle=$	$\displaystyle R(u)^{r}Q_{r}(u).$

This in turn leads to a separable equation for $R(u)$:

$$\frac{R^{\prime}(u)}{R(u)}=-\frac{Q_{r}(u)}{S_{r+1}(u)}.$$

(8)

Let $u_{\alpha}$, $0\leq\alpha\leq r$ be the zeros of $S_{r+1}(u)$, that is:

$$S_{r+1}(u)=\prod_{\alpha=0}^{r}(u-u_{\alpha}).$$

(9)

We decompose the right-hand side of (8) in partial fractions, assuming that none of the $u_{\alpha}$ are double zeros:

$$\frac{Q_{r}(u)}{S_{r+1}(u)}=\frac{1}{r-1}\sum_{\alpha=0}^{r}\frac{\mu_{\alpha}}{u-u_{\alpha}},$$

(10)

where the $1/(r-1)$ prefactor is introduced for future convenience. Matching the $u\to\infty$ behaviours of both sides of (10), remembering that both $Q_{r}(u)$ and $S_{r+1}(u)$ are monic, we obtain

$$\sum_{\alpha=0}^{r}\mu_{\alpha}=r-1.$$

(11)

Note that the $\mu_{\alpha}$ and the $u_{\alpha}$ are altogether independent of the initial conditions and instead characterise the system (2) itself.

(8) is now immediately integrated to yield

$$R(u)=C\prod_{\alpha=0}^{r}\left(u-u_{\alpha}\right)^{-\mu_{\alpha}/(r-1)}.$$

(12)

Here $C$ is an integration constant determined by the relation

$$x_{2}(0)^{2}=C\prod_{\alpha=0}^{r}\left[x_{1}(0)-u_{\alpha}x_{2}(0)\right]^{-\mu_{\alpha}/(r-1)}.$$

(13)

We now proceed to a final normalisation step: the solutions of (2) can always be scaled by a fixed real factor $\lambda$, which corresponds to a scaling of $t$ by the factor $\lambda^{r-1}$. We may hence, without loss of generality, scale the initial conditions accordingly and therefore fix the norm of $C$. We thus set $|C|=1$ and

$$C=-e^{i{\chi_{0}}}.$$

(14)

We may now determine the time-dependence of $u$ using (6):

$$\dot{u}=e^{i{\chi_{0}}}\prod_{\alpha=0}^{r}\left(u-u_{\alpha}\right)^{-\mu_{\alpha}+1}.$$

(15)

which leads to the expression via quadratures

$$t=e^{-i{\chi_{0}}}\int_{u(0)}^{u}\prod_{\alpha=0}^{r}\left(u^{\prime}-u_{\alpha}\right)^{\mu_{\alpha}-1}du^{\prime}$$

(16)

where $u(0)=x_{1}(0)/x_{2}(0)$.

We now limit ourselves to the subclass of systems in which the coefficients of the polynomials $P_{r}(u)$ and $Q_{r}(u)$ are all real. Under these conditions, the fact that all $u_{\alpha}$ should be real and simple, is no more exceptional. Indeed, given a system with that property, all other systems that are sufficiently close also have this property, so that we are in a generic case.

The essential observation we now make is the following: if all $u_{\alpha}\in\mathbb{R}$, all $\mu_{\alpha}\in\mathbb{R}$ as well. It then follows that, for appropriate values of the $\mu_{\alpha}$, specifically for $0\leq\mu_{\alpha}\leq 1$, the transformation defined by (16) is the conformal map from the upper half-plane to a finite, convex polygon ${\cal P}$, having $r+1$ sides and interior angles $\mu_{\alpha}\pi$, the well-known Schwarz–Christoffel transformation schwarz1 ; schwarz2 ; schwarz . It follows immediately from (11) that the sum of the interior angles of the polygon is, as it must be, equal to $(r-1)\pi$. Further note that the polygon’s shape, which is the main object of our consideration, is determined both by the interior angles given by the $\mu_{\alpha}$, and by the relative lengths of the sides, determined by $r-3$ values of the $u_{\alpha}$. We denote by $v_{\alpha}$ the vertices of ${\cal P}$ corresponding to $u_{\alpha}$, by ${\cal S}_{\alpha}$ the side of ${\cal P}$ connecting $v_{\alpha}$ to $v_{\alpha+1}$, where $\alpha+1$ is computed modulo $r+1$. To ${\cal S}_{\alpha}$ corresponds in the boundary of the upper half-plane, the interval $I_{\alpha}=[u_{\alpha},u_{\alpha+1}]$.

Note in passing that the task we face here is different from, and in many ways easier than, the one usually solved by the Schwarz–Christoffel transformation: normally one is given a polygonal domain and looks for a conformal transformation. In that case, the determination of the $u_{\alpha}$ can be challenging schwarz1 . In our case, we are given the $u_{\alpha}$ and the angles $\mu_{\alpha}$, and our task is merely to determine the image of the upper half-plane under this transformation.

For definiteness’s sake, let us assume the initial condition $u(0)$ to be in the upper half-plane (the opposite case is similar). The map

$$\Phi(u)=\int_{u(0)}^{u}du^{\prime}\,\prod_{\alpha=0}^{r}\left(u^{\prime}-u_{\alpha}\right)^{\mu_{\alpha}-1}.$$

(17)

maps the upper half-plane onto the inside of a convex $r$-sided polygon $\cal{P}$ containing the origin. We have in particular

$$\Phi(u_{\alpha})=v_{\alpha}\qquad(0\leq\alpha\leq r).$$

(18)

The equation (16) means that the straight line $\cal{L}$ defined by $e^{i{\chi_{0}}}t$, for all real $t$, is the image of the orbit $u(t)$ under the map $\Phi$.

We must therefore determine the inverse image of $\cal{L}$ under $\Phi$. For the segment of $\cal{L}$ that lies entirely in $\cal{P}$, the corresponding part of the trajectory lies wholly in the upper half-plane. As the line leaves $\cal{P}$ by the side $\cal{S_{\alpha}}$ corresponding to the real line interval $I_{\alpha}=[u_{\alpha},u_{\alpha+1}]$, the corresponding orbit of $u$ leaves the upper half-plane by the interval $I_{\alpha}$. Due to the Schwarz reflection principleconway , the image under $\Phi$ of the sheet which the orbit $u(t)$ enters, is the polygon described by the reflection of $\cal{P}$ on the side $\cal{S_{\alpha}}$. We may therefore keep the orbit of $u$ inside the upper half-plane by specularly reflecting it with respect to the real axis, and correspondingly keep the line $\cal{L}$ inside the polygon $\cal{P}$, also by reflection. Indefinite repetition of this procedure, for both positive and negative times, leads to a billiard orbit inside ${\cal P}$. For an illustration of the way this proceeds, see Figure 2. Note that this construction is rather similar to one used in GS for a somewhat related problem. If we now take the inverse image under $\Phi$ of this billiard orbit, we obtain the orbit $u(t)$ specularly reflected each time it crosses the real axis, from which the actual $u(t)$ orbit is readily reconstructed.

At this stage let us define some additional notation: the initial segment of the billiard trajectory is the straight line segment $e^{i\chi_{0}}t$. After $n$ bounces we define the corresponding straight line segment to be $e^{i\chi_{n}}(t-\tau_{n})$. Here $\tau_{n}$ is the time at which the orbit hits ${\cal P}$ and begins the $n$-th bounce. As the billiard orbit is successively followed, the connection between $t$ and $u$ given by (16) is modified to

$$t-\tau_{n}=e^{-i{\chi_{n}}}\int_{u(0)}^{u}\prod_{\alpha=0}^{r}\left(u^{\prime}-u_{\alpha}\right)^{\mu_{\alpha}-1}du^{\prime}.$$

(19)

We therefore see that the inverse image under $\Phi$ of a billiard orbit of $\cal{P}$ is the orbit of $u$ reflected back into the upper half-plane each time it hits the real axis. Since $\Phi$ is not an easily determined map, this does not represent an exact solution, but remembering the many results known about polygonal billiards, the correspondence yields several non-trivial results concerning the solutions’ qualitative behaviour.

Before we proceed to describe these, however, it is of some importance to extend the validity of the correspondence as far as possible. In the case $0\leq\mu_{\alpha}<1$, the image of the upper half-plane is an $r$-sided convex polygon. If we generalise this to $0\leq\mu_{\alpha}<2$, we obtain arbitrary bounded polygons, whether convex or not: the polygon’s interior angles are then $\mu_{\alpha}\pi$, and they still add up to $(r-1)\pi$. The extension to negative values of $\mu_{\alpha}$ leads to unbounded polygons. Due to (11), negative values of $\mu_{\alpha}$ must always coexist with positive ones. If $\mu_{\alpha}<0$, the integral describing $\Phi$ diverges as $u\to u_{\alpha}$ on $\mathbb{R}$, both from the right and the left. The polygon’s boundary thus contains two lines diverging to infinity and forming an angle $\mu_{\alpha}\pi$. We may thus draw a polygon corresponding to all real values of $\mu_{\alpha}$ satisfying both (11) and $-2<\mu_{\alpha}<2$.

Further extensions, whether to values of $\mu_{\alpha}$ with $|\mu_{\alpha}|\geq 2$ or complex values of $\mu_{\alpha}$ may well be possible, but it is not obvious how to extend the above construction to such cases, and more generally speaking, how to obtain meaningful results from them.

In the following we illustrate using numerical simulations some of the findings described in Section IV. All the simulations are performed directly on the system (2), without using the results of Section II.

The system is constructed from the given data $\mu_{\alpha}$, $1\leq\alpha\leq r+1$, which vary from system to system, as follows: the $u_{\alpha}$ are always conventionally taken to be

$$u_{\alpha}=\alpha-1/2-\left\lfloor{\frac{r}{2}}\right\rfloor$$

(33)

and the polynomial $S_{r+1}(u)$ is computed from (9), from which $Q_{r}(u)$ and from that eventually $P_{r}(u)$ are computed using (10). The initial conditions are taken with a random, or generic, value of $\chi_{0}$ and, if not otherwise stated, with a value of $u(0)=1/4$ always different rom $u_{\alpha}$. Since $|C|=1$, we can fully determine the initial conditions. If not stated otherwise, we shall always be dealing with the case $r=2$.

First let us show periodic orbits. As we saw, whenever the $u_{\alpha}$ are real and $\chi_{0}=\pi/2$, the resulting orbit is almost surely periodic. Further, they vary continuously as $u(0)$ varies, apart from an obvious discontinuity when $u(0)$ crosses a $u_{\alpha}$, since this corresponds to a corner. We show this in Figure 4, where we look at a a rational case $\mu_{0}=1/2,\mu_{1}=3/7$ and $\mu_{2}=3/14$ which yields the equations

	$\displaystyle\dot{x}_{1}$	$\displaystyle=$	$\displaystyle\frac{3}{14}x_{1}^{2}+\frac{5}{14}x_{1}x_{2}-\frac{3}{8}x_{2}^{2},$		(34)
	$\displaystyle\dot{x}_{2}$	$\displaystyle=$	$\displaystyle x_{1}^{2}-\frac{9}{7}x_{1}x_{2}+\frac{3}{28}x_{2}^{2}.$		(35)

We proceed to display ergodicity. If the triangle’s angles, that is the $\mu_{\alpha}$, are rational, then almost every direction is ergodic. The places where the orbit is reflected on the triangle ${\cal P}$ are thus uniformly distributed. Translating this to the $u$ variables, this means that the values of $u$ where the orbit crosses the real axis have the probability distribution

$$p(u)=\frac{1}{\cal N}\prod_{k=0}^{r}(u-u_{\alpha})^{\mu_{\alpha}-1},$$

(36)

where $\cal N$ is the normalisation. An example of the histogram for the $u$ values of the real crossings of a single orbit over a time of $5\cdot 10^{4}$ is given, together with the predicted distribution (36). We see a good agreement in the case described by Figure 5 in the rational case shown in (35). A similarly good agreement (not reported) is found for the irrational case $\mu_{0}=1/\sqrt{5}$, $\mu_{0}=1/\sqrt{7}$ and $\mu_{2}=1-1/\sqrt{5}-1/\sqrt{7}$.

On the other hand, in Figure 6, we display evidence for one of the basic differences between rational and irrational angles: in Figure 6, we display the values of $x_{2}(t)$ on the complex plane at those times in which $u(t)=x_{1}(t)/x_{2}(t)$ crosses the real axis for one single long orbit. Indeed, an arbitrary orbit lies in the 3-dimensional subspace of possible values of $x_{1}$ and $x_{2}$ defined by the equation $|C|=1$. The intersection of the orbit with the 2-dimensional space defined by imposing the additional condition $\mbox{Im}\,x_{1}(t)/x_{2}(t)=0$ are thus isolated points. Apart from this we have no further indication. In the two plots of this nature shown in Figure 6, the lower one, corresponding to the case of irrational values of $\mu_{\alpha}$, indeed shows a set of points more or less randomly scattered on the plane. However, this is definitely not the case for the upper diagram, which corresponds to simple rational values of the $\mu_{\alpha}$. There the set of points is essentially a set of curves, in other words, it is one-dimensional. This corresponds, of course, to the fact that in the corresponding orbit in the triangular billiard, the direction of the orbit can only take a finite number of values trian1 .

Finally, let us verify the validity of the remarks made at the end of Section IV concerning large peaks in the values of $x_{1,2}(t)$. We consider the case $r=2$, and the rational case discussed above. We find as a histogram of the absolute value of $x_{1}(t)$ taken at unit time intervals for an orbit of duration $5\cdot 10^{4}$. The existence of a power-law is undeniable, and the agreement with the theoretical prediction of an exponent $-3$ is fairly convincing.

We identify a class of systems of 2 complex ODE’s with remarkable properties which follow from the fact that we can associate to every orbit of the system a unique orbit of a corresponding bounded polygonal billiard. From this identification follow various remarkable qualitative properties: the Lyapunov exponents are all zero, the motion almost surely never diverges and remains bounded by a constant $B$ for a fraction of the time that goes to one as $B\to\infty$.

While the results are rather special, being limited to systems of 4 real homogeneous ODE’s derived from a complex analytic system, they can be significantly extended: since the properties here described are qualitative in nature, they extend to every system that can be obtained from (2) via a change of variables. As a trivial example, using real linear transformations, it is possible to obtain homogeneous systems of 4 ODE’s for which the analyticity property is hidden. Similarly, all non-linear transformations which preserve the homogeneity property can also be used to extend the relevant class.

Similarly, starting from the complex Newtonian equation

$$\ddot{z}=z^{k}$$

(37)

we obtain by the transformationCCL

	$\displaystyle x_{1}$	$\displaystyle=$	$\displaystyle z^{(k-1)/2},$		(38)
	$\displaystyle x_{2}$	$\displaystyle=$	$\displaystyle\frac{\dot{x}_{1}}{x_{1}}.$		(39)

the set of complex equations

	$\displaystyle\dot{x}_{1}$	$\displaystyle=$	$\displaystyle x_{1}x_{2},$		(40)
	$\displaystyle\dot{x}_{2}$	$\displaystyle=$	$\displaystyle x_{1}^{2}+\frac{2}{1-k}x_{2}^{2}.$		(41)

These belong to our class for all real values of $k\neq 1$, so the various results derived above, concerning the existence of periodic orbits, the vanishing of the Lyapunov exponent and so on, all follow for this Newtonian equation. Note that, in this case, the existence of a solution in terms of quadratures follows trivially from energy conservation, but this solution leads to hyperelliptic integrals for which it is not straightforward to obtain the various results stated above.

Other extensions are possible. In particular, it is possible to extend the solution by quadratures to the case of a set of 2 homogeneous complex ODE’s with a linear term of the following form

	$\displaystyle\dot{x}_{1}$	$\displaystyle=$	$\displaystyle-\alpha x_{1}+p_{r}(x_{1},x_{2})$		(42a)
	$\displaystyle\dot{x}_{2}$	$\displaystyle=$	$\displaystyle-\alpha x_{2}+q_{r}(x_{1},x_{2}).$		(42b)

However, the nature of the billiard motion is significantly different and its study is left for future work.

I would like to acknowledge financial support by UNAM PAPIIT-DGAPA-IN113620 as well as by CONACyT Ciencias Básicas 254515.

Data sharing not applicable—no new data generated.