Empirical generating partitions of driven oscillators using optimized symbolic shadowing

We present an optimized version of the symbolic shadowing algorithm for coarse graining a continuous state dynamical system, originally due to Hirata and co-workers [Phys. Rev. E 70, 016215 (2004)]. We validate our algorithm by finding generating partitions presented previously in the literature. We show that, unlike the original, the optimized algorithm can approximate generating partitions for periodically driven continuous-time nonlinear oscillators. We recover known generating partitions for the driven Duffing oscillator and compute generating partitions for the driven van der Pol oscillator. We also examine the problem of how algorithms such as ours can be applied “objectively,” that is, by starting from arbitrary initial partition guesses. By applying our algorithm to large ensembles of initial random partitions, we show that symbolic shadowing leads to a multiplicity of candidate generating partitions that localize points in phase space to a high degree, thus making it difficult to select the best choice(s). We thus propose using the Lempel-Ziv complexity to identify partitions from this set of candidates that are, in a specific sense, “minimal,” i.e., those with contiguous cells, fewer cell boundaries, and a smaller number of cells compared to their rivals. We also show how our methods can be used to indicate the appropriate number of symbols needed to approximate a generating partition.

Figure 1

Generating partition consisting of three cells labeled using symbols $\{0,1,2\}$ for the strange attractor of the stochastically unperturbed, two-well, driven Duffing oscillator given by a second-order ordinary differential equation (ODE): $\ddot{x}+0.25\dot{x}-x+x^3=0.4cost$. Also shown in the phase space ${\mathbb{R}}^2$ are iterates of the Poincaré map, each obtained by numerical integration for half the forcing period (i.e., from 0 to $\pi$), followed by inversion of coordinates: $x\to -x$ and $\dot{x}\to -\dot{x}$. Although the partition has been obtained for the strange attractor, the cell boundaries have been extrapolated outside the attractor for clarity. This partition, obtained from time series data using methods discussed later in the paper, matches well with the previously published generating partition found using independent methods [12, 19].

Figure 2

Comparison of cylinder diameter ${\mathcal{D}}_{\text{max}}^l$ [Eq. (6)] vs string length $l$ for four partitions (with different alphabet sizes $k$) of the Duffing attractor, empirically estimated from a time series of four million points. Note the base-10 logarithmic scale on the $y$ axis. For the $k=3$ (three symbol) generating partition of Fig. 1, ${\mathcal{D}}_{\text{max}}^l$ shows a sustained decrease until about $l=10$. For the $k=16$ “maximum entropy” (equiprobable) and $k=21$ equisized (same area) partitions, ${\mathcal{D}}_{\text{max}}^l$ decreases by a factor of one-half until $l=3$, without further substantial decrease for longer strings. For a $k=3$ random coloring every data point is uniformly randomly assigned a symbol from $\mathcal{A}=\{0,1,2\}$. Symbol sequences from this partition correspond to sets of points spread over all of the attractor so that ${\mathcal{D}}_{\text{max}}^l$ does not decrease with $l$ and equals the diameter of the attractor itself. As expected, unlike other partitions, the generating partition consistently increases the degree of localization of points on the Duffing attractor, even at string lengths larger than $l=5$.

Figure 3

Variation of maximum cylinder diameter ${\mathcal{D}}_{\text{max}}^{l,m}$ [Eq. (5)] for three-symbol generating partition of the Duffing attractor (Fig. 1). For each $l$, the $l$ largest diameters ${\mathcal{D}}_{\text{max}}^{l,m}$ are shown by $\times$ and solid line, corresponding to values of $m$ from 0 to $l-1$ from left to right within each curve. We see ${\mathcal{D}}_{\text{max}}^{l,m}$ varies significantly with $m$, with minimum ${\mathcal{D}}_{\text{max}}^l$ shown by $\square$ and dashed line: these minima are the same values as in Fig. 2, with $m_{\text{opt}}$ as in Eq. (7).

Figure 4

The most frequently obtained three-cell partition of the Duffing attractor, generated by application of the symbolic shadowing algorithm to thousands of random initial symbol sequences. Note that this does not match the generating partition of Fig. 1, a limitation addressed by our modification to the original algorithm. The cells of the partition are a union of multiple, apparently disjoint regions in the phase space. Inset: An enlarged view of one area in the phase space where the cells for symbols {0} and {1} are intermingled in a complicated way.

Figure 5

