BURSTING SOLUTIONS OF THE RÖSSLER EQUATIONS

Abstract

We provide an analytic solution of the Rössler equations based on the asymptotic limit $c\to \infty $ and we show in this limit that the solution takes the form of multiple pulses, similar to “burst” firing of neurons. We are able to derive an approximate Poincaré map for the solutions, which compares reasonably with a numerically derived map.

1 Introduction

Otto Rössler introduced his set of equations:

(1.1) $$ \begin{align} \begin{aligned} \dot{x}&=-y-z,\\ \dot{y}&=x+ay,\\ \dot{z}&=b+(x-c)z, \end{aligned} \end{align} $$

where a, b and c are positive constants, as a simple prototype system which has chaotic solutions despite the presence of only a single nonlinearity [15]. The system has a strange attractor which, as shown in Figure 1, has a pleasingly simple form and, in particular, is almost sheet-like. Our purpose in this note is to provide an explanation of this fact and thereby provide an analytic solution of the equations. In the course of this analysis, we will show that the solutions take on the characteristics of “bursting” action potentials, which are familiar from models of signal propagation in neurons (for example, see [12, 19]).

Figure 1 The Rössler attractor obtained by solving equation (1.1) at parameter values $a=0.2$ , $b=0.2$ , $c=10$ .

The Rössler system has generated huge interest which continues apace; much of this work is numerically based (for example, [1, 10, 18]), but there have been several attempts to provide approximate methods for solving the equations. The first of these is the paper by Nauenberg [14]; in his study, he considers an asymptotic limit for equation (1.1) in which $c=O(1)$ , but $a,b\ll 1$ . After rescaling z with b, the equations for x and y take the form of a weakly perturbed oscillator, and so long as (the rescaled) $z=O(1)$ , the model can be solved by averaging, which produces a one-dimensional unimodal map of the type found numerically. Nauenberg emphasises that his method only works for the relatively flat parts of Figure 1 and not for what we will call the pulses, where z becomes large. Maris and Goussis [11] also discuss asymptotic methods based on the discrepancy between slow and fast time scales, but without basing their discussion on specific asymptotic limits. Cheng and Chen [2] provide approximate solutions for the pulses (and some of their results mirror our own), but again without reference to specific asymptotic limits: the approximations are made intuitively.

In anticipation of our analysis below, we begin by rescaling the variables $x,y,z\sim c$ , so that equation (1.1) takes the form

(1.2) $$ \begin{align} \begin{aligned} \dot{x}&=-y-z,\\[-2pt] \dot{y}&=x+ay,\\[-2pt] \varepsilon\dot{z}&=\varepsilon^2b+(x-1)z, \end{aligned} \end{align} $$

where $\varepsilon =1/c$ . The system in equation (1.2) has two fixed points $P_\pm $ , where

(1.3) $$ \begin{align}x=az_\pm,\quad y=-z_\pm,\quad z_\pm=\frac{1}{2a}[1\pm(1-4a\varepsilon^2)^{1/2}],\end{align} $$

so that when c is large,

(1.4) $$ \begin{align}z_+\approx \frac{1}{a}-\varepsilon^2,\quad z_-\approx\varepsilon^2.\end{align} $$

The values of the parameters chosen by Rössler [15] were $a=0.2$ , $b=0.2$ , $c=5.7$ and, for these parameters, equation (1.4) provides a useful approximation to the fixed points. Linearising the system about each fixed point, we find that the eigenvalues of the solutions $e^{\lambda t}$ are approximately (with $\varepsilon \ll 1$ )

$$ \begin{align*} \lambda_+=\frac{1}{a},\quad -\varepsilon^2a,\quad \frac{1}{\varepsilon} \end{align*} $$

for the fixed point $P_+$ corresponding to $z_+$ and

$$ \begin{align*} \lambda_-=-\frac{1}{\varepsilon},\quad \pm i+\frac{1}{2} a \end{align*} $$

for $P_-$ , corresponding to $z_-$ . The large fixed point $P_+$ is an unstable saddle and will be of no further concern, while $P_-$ near the origin is an unstable saddle-focus.

As shown by Barrio et al. [1], the saddle-focus near the origin has a homoclinic orbit near the values $a=0.2$ , $b=0.2$ , $c=10$ , which is the origin of the attracting chaotic set shown in Figure 1. As can be seen, the attractor consists of a spiral sheet of trajectories which winds out in the $(x,y)$ plane before folding up and then back on to the spiral as a Möbius band. One of the noticeable things about this figure is the fact that the attractor has the form of a sheet and is thus nearly two-dimensional. This apparent reduction in dimension is the hallmark of an asymptotic limit; an entirely similar process occurs in the Lorenz equations [9], where Lorenz famously discovered an approximately one-dimensional Poincaré map. In that case, the collapse of the attractor to an approximate sheet is associated with the asymptotic limit $r\sim \sigma \gg 1$ , and this was used to derive the map by Fowler and McGuinness [4, 5]. As a consequence of the sheet-like attractor, an appropriate Poincaré map for the Rössler equations is also approximately one-dimensional. Figure 2 shows an example of such a map, defined on the Poincaré section $y=0$ , $x<0$ . The map relates successive values of $C=-x$ on the section.

Figure 2 Form of the Poincaré map for the system in equation (1.1) on the Poincaré section $y=0$ , $x<0$ . The successive values $C_n=-x_n$ on this section are plotted. The parameter values used are $a=0.2$ , $b=0.2$ , $c=10$ .

2 The limit $\boldsymbol {c\gg 1}$

Our analysis begins with the equations in the scaled form of equation (1.2); we eliminate y to form the set

(2.1) $$ \begin{align} \begin{aligned} &\ddot{x}-a\dot{x}+x=az-\dot{z},\\ &\varepsilon\dot{z}=\varepsilon^2b+(x-1)z, \end{aligned} \end{align} $$

and consider the limiting form of the solution when $\varepsilon \ll 1$ . We will take $a,b\sim O(1)$ , although the small values $a=b=0.2$ which we use certainly help the approximation.

2.1 Slow phase, $\boldsymbol {z\ll 1}$

It is evident from Figure 1 that the solution naturally separates into two distinct phases; in the first of these, which we call the slow phase, the trajectories sit on an unstable spiral close to the plane $z=0$ . In the other phase, which we call a pulse, the spiral sheet lifts up off the z plane before being folded back down on to the spiral. The pulses occur at irregular intervals, as is clearly seen in Figure 3, and they are of relatively short duration; Figure 4 shows a close-up of one of the pulses.

Figure 3 z versus t in the solution of equation (1.2) (the rescaled version). The parameter values used are $a=0.2$ , $b=0.2$ , $c=10$ , and thus $\varepsilon =0.1$ .

The slow phase is described by the regular solution of equation (2.1), where $t\sim 1$ and $x<1$ . Here, z rapidly relaxes to a quasi-equilibrium in which

(2.2) $$ \begin{align}z\sim -\frac{\varepsilon^2b}{x-1},\end{align} $$

and consequently

(2.3) $$ \begin{align}\ddot{x}-a\dot{x}+x\approx 0,\end{align} $$

with solution

(2.4) $$ \begin{align}x=-Ce^{at/2}\cos(\omega t+\chi),\quad \dot{x}=Ce^{at/2}\sin\omega t,\end{align} $$

where we define

$$ \begin{align*} \tfrac{1}{2} a=\sin\chi,\quad \omega=\big(1-\tfrac{1}{4} a^2\big)^{1/2}=\cos\chi. \end{align*} $$

This solution remains valid so long as $x<1$ and provides the increasing part of the Poincaré map seen in Figure 2. Specifically (since $y\approx -\dot {x}$ when $z\ll 1$ ), this part of the map is just

(2.5) $$ \begin{align} C_n\approx e^{\pi a/\omega}C_{n-1}. \end{align} $$

2.2 Transition, $\boldsymbol {x\approx 1}$

The slow phase approximation breaks down when $x\approx 1$ and there is then a transition phase. We denote the value of t when x in equation (2.4) reaches one as $t_c$ , and the value of $\dot {x}$ there as $\beta $ . We thus have

(2.6) $$ \begin{align} x\sim 1+\beta(t-t_c)-\gamma(t-t_c)^2,\quad z\sim -\frac{\varepsilon^2b}{\beta(t-t_c)-\gamma(t-t_c)^2}\quad \text{as } t\to t_c{-}. \end{align} $$

The expression for x is just a local Taylor expansion. Because of equation (2.3),

(2.7) $$ \begin{align} \gamma=\tfrac{1}{2}(1-a\beta). \end{align} $$

A suitable balance of terms follows by defining

(2.8) $$ \begin{align}t=t_c+\sqrt{\varepsilon}T,\quad x=1+\sqrt{\varepsilon}X,\quad z=\varepsilon^{3/2}bZ,\end{align} $$

whence the equations take the form

(2.9) $$ \begin{align}\begin{aligned}\frac{1}{\sqrt{\varepsilon}}\ddot{X}-a\dot{X}+1+\sqrt{\varepsilon}X={\varepsilon}^{3/2}b\bigg[aZ-\frac{1}{\sqrt{\varepsilon}}\dot{Z}\bigg],\dfrac{}{}\\\dot{Z}=1+XZ, \qquad\qquad\quad\qquad\end{aligned}\end{align} $$

and the matching conditions are

(2.10) $$ \begin{align} X\sim \beta T-\sqrt{\varepsilon}\gamma T^2,\quad Z\sim -\frac{1}{\beta T-\sqrt{\varepsilon}\gamma T^2}\quad \text{as } T\to -\infty. \end{align} $$