Cylinder sets $\mathcal{C}\left(s_{-m}^{l-m-1}\right)$ for $m=0,...,9$ and $l=10$, corresponding to the symbolic itinerary $\{0,2,1,1,0,2,1,1,0,2\}$, the most probable itinerary (probability $\approx 0.0209$) among all observed length-10 strings obtained from the generating partition (Fig. 1), overlaid on the Duffing attractor. Sets were estimated from a time series of four million points. As the number of backward symbols $m$ change from 0 to 9, we get cylinders (heavy dots) of different shapes and sizes due to stretching and folding on the Duffing attractor. Also indicated, by the symbol $*$, is each cylinder's representative $\mathit{r}\left(s_{-m}^{l-m-1}\right)$ [Eq. (9)]. Observe that the cylinder size is smallest, and is best approximated by its representative, when $m=2$.

Figure 6

Comparison of shadowing errors (i.e., mean and maximum squared errors) and the maximum cylinder diameter vs number of backward symbols $m$ in delimited strings of length $l=10$ for the generating partition (Fig. 1) of the Duffing attractor. The solid lines show values estimated from all of the cylinder sets over the attractor, while the dotted lines correspond to the cylinder sets (Fig. 5) of the symbolic itinerary $\{0,2,1,1,0,2,1,1,0,2\}$. Both shadowing errors as well as the cylinder diameter are minimized at $m=2$ suggesting an optimal choice for the formation of cylinder sets.

Figure 7

Typical convergence traces of the maximum squared error for symbolic shadowing (SS) and optimized symbolic shadowing (OSS) algorithms starting from 175 initial symbol sequences generated using the RS method. For all these sequences, the optimized algorithm converged to the generating partition of Fig. 1 for which mean and maximum squared errors were $5.0012\times {10}^{-4}$ and $1.6851\times {10}^{-2}$ respectively, while the original algorithm converged to the partition of Fig. 4 with mean and maximum squared errors equal to $7.9951\times {10}^{-4}$ and $7.2883\times {10}^{-2}$ respectively. Also, the original algorithm took about thrice as many iterations as required by the optimized algorithm before it stopped.

Figure 8

Ratio of average execution time per iteration for the optimized symbolic shadowing algorithm (OSS) and the original algorithm (SS) with increasing $l$, estimated using convergence traces of Fig. 7. In the OSS algorithm, although $l$ symbol sequences are found at each iteration, it does not run $l$ times slower as compared to the SS algorithm. For the Duffing oscillator, at $6<l\le 10$, we have $T_{\mathrm{iter}}^{\text{OSS}}/T_{\mathrm{iter}}^{\text{SS}}\approx \epsilon l$, where $\epsilon \in (0.42,0.53)$. A simple extrapolation, assuming $\epsilon$ increases linearly with $l$, leads to a prediction: ${\epsilon }_{\text{max}}\approx 0.61$ at $l=15$.

Figure 9

Typical variation of values of $m_{\text{opt}}$, estimated at each iteration, and $l$ for one of the traces in Fig. 7 for the optimized algorithm (OSS). For reference, values of $m=\lfloor l/2\rfloor$ for the original algorithm (SS) are also shown. As can be seen, the estimated $m_{\text{opt}}$ is less than or equal to $\lfloor l/2\rfloor$, dependent on $l$, and seldom changes from iteration to iteration at fixed $l$.

Figure 10

Left: Two-cell generating partitions of the Hénon map (a), (b), Ikeda map (c), and the driven, two-well Duffing oscillator at an alternate set of parameters (d), respectively, found with optimized symbolic shadowing (cell boundaries shown with solid lines). For Heńon and Ikeda maps, cell boundaries for partitions obtained from the original algorithm (broken lines) lie almost exactly on top of those found using the optimized algorithm. As a result, the cell boundaries of these maps obtained using both the algorithms are virtually indistinguishable from each other in these figures. Right: Maximum cylinder diameter ${\mathcal{D}}_{\text{max}}^{l,m}$ normalized by attractor diameter ${\mathcal{D}}^0$ vs number of backward symbols $m$ in delimited strings of length $l=15$, compared for all systems on left. The sampling of strings was done using $N={10}^6$ iterates on the attractor for each system. The vertical broken line corresponds to $m=\lfloor l/2\rfloor =7$. For the Hénon and Ikeda maps, there is only about 2–4% of the attractor diameter increase in ${\mathcal{D}}_{\text{max}}^{l,m}$ for $m$ values between $m=m_{\text{opt}}$ and $m=\lfloor l/2\rfloor$, as opposed to an 18% increase for the Duffing oscillator.

Figure 11