We frequently refer to matching conditions in this paper: for example, here and at equations (2.13) and (2.35). The formal principle of matching uses the concept of an intermediate matching region [8], but the essence of the matter is that two approximations in adjoining regions should be the same in an overlap region. In practice, this can be effected as follows. We have an outer approximation given by equation (2.6) and an inner region which satisfies equation (2.9). We take the outer approximation and rewrite it in terms of the inner variables described by equation (2.8), and this yields equation (2.10). More subtle effects commonly occur in higher order approximations, including the possibility of transition layers when terms appear which cannot be matched, and the use of (small) translations in the time variable. At leading order, such effects can be ignored; for a classic treatment of such a problem, see Kevorkian and Cole’s treatment of the relaxation van der Pol oscillator.

In the present situation, it is clear from equation (2.9) that the solution has the form of an expansion $X\sim X_0+\sqrt {\varepsilon }X_1+\cdots $ . At leading order, $X=\beta T$ and thus Z satisfies

$$ \begin{align*} \dot{Z}=1+\beta TZ, \end{align*} $$

with solution

$$ \begin{align*} Z=\int_0^\infty\exp\big[\beta Ts-\tfrac{1}{2} s^2\big]\,ds, \end{align*} $$

which satisfies the matching condition. In turn, this solution breaks down, because

$$ \begin{align*} Z\sim\sqrt{\frac{2\pi}{\beta}}e^{\beta T^2/2}\quad \text{as } T\to\infty, \end{align*} $$

and the solution enters the pulse phase.

2.3 Pulse, $\boldsymbol {z=O(1)}$

The pulse phase is rather awkward to deal with, as we shall see, because it involves logarithmic terms in the asymptotic expansions. It is convenient to keep the same time scale T. We seek a balance in which $z\sim O(1)$ , and equation (2.9) becomes

(2.11) $$ \begin{align} \frac{1}{\sqrt{\varepsilon}}\ddot{X}-a\dot{X}+1+\sqrt{\varepsilon}X=az-\frac{1}{\sqrt{\varepsilon}}\dot{z},\nonumber\\ \dot{z}=\varepsilon^{3/2}b+Xz. \end{align} $$

It is convenient to define

(2.12) $$ \begin{align}z=e^\theta,\end{align} $$

so that correct to $O(\varepsilon ^{3/2})$ , the second equation in equation (2.11) implies

$$ \begin{align*} X=\dot{\theta}, \end{align*} $$

and the matching conditions as $T\to 0$ are ¹

(2.13) $$ \begin{align}X\sim\beta T,\quad \theta\sim -\theta_m+\tfrac{1}{2}\beta T^2,\end{align} $$

where

(2.14) $$ \begin{align}\theta_m=\ln\bigg\{\frac{1}{\varepsilon^{3/2}b}\sqrt{\frac{\beta}{2\pi}}\bigg\},\end{align} $$

and is thus formally large. Equivalently, $\theta $ satisfies the initial conditions

(2.15) $$ \begin{align}\theta=-\theta_m,\quad \dot{\theta}=0\quad \text{at } T=0.\end{align} $$

At leading order, $\ddot {X}=-\dot {z}$ , whence

$$ \begin{align*} \dot{X}=-z+\beta \end{align*} $$

to satisfy the matching condition, and thus $\theta $ satisfies

(2.16) $$ \begin{align}\ddot{\theta}+e^\theta-\beta=0,\end{align} $$

subject to the initial conditions (2.15). This equation is that of a conservative nonlinear oscillator and has a first integral

$$ \begin{align*} \tfrac{1}{2}\dot{\theta}^2+V(\theta)=E=V(\theta_m),\quad V=e^\theta-\beta\theta, \end{align*} $$

where the value of E is obtained from the matching condition. Since E is large, asymptotic solutions of this equation can be found and these were given by Fowler [3]. The analysis below closely follows this and is therefore only summarised here.

2.3.1 The rise

If we consider the graph of $V(\theta )$ when E is large, it is fairly clear that for $\theta $ near $-\theta _m$ , $V\approx -\beta \theta $ and the exponential is negligible, whereas at the maximum $\theta _M$ of $\theta $ , the opposite is true and, in particular,

(2.17) $$ \begin{align}\theta_M\approx \ln E,\quad \theta_m\approx \frac{E}{\beta}.\end{align} $$

For sufficiently small T, we can suppose the exponential is negligible and thus

(2.18) $$ \begin{align}\theta\approx -\theta_m+\tfrac{1}{2}\beta T^2,\end{align} $$

and this approximation becomes invalid when $\theta $ becomes of $O(1)$ and thus at $T\sim \sqrt {2\theta _m/\beta }$ , near where $\theta $ has a sharp peak.

2.3.2 The peak

In the vicinity of the maximum,

(2.19) $$ \begin{align}\theta=\theta_M+\phi,\quad T=T_M+e^{-\theta_M/2}\eta,\end{align} $$

which yields

(2.20) $$ \begin{align}\phi_{\eta\eta}+e^\phi=\beta e^{-\theta_M}\approx 0.\end{align} $$

The solution must match to equation (2.18) as $\eta \to -\infty $ , and if we choose to define $T_M$ and $\theta _M$ precisely via

$$ \begin{align*} T_M &=\bigg[\frac{2}{\beta}(\theta_M+\theta_m+2\ln 2)\bigg]^{1/2}\approx \sqrt{\dfrac{2\theta_m}{\beta}},\\ \theta_M&=2\ln\bigg(\frac{\beta T_M}{\sqrt{2}}\bigg), \end{align*} $$

then this is (using also equation (2.17))

(2.21) $$ \begin{align}\phi\sim 2\ln 2+\sqrt{2}\eta\quad \text{as } \eta\to -\infty.\end{align} $$

The solution of equation (2.20) with the matching condition (2.21) is

(2.22) $$ \begin{align}\phi=-2\ln \cosh\frac{\eta}{\sqrt{2}},\end{align} $$

and on the other side of the (sharp) peak, $\phi \sim 2\ln 2-\sqrt {2}\eta $ as $\eta \to \infty $ .

2.3.3 The fall

To the right of the peak, the exponential is again negligible and if the next minimum is at $T=T_m$ , then

$$ \begin{align*} \theta\approx -\theta_m+\tfrac{1}{2}\beta(T-T_m)^2, \end{align*} $$

and adopting equation (2.19), this matches to the peak if $T_m=2T_M$ , as we might expect. The period of the oscillation is thus $2T_M\approx 2\sqrt {2\theta _m/\beta }$ .

We have compared the various analytic approximations we have constructed here to the actual numerical solution and they appear, even at this leading order level of accuracy, to give a relatively good account of the analytic structure of the solution.

2.4 Multiple pulses

A difficulty with the approximation now becomes apparent. Although the pulse shown in Figure 4 decays smoothly back to zero, which we would anticipate as occurring through another transition region, efforts to compute this fail (for the moment), as the solution apparently continues to oscillate.

Figure 4 A close-up of the pulse at $t=287$ in Figure 3.

The clue to the resolution of this conundrum lies in Figure 5, which shows part of a solution at the much smaller value $\varepsilon =0.001$ . The figure plots $\theta $ and a scaled value of $x-1$ (or X). It can be seen that at $t=30$ , $\theta \approx -15$ so that $z\approx 0.3\times 10^{-6}\sim \varepsilon ^2b$ and the solution is in the slow phase. As $x\to 1$ ( $5(x-1)\to 0$ in the figure), $\theta $ increases and there is a transition to a pulse at approximately $t=30.7$ . And not just one pulse, but a sequence of four, with a late aborted one, before the trajectory recovers back to the slow phase. These pulses manifest themselves as “bursts” in the time series for z.

Figure 5 A series of bursts. Shown are $\theta =\ln z$ and $5(x-1)$ . The parameter values used are $a=0.2$ , $b=0.2$ , $c=1,\!000$ .

The clue consists of the observation that the slope of X clearly decays between pulses, and this causes a gradual waning of the amplitude of the oscillations until, when $X<0$ , no further oscillations take place. Inspection of the graph of $\dot {x}$ suggests that this change in slope is just part of the secular change due to the evolution of X. In more detail, equation (2.11) can be written in the form

(2.23) $$ \begin{align} \frac{d {}}{d {T}}(\dot{X}+z)=\sqrt{\varepsilon}\{a(\dot{X}+z)-1\}+O(\varepsilon). \end{align} $$

At leading order, $\dot {X}+z$ is constant, and in view of the matching conditions in equation (2.13) which imply $\dot {X}\to \beta $ and (using equation (2.14)) $z\sim \varepsilon ^{3/2}$ as $T\to 0$ , this constant is $\beta $ . Substituting this into the right-hand side of equation (2.23) allows an improved approximation in the form

(2.24) $$ \begin{align} \dot{X}+z=\beta-2\sqrt{\varepsilon}\gamma T+\cdots{}, \end{align} $$

bearing in mind the definition of $\gamma $ in equation (2.7), and therefore a correction to equation (2.16) is

(2.25) $$ \begin{align} \ddot{\theta}+e^\theta=\bar{\beta}=\beta-2\sqrt{\varepsilon}\gamma T, \end{align} $$

again subject to the initial conditions in equation (2.15).

For sufficiently small $\varepsilon $ , equation (2.24) suggests that there will be a series of pulses ¹ , which terminate approximately when $\bar {\beta }$ in equation (2.25) becomes negative. The reason then that Figure 4 exhibits a single pulse is because the value of $\beta $ is also formally small, or the value of $\sqrt {\varepsilon }$ is not so small. We will explore this in the following section.