Six examples of partitions of the Duffing attractor, with two cells [partitions (a)–(c)] and three cells [partitions (d)–(f)]. Symbols $\{0,1,2\}$ label partition cells. Cell boundaries for partitions (a), (b) and (d), (e) were chosen so that each cell for a given partition had an equal number of data points falling in it. Furthermore, by construction, the boundaries of partition (c) are the union of those of partitions (a) and (b). Similarly, cell boundaries of (f) consist of boundaries of (d) and (e). Partitions (c) and (f) have noncontiguous cells which are disjoint unions of smaller boxes.

Figure 12

Comparison of Lempel-Ziv (LZ) complexity $c_E$ of partitions of Fig. 11. A time series of four million points of the Duffing oscillator was divided into 20 consecutive segments of length $N=2\times {10}^5$ each. For each partition, these segments were encoded and their $c_E$ values are as shown in the bar plot. Statistical fluctuations in $c_E$ estimates are about 1–$2\%$. Evidently, partitions (c) and (f), which have noncontiguous cells but use same number of symbols as partitions (a), (b) and (d), (e) respectively, are more complex in the LZ sense. Furthermore, three-cell partitions (d)–(f) are also more complex than two-cell partitions (a)–(c).

Figure 13

For Duffing oscillator with parameters as for Figs. 1 and 4: comparison of three-cell partitions obtained from symbolic shadowing, presented in normalized ${\mathcal{D}}_{\text{max}}^l$ vs $c_E$ plane. All distinct solutions from four data sets obtained after applying both original (SS) and optimized symbolic shadowing (OSS). Initial symbol sequences generated using both RR and RS methods. Heavy lines indicate piecewise linear greatest lower envelope for each of the four sets of solutions. Up arrow and left arrow point to most frequently appearing partitions, as in Figs. 1 and 4, respectively.

Figure 14

For Duffing oscillator with parameters as in Fig. 13: comparison of partitions in the normalized $c_E-{\mathcal{D}}_{\text{max}}^l$ plane, obtained with optimized symbolic shadowing applied to initial symbol sequences prepared using the RR method. All distinct solutions from three sets corresponding to alphabet sizes $k\in \{2,3,4\}$ of the initial sequences are shown. The up arrow points to the generating partition of Fig. 1.

Figure 15

van der Pol Oscillator: comparison of partitions obtained from symbolic shadowing algorithms in the normalized $c_E-{\mathcal{D}}_{\text{max}}^l$ plane. All distinct solutions from set 1 to 4 (see Table 1) are shown. All are three-cell partitions except for the encircled ones which have two cells. A time series of two million points of the van der Pol oscillator was divided into ten consecutive segments of length $N=2\times {10}^5$ each. After encoding each of these segments using a given partition, values of $c_E$ and ${\mathcal{D}}_{\text{max}}^l$ (at $l=20$) were calculated. The points shown correspond to minimum of $c_E$ and maximum of ${\mathcal{D}}_{\text{max}}^l$ over these ten values for each partition. The cylinder diameter ${\mathcal{D}}_{\text{max}}^l$ was normalized using the diameter of the van der Pol attractor, ${\mathcal{D}}^0=1.6026$, while the LZ complexity $c_E$ was normalized using its asymptotic value for binary random sequences, i.e., $N/{log}_{2}N\approx 11357$ phrases. Heavy lines constitute piecewise linear, nonincreasing, greatest lower envelope for each of the four sets of solutions. Up arrows, with increasing values of LZ complexity, point to partitions (a) and (b) respectively of Fig. 17, while the down arrow points to the partition (i) of Fig. 18.

Figure 16

van der Pol Oscillator: comparison of partitions in the normalized $c_E-{\mathcal{D}}_{\text{max}}^l$ plane (similar to Fig. 15) obtained after applying the optimized symbolic shadowing algorithm. All distinct solutions from sets 3, 5, and 6 (see Table 1) are shown. See caption of Fig. 15 for details on the calculation of $c_E$ and ${\mathcal{D}}_{\text{max}}^l$ values.

Figure 17

Partitions of the van der Pol attractor obtained using optimized symbolic shadowing algorithm consisting of three contiguous cells labeled $\{0,1,2\}$.

Figure 18

Three-cell partitions of the van der Pol attractor found using optimized symbolic shadowing algorithm. Unlike partitions of Fig. 17, these solutions are sensitive to parameters of the numerical integration algorithm used to prepare the data set and also to some extent the initial string length $l_{\text{initial}}$. Partitions (i) and (ii) are very similar except in a small area around the boundary between cells {0} and {1}. An enlarged view of this area for partition (i) is shown in its inset.