2.5 The last pulse

Formally, for sufficiently small $\varepsilon $ , we would treat the series of pulses as a slowly varying nonlinear oscillator, and use the method of Kuzmak and Luke (see [8]). The final pulse then occurs when $\bar {\beta }$ approaches zero. However, for the more commonly viewed lower values of c wherein a single pulse occurs, it is of interest to address that case directly. Much of our discussion below is closely related to the work of Cheng and Chen [2], except that the approximations we make are based on explicit asymptotic limits.

A single pulse may occur if $\beta $ is small. This occurs if a is small. Specifically, if we expand equation (2.4) for small a (and hence suppose ), we find that, starting with $\dot {x}=0$ at $t=0$ , then $x=1$ at $t\approx \pi -\{\pi a-2(1-C)\}^{1/2}$ , and consequently

(2.26) $$ \begin{align} \beta\approx \{\pi a-2(1-C)\}^{1/2}. \end{align} $$

A puritan approach would then start again with the additional assumption that $a\ll 1$ , but in the interests of brevity, we will go off-piste and anticipate that equation (2.25) applies, but we will suppose that $\beta $ as well as $\varepsilon $ is small. With both a and $\beta $ being small, we can take $\gamma =1/2$ , and we presume the initial conditions in equation (2.15) still apply. (Note that properly speaking, the solution would require a multiple scales expansion.)

An appropriate distinguished limit is to assume that $\beta \sim \varepsilon ^{1/3}$ and, to this end, we define the quantity B by the relation

(2.27) $$ \begin{align} \beta=\varepsilon^{1/3}B, \end{align} $$

and we will assume that $B=O(1)$ . We then define

(2.28) $$ \begin{align} T=\frac{\tau}{\varepsilon^{1/6}},\quad\theta=\frac{1}{3}\ln\varepsilon+\Theta, \end{align} $$

so that equation (2.25) becomes

(2.29) $$ \begin{align} {\Theta}"+e^\Theta=B-\tau, \end{align} $$

where the dashes denote differentiation with respect to $\tau $ and the initial conditions are

$$ \begin{align*} \Theta=-\Theta_m,\quad {\Theta}'=0\quad\text{at } \tau=0, \end{align*} $$

where

(2.30) $$ \begin{align} \Theta_m=\ln\bigg[\frac{1}{\varepsilon b}\sqrt{\frac{B}{2\pi}}\bigg]. \end{align} $$

To solve this, we use the same asymptotic procedure as before, since $\Theta _m$ is large. Initially, we can neglect the exponential, so that

(2.31) $$ \begin{align} \Theta\approx -\Theta_m+\tfrac{1}{2} B\tau^2-\tfrac{1}{6} \tau^3. \end{align} $$

This has a maximum at $\tau =2B$ , where

$$ \begin{align*} \Theta=-\Theta_m+\tfrac{2}{3} B^3. \end{align*} $$

Formally, if $B=O(1)$ , this solution is uniformly valid and provides the aborted pulse seen in Figure 5. If B is slightly larger,

(2.32)

there will be another sharp peak.

Let us consider this in more detail. With the definitions in equations (2.6), (2.8), (2.12) and (2.28), we can define

(2.33) $$ \begin{align} x=1+\nu^2\xi,\quad z=\nu e^\Theta,\quad t=t_c+\nu\tau, \end{align} $$

where we write

$$ \begin{align*} \nu=\varepsilon^{1/3}. \end{align*} $$

In addition, bearing in mind equations (2.26) and (2.27), it is appropriate to define

$$ \begin{align*} a=\nu^2\alpha \end{align*} $$

and to suppose that $\alpha =O(1)$ . (Actually, for the example we have used with $c=10$ and $a=0.2$ , one has $\alpha \approx 0.93$ .) The equations then take the exact form

(2.34) $$ \begin{align}{\xi}"+1+\nu^2\xi+e^\Theta{\Theta}'&=\nu^3\alpha({\xi}'+e^\Theta),\dfrac{}{}\nonumber\\ {\Theta}'&=\nu^3be^{-\Theta}+\xi,\end{align} $$

with matching conditions as $\tau \to 0$ ,

(2.35) $$ \begin{align}\xi\sim B\tau-\tfrac{1}{2}\tau^2,\quad \Theta\sim -\Theta_m+\tfrac{1}{2}B\tau^2-\tfrac{1}{6}\tau^3.\end{align} $$

At leading order, this gives $\xi ={\Theta }'$ and therefore also

(2.36) $$ \begin{align}{\Theta}"'+1+e^\Theta{\Theta}'=0,\end{align} $$

with the first integral

(2.37) $$ \begin{align}{\Theta}"+e^\Theta =B-\tau,\end{align} $$

which regains equation (2.29).

As for the earlier pulse analysis and supposing that , there is a rise phase in which equation (2.31) applies:

$$ \begin{align*} \Theta\sim-\Theta_m+\tfrac{1}{2} B\tau^2-\tfrac{1}{6}\tau^3, \end{align*} $$

followed by a peak in which, as in equation (2.19),

(2.38) $$ \begin{align}\Theta=\Theta_M+\phi,\quad \tau=\tau_M+e^{-\Theta_M/2}\eta,\end{align} $$

so that the solution is again given by equation (2.22), provided we choose

(2.39) $$ \begin{align} &\Theta_M=-\Theta_m+\tfrac{1}{2} B\tau_M^2-\tfrac{1}{6}\tau_M^3-2\ln 2,\nonumber\\ &\big(B\tau_M-\tfrac{1}{2}\tau_M^2\big)e^{-\Theta_M/2}=\sqrt{2}, \end{align} $$

with order of magnitude solutions

(2.40) $$ \begin{align}\tau_M\sim\sqrt{\frac{2\Theta_m}{B}},\quad \Theta_M\sim 2\ln\bigg(\frac{B\tau_M}{\sqrt{2}}\bigg),\end{align} $$

consistent with equation (2.32). Following the sharp peak, there is a decline in which

(2.41) $$ \begin{align}\Theta=K-L(\tau-\tau_M)+\tfrac{1}{2} (B-\tau_M)(\tau-\tau_M)^2-\tfrac{1}{6}(\tau-\tau_M)^3,\end{align} $$

and matching to the peak solution given by equation (2.38) and also equation (2.22) implies

$$ \begin{align*} K=\Theta_M+2\ln 2,\quad L=B\tau_M-\tfrac{1}{2}\tau_M^2=\sqrt{2}e^{\Theta_M/2}. \end{align*} $$

In seeking to determine the downward sloping part of the map in Figure 2, the key observation is in Figure 6, which shows a close-up of the passage of x and $\ln z$ (and thus $\Theta $ ) through a pulse when $c=10$ . Clearly, the sharp peak in $\Theta $ causes a rapid decrease in ${\Theta }'=\xi $ and thus also x. This can also be seen in Figure 5. However, whereas the drop in Figure 5 is relatively small, it can be seen in Figure 6 that it is actually of numerical $O(1)$ .

Figure 6 A detail of the solution in Figure 3 in the vicinity of a near-homoclinic burst: x (and y) reach a value close to zero before z reaches the transition to the slow equilibrium in equation (2.2).

From equation (2.39), using the fact that the jump in slope of $\phi (\eta )$ across the pulse is $-2\sqrt {2}$ , it follows that the jump in x across the pulse is

$$ \begin{align*} -\lambda=-\nu^2L\approx -2\sqrt{2}\nu^2e^{\Theta_M/2}\sim -2\nu^2B\tau_M\sim - 2\nu^2\sqrt{2B\Theta_m}, \end{align*} $$

where we have used the crude estimates in equation (2.40). Using the definition of $\Theta _m$ in equation (2.30),

$$ \begin{align*} \lambda\sim 2\varepsilon^{2/3}\bigg\{2B\ln\bigg(\frac{1}{\varepsilon b}\sqrt{\frac{B}{2\pi}}\bigg)\bigg\}^{1/2}. \end{align*} $$

If we use values $B=1$ , $b=0.2$ , then for $c=10$ , we find $\lambda \approx 1.05$ , whereas for $c=1,\!000$ , it is $0.078$ . If we want to be fastidious, we could then consider a distinguished limit in which

$$ \begin{align*} b=\frac{b^*}{\varepsilon\sqrt{2\pi}}\exp\bigg[-\frac{1}{8\varepsilon^{4/3}}\bigg],\quad b^*=O(1), \end{align*} $$

but in the present note, we will avoid such delicacy. The quotation from Montaigne’s essay “Of usefulness and honesty” (De l’utile et de l’honnête) at the head of Murray’s [13] book is wonderfully apt and provides our watchword here.

Figure 7 illustrates what happens when $c=10$ . Relatively low values of C and thus B cause low pulses, which return to the plane approximately on the y axis. Higher values of C and thus B and hence also $\lambda $ give higher peaks of z (as can be seen in Figure 1) and also x, but the subsequent drop is larger, and brings the trajectory almost back to the origin, and thus lower values of C. This will explain the descending part of the map in Figure 2.

Figure 7 A view looking directly down on the attractor of equation (1.1) for the values $a=0.18$ , $b=0.2$ , $c=10$ . Note that the coordinates are unscaled.

2.6 A Poincaré map

We can use these approximations to construct an approximate Poincaré map for the solution. We use the plane $y=0$ , $x<0$ as the Poincaré surface. From equation (1.2), this is $\dot {x}+z=0$ , and in the slow phase where $z=O(\varepsilon ^2)$ (equation (2.2)), this is approximately $\dot {x}=0$ , $x<0$ , which we used to derive equation (2.5), which gives the rising limb of the map in Figure 2.

We restrict attention to the single-pulse case described in the previous section. Working our way through the various scalings, we find from equation (2.33) (bearing in mind that $\xi \approx {\Theta }'$ ) that

$$ \begin{align*} y=-\nu(\Theta"+e^\Theta),\quad x=1+\nu^2{\Theta}'. \end{align*} $$

Since ${\Theta }"+e^\Theta =B-\tau $ in the pulse, we see that $y=0$ when $\tau =B$ , but this is when $x>0$ , as is also clear in Figure 7. The next crossing of $y=0$ when $x<0$ is back on the slow phase, and for x, this begins immediately after the sharp peak, where equation (2.41) implies

(2.42) $$ \begin{align}x\sim 1-\lambda+\beta_M(t-t_M)-\tfrac{1}{2}(t-t_M)^2,\end{align} $$

where we define

$$ \begin{align*} \lambda=\nu^2L, \end{align*} $$

as earlier, and also

$$ \begin{align*} t_M=t_c+\nu\tau_M,\quad \beta_M=\beta-\nu\tau_M. \end{align*} $$

In the slow phase, the exponential in equation (2.37) is negligible and x satisfies equation (2.3) again, and thus also equation (2.4), which we write in the form

$$ \begin{align*} x=-C'e^{at'/2}\cos(\omega t'+\chi),\quad \dot{x}=C'e^{at'/2}\sin\omega t', \end{align*} $$

where, if the next intersection with $\dot {x}=0$ , $x<0$ , is at $t=t_2>t_M$ , we define $t'=t-t_2$ , so that equation (2.42) implies

(2.43) $$ \begin{align}x=1-\lambda,\quad \dot{x}=\beta_M\quad \text{at } t'=-s,\quad s=t_2-t_M. \end{align} $$

The value of $C'$ is the next iterate of the approximate Poincaré map. We retain the $\lambda =\nu ^2L$ term in view of its significance at lower values of c, but note that equation (2.42) then appears to be inconsistent with the fact that $\ddot {x}+x\approx 0$ ; the resolution is that the quadratic term coefficient also has a correction to $(1-\nu ^2L)/2$ , as can easily be checked by including the $O(\nu ^2)$ correction in equation (2.34); but choosing to satisfy equation (2.43) automatically caters for this.

The value of $C'$ is then determined by eliminating s from the pair of equations

$$ \begin{align*} 1-\lambda=-C'e^{-as/2}\cos(\omega s-\chi),\quad \beta_M=-C'e^{-as/2}\sin\omega s. \end{align*} $$

Putting all of this together and defining $\omega s=\pi +\sigma $ , the Poincaré map can be approximated as a one-dimensional map $C\to C'$ , where

(2.44) $$ \begin{align}C'=\frac{\beta_M}{\sin\sigma}\exp\bigg[\frac{a(\pi+\sigma)}{2\omega}\bigg],\quad \sigma=\tan^{-1}\bigg[\frac{\omega\beta_M}{1-\lambda-a\beta_M/2}\bigg],\end{align} $$

where, from equations (2.30) and (2.39),

$$ \begin{align*}\beta_M=\beta-\varepsilon^{1/3}\tau_M,\quad \beta=\varepsilon^{1/3}B,\quad \omega=\big(1-\tfrac{1}{4} a^2\big)^{1/2},\end{align*} $$

(2.45) $$ \begin{align}\lambda=\varepsilon^{2/3}B^2\big(S-\tfrac{1}{2} S^2\big),\quad \tau_M=BS,\end{align} $$

S is given in terms of B by the zero contour of

(2.46) $$ \begin{align}F(S,B)=B^3\bigg(\frac{1}{2} S^2-\frac{1}{6} S^3\bigg)-2\ln\bigg[B^2\bigg(S-\frac{1}{2} S^2\bigg)\bigg]-\ln\bigg[\frac{2}{\varepsilon b}\sqrt{\frac{B}{2\pi}}\bigg],\end{align} $$

and finally B is given in terms of C by equation (2.26):

(2.47) $$ \begin{align}B=\frac{\{\pi a+2(C-1)\}^{1/2}}{\varepsilon^{1/3}}.\end{align} $$

We assume that the drop in x across the pulse, $\lambda $ , is less than $1- a\beta _M/2$ , otherwise equation (2.44) implies that $C'<0$ . In fact, what happens in this latter case is that the bottom-most (lowest y) strand of the collapsing pulse in Figure 7 lands back on the slow phase plane in $y<0$ , so that its next traverse through $y=0$ is in $x>0$ and a further half-turn is necessary to reach the Poincaré surface; this occurs for $a=0.2$ , which is why there is an upturn at the right in Figure 2 and also why Figure 7 uses the value $a=0.18$ , for which this does not occur (and the corresponding map is simply unimodal).

Given C, equation (2.47) determines B and then equation (2.46) determines the function $S(B)$ implicitly. After that, the other definitions in equations (2.44) and (2.45) provide an explicit recipe for $C'$ . Just as the end is in sight, a final complication occurs: the function $S(B)$ is multi-valued, as shown in Figure 8. It is easy to explain this, as there are two distinct asymptotic limits of $F=0$ as $B\to \infty $ or $S\to 0$ , these being

$$ \begin{align*} S\sim \frac{1}{B^{9/4}}\mbox{ (lower)},\quad S\sim\bigg(\frac{3\ln B}{B^3}\bigg)^{1/2}\mbox{ (upper)}\quad \text{as } B\to\infty. \nonumber \end{align*} $$

Figure 8 The zero contour of $F(S,B)$ as defined in equation (2.46).

Which should we choose? The answer to this lies in the dynamics of the pulse and can be seen in Figure 6, for example. It is clear that x reaches its maximum before $\Theta $ does. The latter is at $\tau =\tau _M$ , while the former is at $\tau \approx B$ (where $e^\Theta $ can be ignored). Therefore, $\tau _M>B$ , or equivalently (see equation (2.45)), $S>1$ . It is perhaps significant that the graph in Figure 8 has its sole (that is, the left-hand limit of the curve, or the sole of the foot of the sleeping giant) at $S\approx 1$ . Therefore, we choose the upper branch of the figure to provide a suitable uniquely defined function $S(B)$ . It should be noted, nevertheless, that $S<1$ on the upper branch for , so the reliability of this argument may require further attention. However, corresponds to , which is where the approximate map starts to break down anyway. It is tempting to suppose that the lower branch corresponds to the minimum in x after the pulse maximum.

Figure 9 shows the map we obtain from the combination of the approximations in equation (2.5) for $C<1-\pi a/2\approx 0.686$ and equation (2.44) for C values larger than this. In fact, we have simply plotted the minimum of the two expressions, where the join occurs at approximately 0.715. The approximation in equation (2.44) assumes a sharp peak of $\Theta $ and it is reasonable to suppose that the region near the join is in fact smooth, accommodating the weak pulses in that part of the graph. The approximate map also breaks down near $C=1.6$ , since the increasing value of $\lambda $ causes $\sigma $ to rapidly increase towards $\pi /2$ , causing the upturn in $C'$ . In fact, Figure 7 illustrates the issue. At these large values of $\lambda $ , the end-pulse trajectory lands very close to the origin and thus also the fixed point of the equations. Our approximate map is actually mimicking the upturn in Figure 2, but it is unable to deal with the homoclinic approach to the fixed point, which in these kinds of approximation requires separate approximation [7]).

Figure 9 One-dimensional approximate map for the rescaled Rössler equations, at $a=0.2$ , $b=0.2$ , $c=10$ . The upper branch of the curve in Figure 8 has been approximated by a function which fits well until the sole in Figure 8 is approached. Note in comparing this with Figure 2 that the present figure is in scaled units, whereas Figure 2 is in the original variables, a factor of ten larger.

3 Conclusions

We have sketched a way in which the Rössler equations can be asymptotically solved in the limit $c\gg 1$ . When dealing with chaotic systems, it is natural to question how approximation methods can give useful information when the solutions display sensitive dependence on initial conditions. Although we will not attempt to elaborate the argument here, the resolution of this conundrum is similar to that involved in the use of Shilnikov’s method to construct strange attractors (see, for example, [6, Ch. 4]), which is, in fact, a procedure based on matched asymptotic expansions as we use here. The key point is that one constructs an approximate Poincaré map over a finite time interval: precisely the context in which asymptotic expansions “converge”. The two key ingredients of chaos – the existence of an infinite number of periodic orbits and the existence of a Smale horseshoe – are guaranteed by the smoothness of the approximation and the implicit function theorem. The fact that an individual trajectory will deviate from its approximation over long times does not alter this essential fact.

The basic approximation is relatively straightforward and leads to the prediction of the observed structure at values $c\sim 10^3$ in which more or less periodic sequences of bursts occur. The asymptotic limit is manifested by the almost two-dimensional structure of the attractor and this feature is maintained at lower values of $c\sim 10$ . We have taken the risky step of endeavouring to explain the chaotic solutions which occur at such values. The investigation is clouded with a number of difficulties, not the least of which is the blurring of the distinction between the sizes of $\varepsilon ^{1/3}$ and $\ln (1/\varepsilon )$ compared to one. Nevertheless, we are able to produce a credible approximation to the unimodal map which can occur.

That said, the study is one that only dips an asymptotic toe into the analysis of these equations and their further investigation warrants a much more thorough treatment. In the same way that Sparrow sorted out the Lorenz equations in his thesis [16] and published the results as a book [17], we hope that our exploration of the Rössler equations will inspire a future examination of these equations in much greater depth than we have attempted here.