Wind Turbines - Advances and Challenges in Design, Manufacture and Operation

*Edited by Karam Maalawi*

Published in London, United Kingdom

Wind Turbines - Advances and Challenges in Design, Manufacture and Operation http://dx.doi.org/10.5772/intechopen.97967 Edited by Karam Maalawi

#### Contributors

Lijun Cai, Karthik Jayanarasimhan, Vignesh Subramani-Mahalakshmi, Mohamed Mahran Kasem, Tefera Kitaba, Karam Youssef Maalawi, Vitaly Beresnevich, Janis Viba, Martins Irbe, Anubhav Jain, Enas Raafat Maamoun Shouman, Girmaw Teshager Bitew, Minxiao Han, Laurent Barthélemy, Alice Elizabeth Gonzalez, José Cataldo-Ottieri, John Campbell, Balaraman Kannan, Bastin Jeyaraj, Jigneshkumar P. Pramodbhai Desai, Ramesh Kumar Behara, Kavita Behara

#### © The Editor(s) and the Author(s) 2022

The rights of the editor(s) and the author(s) have been asserted in accordance with the Copyright, Designs and Patents Act 1988. All rights to the book as a whole are reserved by INTECHOPEN LIMITED. The book as a whole (compilation) cannot be reproduced, distributed or used for commercial or non-commercial purposes without INTECHOPEN LIMITED's written permission. Enquiries concerning the use of the book should be directed to INTECHOPEN LIMITED rights and permissions department (permissions@intechopen.com).

Violations are liable to prosecution under the governing Copyright Law.

Individual chapters of this publication are distributed under the terms of the Creative Commons Attribution 3.0 Unported License which permits commercial use, distribution and reproduction of the individual chapters, provided the original author(s) and source publication are appropriately acknowledged. If so indicated, certain images may not be included under the Creative Commons license. In such cases users will need to obtain permission from the license holder to reproduce the material. More details and guidelines concerning content reuse and adaptation can be found at http://www.intechopen.com/copyright-policy.html.

#### Notice

Statements and opinions expressed in the chapters are these of the individual contributors and not necessarily those of the editors or publisher. No responsibility is accepted for the accuracy of information contained in the published chapters. The publisher assumes no responsibility for any damage or injury to persons or property arising out of the use of any materials, instructions, methods or ideas contained in the book.

First published in London, United Kingdom, 2022 by IntechOpen IntechOpen is the global imprint of INTECHOPEN LIMITED, registered in England and Wales, registration number: 11086078, 5 Princes Gate Court, London, SW7 2QJ, United Kingdom

British Library Cataloguing-in-Publication Data A catalogue record for this book is available from the British Library

Additional hard and PDF copies can be obtained from<EMAIL_ADDRESS>
Wind Turbines - Advances and Challenges in Design, Manufacture and Operation Edited by Karam Maalawi p. cm. Print ISBN 978-1-80355-492-1 Online ISBN 978-1-80355-493-8 eBook (PDF) ISBN 978-1-80355-494-5


*Article* **Alternate Entropy Computations by Applying Recurrence Matrix Masking**

**Charles L. Webber, Jr.**

Department of Cell and Molecular Physiology, Health Sciences Campus, Loyola University Chicago, 2160 South First Avenue, Maywood, IL 60153, USA<EMAIL_ADDRESS>Tel.: +1-708-638-7497

**Abstract:** In practicality, recurrence analyses of dynamical systems can only process short sections of signals that may be infinitely long. By necessity, the recurrence plot and its quantifications are constrained within a truncated triangle that clips the signals at its borders. Recurrence variables defined within these confining borders can be influenced more or less by truncation effects depending upon the system under evaluation. In this study, the question being asked is what if the boundary borders were tilted, what would be the effect on all recurrence variables? This question was prompted by the observation that line entropy values are maximized for highly periodic systems in which the infinitely long line elements are truncated to different unique lengths. However, by redefining the recurrence plot area to a 45-degree tilted box within the triangular area, the diagonal lines would consequently be truncated to identical lengths. Such masking would minimize the line entropy to 0.000 bits/bin. However, what new truncation influences would be imposed on the other recurrence variables? This question is examined by comparing recurrence variables computed with the triangular recurrence area versus boxed recurrence area. Examples include the logistic equation (mathematical series), the Dow Jones Industrial Average over a decade (real-word data), and a square wave pulse (toy series). Good agreement among the variables in terms of timing and amplitude was found for most, but not all variables. These important results are discussed.

**Keywords:** nonlinear dynamics; recurrence quantifications; line entropy; recurrence matrix masking

## **1. Introduction**

As an extension of recurrence plots [1], recurrence quantification analysis (RQA) was introduced by Zbilut and Webber [2,3] almost three decades ago and Marwan et al. [4] two decades ago. These quantifications include eight recurrence variables extracted from recurrence plots [5] which have proven to have utility in general-purpose data analyses for linear and nonlinear systems alike [6]. The most challenging aspect of recurrence analyses is the setting of the multiple recurrence parameters [7]. What may be less appreciated is the effect border truncations have on all recurrence variables.

The fundamental feature of recurrence plots, which distinguishes deterministic from stochastic systems, is the presence of diagonal line structures. Since the traditional recurrence plot is symmetrical on either side of the line of identity (LOI), only the uppertriangular half of the plot is utilized. The ubiquitous LOI is excluded. Three recurrence quantification variables are directly derived from these diagonals. First, percent determinism (DETERM) is defined as the ratio of points forming diagonal line structures to the total number of recurrent points. Second, the diagonal maximum (DMAX) is defined as the integer number of points constructing the longest diagonal line. Third, information or line entropy (ENT) is defined as the Shannon entropy [8] of the histogram distribution of all diagonal line lengths within the triangular recurrence window. Because long diagonal lines are necessarily truncated at the borders of the recurrence plot, the question arises how much the truncation effect influences the accuracy of DETERM, DMAX and ENT computations, not to mention the remaining recurrence variables. This is the same question formerly asked by Kraemer and Marwan [9].

**Citation:** Webber, C.L., Jr. Alternate Entropy Computations by Applying Recurrence Matrix Masking. *Entropy* **2022**, *24*, 16. https://doi.org/ 10.3390/e24010016

Academic Editors: Franco Orsucci and Wolfgang Tschacher

Received: 31 August 2021 Accepted: 20 December 2021 Published: 23 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

This paper will focus on the effect recurrence borders have on the computation of line structures and other features of the recurrence plot. Two types of borders are examined: traditional triangular recurrence borders and novel tilted box recurrence borders. The latter represents a smaller area masked off within the larger triangular area. System studied include a mathematical system, a real-world financial system, and a contrived toy system to clearly illustrate time and amplitude shifts in all compared variables.

#### **2. Disparity between Recurrence Line Entropy Values and Lyapunov Exponents**

When this author [10] was studying the logistic equation by recurrence analysis, it was noted that when the equation was in period 1, period 2, period 4, etc., periodic modes, the line entropy values were maximized when the Lyapunov exponents were low. This counterintuitive relationship is easily explained by the differing lengths of diagonal lines which were truncated by the triangular borders into unique lengths. Censi et al. [11] also realized this effect and introduced a correction by assigning all diagonals to the same length as the central LOI. This method effectively minimized the line entropy values but carried with it the assumption that all lines are of equal length. Eroglu et al. [12] took a different approach using weighted recurrence plots. Entropic measures were redefined according to the distribution of return times, not line lengths. Indeed, Shannon entropy values became positively correlated with the largest Lyapunov exponent. More sophisticated approaches for quantifying recurrence entropy are based on matrix microstates [13] and categorical time series [14]. Finally, Kraemer and Marwan [9] described another method similar to the one introduced in this paper. That is they recomputed entropy by masking the recurrence plot such that the square recurrence matrix was windowed within a diamond pattern overlaying the traditional plot. Thus, for fully periodic systems, all diagonal lines were truncated to identical lengths, minimizing line entropy values to 0.000 bits/bin.

Because the auto-recurrence plot is symmetrical, in this paper, only the upper triangle was masked with the best fitting square (not half-diamond or rectangle as in the Kraemer-Marwan case). The box-masking will be fully described. The masking gives low entropy values for periodic systems and high entropy values for chaotic systems. It was necessary to compare these boxed entropy values against the traditional triangular entropy values for several different cases and situations. Indeed, performance testing was accomplished for all eight recurrence variables. The necessary comparisons include amplitude and timing characteristics of the recurrence variables for both triangular and boxed recurrences.

#### **3. Redefinition of Recurrence Plot Boundary Conditions**

The boxed recurrence masking can be best be illustrated graphically. Figure 1 presents the 400-point recurrence plot of a simple repeating sinewave consisting of 16 cycles or waves as computed by traditional recurrence programs RQD and RQC, and new boxed programs RQDB and RQDC (see Appendix A). Because all 16 sinewaves are identical across the time series (green), 15 border-to-border lines excluding the LOI (black) are inscribed within the triangular recurrence plot (blue and red). Each one of these deterministic lines is truncated to a unique (different) length spanning from the west border to the north border. That is, the further the line is removed from the LOI, the shorter that line is. This has ramifications for line entropy calculations as described below.

Another approach to this sinewave truncation issue is to fashion a tilted box within the recurrence triangle parallel to the LOI (45-degree, clockwise tilt) as is illustrated by the black square (Figure 1). All recurrent points within this smaller boxed area (red) are simply shorter parts of the longer lines within the larger triangular area (blue). The advantage of this system is that the 10 diagonal line structures within the box are each truncated to identical lengths.

**Figure 1.** Recurrence plot of a sinewave with no added noise. The 16 noise-free sinewave cycles represent a 100 Hz sinewave digitized at 2500 Hz (25 points/cycle). The boxed recurrences (red lines) are constrained within the tilted box. The triangular recurrences (blue lines + red lines) are constrained with the upper half of the recurrence plot (above the LOI). Parameter settings: DELAY = 1; EMBED = 2; NORM = Euclid; WINDOW = 400; RESCALE = max distance; RADIUS = 1%; LINE = 2. **Figure 1.** Recurrence plot of a sinewave with no added noise. The 16 noise-free sinewave cycles represent a 100 Hz sinewave digitized at 2500 Hz (25 points/cycle). The boxed recurrences (red lines) are constrained within the tilted box. The triangular recurrences (blue lines + red lines) are constrained with the upper half of the recurrence plot (above the LOI). Parameter settings: DELAY = 1; EMBED = 2; NORM = Euclid; WINDOW = 400; RESCALE = max distance; RADIUS = 1%; LINE = 2.

Another approach to this sinewave truncation issue is to fashion a tilted box within the recurrence triangle parallel to the LOI (45-degree, clockwise tilt) as is illustrated by the black square (Figure 1). All recurrent points within this smaller boxed area (red) are simply shorter parts of the longer lines within the larger triangular area (blue). The advantage of this system is that the 10 diagonal line structures within the box are each truncated to identical lengths. Then, when 5% random noise is added to the pure sinewave, the long diagonal line structures are chopped up into shorter segments as shown in Figure 2. In this case, the noise mitigates the influence of the truncation effects of the triangular borders insofar that the line structures within the tilted box and triangle appear very similar. Are they in reality? By the way, the embedding dimension was set to 2 for both the pure sinewave with a radius of 1% of the maximum distance. Using these parameter settings there was no false thickening

of the diagonal lines as discussed by Thiel et al. [15]. That is, the selected radius was high enough to include the noisy points, but not so high as to thicken the lines. thickening of the diagonal lines as discussed by Thiel et al. [15]. That is, the selected radius was high enough to include the noisy points, but not so high as to thicken the lines.

Then, when 5% random noise is added to the pure sinewave, the long diagonal line structures are chopped up into shorter segments as shown in Figure 2. In this case, the noise mitigates the influence of the truncation effects of the triangular borders insofar that the line structures within the tilted box and triangle appear very similar. Are they in reality? By the way, the embedding dimension was set to 2 for both the pure sinewave with a radius of 1% of the maximum distance. Using these parameter settings there was no false

*Entropy* **2022**, *24*, x FOR PEER REVIEW 4 of 18

**Figure 2.** Recurrence plot of a sinewave with 5% added noise. The 16 noise-free sinewave cycles represent a 100 Hz sinewave digitized at 2500 Hz (25 points/cycle). The boxed recurrences (red lines) are constrained within the tilted box. The triangular recurrences (blue lines + red lines) are **Figure 2.** Recurrence plot of a sinewave with 5% added noise. The 16 noise-free sinewave cycles represent a 100 Hz sinewave digitized at 2500 Hz (25 points/cycle). The boxed recurrences (red lines) are constrained within the tilted box. The triangular recurrences (blue lines + red lines) are constrained with the upper half of the recurrence plot (above the LOI). DELAY = 1; EMBED = 2; NORM = Euclid; WINDOW = 400; RESCALE = max distance; RADIUS = 1%; LINE = 2.

#### **4. Tilted Box Boundaries and Line Entropy Modifications 4. Tilted Box Boundaries and Line Entropy Modifications**

*Entropy* **2022**, *24*, x FOR PEER REVIEW 5 of 18

Figure 3 is key to understanding the differences of the boxed entropies over the triangular entropies computed from the histogram distribution of diagonal lines within the recurrence plot. The upper histogram is the distribution of diagonal lines within the noise-free sinewave (Figure 1) and the lower histogram is the distribution of diagonal lines within the sinewave with 5% added noise (Figure 2). Note that both horizontal and vertical scales are base-10 logarithmic. In the first case (blue lines), all 15 line lengths are of differing lengths (375, 350, 325, 300, . . . , 100, 75, 50 25 points) and there is but one unique line length of each. However, for the second case (red line), the 10 diagonal lines are all of identical length (133 points) and fall into a single histogram bin. Consequently, the triangular entropy computes as 3.907 bits/bin, whereas the box entropy computes as 0.000 bits/bin (see below). Figure 3 is key to understanding the differences of the boxed entropies over the triangular entropies computed from the histogram distribution of diagonal lines within the recurrence plot. The upper histogram is the distribution of diagonal lines within the noisefree sinewave (Figure 1) and the lower histogram is the distribution of diagonal lines within the sinewave with 5% added noise (Figure 2). Note that both horizontal and vertical scales are base-10 logarithmic. In the first case (blue lines), all 15 line lengths are of differing lengths (375, 350, 325, 300,…, 100, 75, 50 25 points) and there is but one unique line length of each. However, for the second case (red line), the 10 diagonal lines are all of identical length (133 points) and fall into a single histogram bin. Consequently, the triangular entropy computes as 3.907 bits/bin, whereas the box entropy computes as 0.000 bits/bin (see below).

constrained with the upper half of the recurrence plot (above the LOI). DELAY = 1; EMBED = 2;

NORM = Euclid; WINDOW = 400; RESCALE = max distance; RADIUS = 1%; LINE = 2.

**Figure 3.** Histograms of diagonal line length distributions for noise-free sinewave (upper panel) and sinewave with 5% added noise (lower panel). Diagonal line lengths and their counts are shown for triangular recurrences (blue) versus boxed recurrences (red). Corresponding, color-coded entropy values are also given. Note that entropy values in the top panel are very different, but entropy values in the lower panel are very similar. **Figure 3.** Histograms of diagonal line length distributions for noise-free sinewave (upper panel) and sinewave with 5% added noise (lower panel). Diagonal line lengths and their counts are shown for triangular recurrences (blue) versus boxed recurrences (red). Corresponding, color-coded entropy values are also given. Note that entropy values in the top panel are very different, but entropy values in the lower panel are very similar.

Computation of line entropy follows the generalized formula used for Shannon information entropy [8] as shown in Equation (1). The entropy value is maximized when each non-zero bin has the same number of counts (identical probabilities) as computed by Equation (2). Conversely, the entropy value is minimized if bin counts are restricted to a single bin. In this case, the entropy falls out as 0.000 bits/bin as given by Equation (3). So from Figure 3, where 15 bins are each filled with the count of one (blue) the entropy maximized (−log2(1/15) = 3.907 bits/bin). However, for the 10 boxed lines (red), each line is 133 units long and the entropy is minimized (−log2(1/1) = 0.000 bits/bin).

$$\text{ENT}\_{\text{gen}} = -\Sigma(\text{P}\_{\text{bin}}) \log\_2(\text{P}\_{\text{bin}}) \tag{1}$$

$$\text{ENT}\_{\text{max}} = -\log\_2(1/\text{N}\_{\text{bin}}) \tag{2}$$

$$\text{ENT}\_{\text{min}} = -\log\_2(1/1) = 0.000 \text{ bits/bin} \tag{3}$$

So for noise-free sinewaves, at least, the computed values for line entropy range from maximal entropy (triangular recurrence area) to minimal entropy (tiled box recurrence area). It is the tilted box masking that makes all the difference. Both entropy values are mathematically correct, but the boxed entropy makes more sense with respect to the complexity of the signal. This is not globally true for all situations and signals.

With the jostling of the pure sinewave with 5% random noise, the long diagonal lines both within the boxed area and within the triangular area get parceled (chopped up) into shorter segments (Figure 2). Quantitatively, the distributions of these shortened line segments are very similar (Figure 3), which is why the information entropy values for both the triangular recurrences and boxed recurrences are likewise very similar (3.000 versus 2.807 bits/bin, respectively). Since noise is ubiquitous in real-world systems, possibly the two methods of entropy computations are not that much different after all. However, such a (good) conclusion has yet to be verified using experimental data.

#### **5. Boundary Details of the Tilted Recurrence Box within the Recurrence Triangle**

The edges of the recurrence box are not smooth, and the shape of the box is not perfectly square as implied by Figures 1 and 2. For example, let us construct the largest 5 × 5 tilted square within a 14 × 14 recurrence matrix, half of which is shown in the upper panel of Figure 4. Here, the edges of the tilted square box are marked with black dots. As can be seen, each side of the box is 5 units in length. However, the area of this square is not 25 square units (dark pink pixels) but must include all gap spaces as well (light pink pixels). Secondly, the length of the diagonal must all be the same unit length (5 pixels here) increasing the box area to 45 square units (Equation (4)). Thirdly, the area of the triangle includes the box area plus all other empty pixels (white pixels) excluding the LOI (green pixels) (Equation (5)). Taking the ratio of box to rectangle areas computes the percentage of the recurrence plot occupied by the tilted box (Equation (6)).

$$\text{Area}\_{\text{box}} = \text{Side}\_{\text{box}} \cdot \text{Side}\_{\text{box}} + \text{Side}\_{\text{box}} \cdot (\text{Side}\_{\text{box}} - 1) \tag{4}$$

$$\text{Area}\_{\text{triangle}} = (\text{Side}\_{\text{triangle}} \cdot \text{Side}\_{\text{triangle}} - \text{Side}\_{\text{triangle}}) / 2 \tag{5}$$

$$\% \text{Box} = 100 \cdot \text{Area}\_{\text{box}} / \text{Area}\_{\text{triangle}} \tag{6}$$

Now, if the 14 × 14 recurrence matrix size is increased to 15 × 15 square units or even to 16 × 16 square units, the shape and area of the tilted recurrence box as defined above remain the same. This is proven graphically in panels 2 and 3, respectively, in Figure 4. At the same time, however, the area of the triangles increases as more empty pixels are added. Consequently, the ratios of box area to triangle area (%Box) must necessarily decrease. For this trio of pairs (Figure 4) the ratios decrease from 49.5% to 42.9% to 37.5% as indicated. However, by incrementing the size of the recurrence matrix up by one step to 17 × 17 square units, the legal tilted-box size likewise increases from 5 × 5 to 6 × 6 (not shown).

**Figure 4.** Tilted recurrence box of 5 × 5 units (dark and light pink pixels) fits within 3 recurrence matrices of 14 × 14 units (top), 15 × 15 units (middle), and 16 × 16 units (bottom). %Box is the ratio of the box areas to the triangle areas (excluding the green LOI) which progressively decreases as the size of the recurrence matrix increases. The blue arrow designates the middle third of five input points (e.g., P5–P9) from which the mean and standard deviation are computed. The red arrow indicates the left shift in the boxed recurrence variables by 4 points (P5–P1) for alignment with the triangular recurrence variables. **Figure 4.** Tilted recurrence box of 5 × 5 units (dark and light pink pixels) fits within 3 recurrence matrices of 14 × 14 units (top), 15 × 15 units (middle), and 16 × 16 units (bottom). %Box is the ratio of the box areas to the triangle areas (excluding the green LOI) which progressively decreases as the size of the recurrence matrix increases. The blue arrow designates the middle third of five input points (e.g., P5–P9) from which the mean and standard deviation are computed. The red arrow indicates the left shift in the boxed recurrence variables by 4 points (P5–P1) for alignment with the triangular recurrence variables.

shown).

Clearly, the "square" tilted boxes can fit within 3threeincremental sizes of recurrence matrices. To compute the integer length of the box side one needs only to divide the number of points in the triangular recurrence window by 3. Thus, for our trio example: 14/3 = 4.667 (round up to 5); 15/3 = 5.000 (retain as 5); 16/3 = 5.333 (round down to 5). Now the question arises, what happens to %Box (area ratios) as the number of points in the recurrence window increases? This answer is shown graphically in Figure 5 for multiple sets of triplets. The 5 × 5 example is indicated in red. The slope of each trio decreases as the number of points increases. For example, for set 1001-1002-1003, the area ratios are 45.51%, 44.42% and 44.33%; for set 2000-2001-2002, the area ratios are 44.48%, 44.43% and 44.39%; for set 2999-3000-3001, the area ratios are 44.47%, 44.44% and 44.41%. In the limit for large set numbers, the area ratios converge on 44.44% for all trios. All this means is that the area of the new tilted box is approximately 44% of the area of the traditional recurrence triangle. Clearly, the "square" tilted boxes can fit within 3threeincremental sizes of recurrence matrices. To compute the integer length of the box side one needs only to divide the number of points in the triangular recurrence window by 3. Thus, for our trio example: 14/3 = 4.667 (round up to 5); 15/3 = 5.000 (retain as 5); 16/3 = 5.333 (round down to 5). Now the question arises, what happens to %Box (area ratios) as the number of points in the recurrence window increases? This answer is shown graphically in Figure 5 for multiple sets of triplets. The 5 × 5 example is indicated in red. The slope of each trio decreases as the number of points increases. For example, for set 1001-1002-1003, the area ratios are 45.51%, 44.42% and 44.33%; for set 2000-2001-2002, the area ratios are 44.48%, 44.43% and 44.39%; for set 2999-3000-3001, the area ratios are 44.47%, 44.44% and 44.41%. In the limit for large set numbers, the area ratios converge on 44.44% for all trios. All this means is that the area of the new tilted box is approximately 44% of the area of the traditional recurrence triangle.

%Box = 100 · Areabox/Areatriangle

Now, if the 14 × 14 recurrence matrix size is increased to 15 × 15 square units or even to 16 × 16 square units, the shape and area of the tilted recurrence box as defined above remain the same. This is proven graphically in panels 2 and 3, respectively, in Figure 4. At the same time, however, the area of the triangles increases as more empty pixels are added. Consequently, the ratios of box area to triangle area (%Box) must necessarily decrease. For this trio of pairs (Figure 4) the ratios decrease from 49.5% to 42.9% to 37.5% as indicated. However, by incrementing the size of the recurrence matrix up by one step to 17 × 17 square units, the legal tilted-box size likewise increases from 5 × 5 to 6 × 6 (not

Areabox = Sidebox · Sidebox + Sidebox · (Sidebox <sup>−</sup> 1) (4)

Areatriangle = (Sidetriangle · Sidetriangle <sup>−</sup> Sidetriangle)/2 (5)

(6)

*Entropy* **2022**, *24*, x FOR PEER REVIEW 8 of 18

**Figure 5.** Ratio of tilted-box area to recurrence triangular area (%Box) in sets of 3 for increasing number of points in triangular recurrence window. Box sizes are shown in red integers. The triplet points (14-15-16) for the 5 × 5 box size as described in the text are also shown in red. **Figure 5.** Ratio of tilted-box area to recurrence triangular area (%Box) in sets of 3 for increasing number of points in triangular recurrence window. Box sizes are shown in red integers. The triplet points (14-15-16) for the 5 × 5 box size as described in the text are also shown in red.

#### **6. Logistic Equation 6. Logistic Equation**

At least three studies have examined the logistic map in terms of entropy over a range of control parameter a values [10,12,16] (Equation (7)). Each takes a different perspective on entropy calculations with pros and cons. The present study followed the procedures of the earlier study [10] in which parameter a was incremented on each iterated cycle in steps of 0.00001 from a = 2.8 to 4.0 yielding 120,001 points. This series of points was At least three studies have examined the logistic map in terms of entropy over a range of control parameter a values [10,12,16] (Equation (7)). Each takes a different perspective on entropy calculations with pros and cons. The present study followed the procedures of the earlier study [10] in which parameter a was incremented on each iterated cycle in steps of 0.00001 from a = 2.8 to 4.0 yielding 120,001 points. This series of points was subjected to 800-point moving-window recurrence analyses using the traditional recurrence and boxed recurrence programs (program3 RQE and RQEB). The results for all eight recurrence variables are superimposed in Figure 6 with blue representing the traditional values and red indicating the boxed values. As shown by the tilted blue arrow in Figure 4, the mean and standard deviation values were computed from the middle third of time series points. Additionally, as shown by the horizontal red arrow in Figure 4, the boxed recurrence variables were shifted left by one-third of the window size (267 points = 800 points/3) to better align with the triangular recurrence variables.

$$\mathbf{X\_{n+1}} = \mathbf{a} \cdot \mathbf{X\_n} \cdot (1 - \mathbf{X\_n}) \tag{7}$$

**Figure 6.** Recurrence variables of the logistic map computed by traditional RQA (blue lines) and boxed RQA (red lines). The logistic time series consisted of 120,001 points for a values ranging from 2.8 to 4.0. Each sliding window epoch consisted of 800 points which were shifted by 10 points (98.75% overlap) to yield 11,920 points per variable. Superimposed values are shown in red; divergent values are shown wherever blue appears. Parameter settings: DELAY = 1; EMBED = 3; NORM = Euclid; WINDOW = 800; SHIFT = 10; RESCALE = absolute; RADIUS = 1; LINE = 2.

subjected to 800-point moving-window recurrence analyses using the traditional recurrence and boxed recurrence programs (program3 RQE and RQEB). The results for all eight recurrence variables are superimposed in Figure 6 with blue representing the traditional values and red indicating the boxed values. As shown by the tilted blue arrow in Figure 4, the mean and standard deviation values were computed from the middle third of time series points. Additionally, as shown by the horizontal red arrow in Figure 4, the boxed recurrence variables were shifted left by one-third of the window size (267 points = 800

Xn+1 = a · Xn · (1 − Xn) (7)

points/3) to better align with the triangular recurrence variables.

In most cases, the red traces overwrite the blue traces, proving excellent agreement with both approaches. However, there are three notable exceptions. First, the ENT values are widely divergent when the logistic system was in its fully periodic modes (period 1,

period 2, period 4, period 8, etc.) over the approximate range of a = 2.8 to 3.6 with spikes just over a = 3.2. Additionally, during the period 3 window (a = 3.83) the 2 entropy values diverge. diverge. Second, the DMAX values superimposed nicely, save during the sliding periodic windows. Because the window size was selected as 800 points, the traditional DMAX

**Figure 6.** Recurrence variables of the logistic map computed by traditional RQA (blue lines) and boxed RQA (red lines). The logistic time series consisted of 120,001 points for a values ranging from 2.8 to 4.0. Each sliding window epoch consisted of 800 points which were shifted by 10 points (98.75% overlap) to yield 11,920 points per variable. Superimposed values are shown in red; divergent values are shown wherever blue appears. Parameter settings: DELAY = 1; EMBED = 3; NORM

In most cases, the red traces overwrite the blue traces, proving excellent agreement with both approaches. However, there are three notable exceptions. First, the ENT values are widely divergent when the logistic system was in its fully periodic modes (period 1, period 2, period 4, period 8, etc.) over the approximate range of a = 2.8 to 3.6 with spikes just over a = 3.2. Additionally, during the period 3 window (a = 3.83) the 2 entropy values

= Euclid; WINDOW = 800; SHIFT = 10; RESCALE = absolute; RADIUS = 1; LINE = 2.

*Entropy* **2022**, *24*, x FOR PEER REVIEW 10 of 18

Second, the DMAX values superimposed nicely, save during the sliding periodic windows. Because the window size was selected as 800 points, the traditional DMAX peaks at 799 diagonal points (just next to the LOI). However, the boxed window was only one-third the size of 800 which is why the boxed DMAX is clipped at 267 diagonal points. peaks at 799 diagonal points (just next to the LOI). However, the boxed window was only one-third the size of 800 which is why the boxed DMAX is clipped at 267 diagonal points. **7. Radius versus Embedding**

#### **7. Radius versus Embedding** In RQA, one of the most difficult parameters to select is the radius. If the radius is too

In RQA, one of the most difficult parameters to select is the radius. If the radius is too low, the number of recurrent points will be too sparse; if radius is too high, the number of recurrence points will be too dense. Indeed, when the radius equals or exceeds the maximum distance in the distance matrix, the recurrence plot will completely saturate (RECUR = near 100% and DETER = near 100% for LINE = 2). Taking advantage of this principle, it is possible to increment the radius and see at which point the entropy values for the triangular RQA and boxed RQA diverge. This was accomplished by using the new boxed RQA program RQSB. Secondly, these divergent points must be a function of the embedding dimension. To quantify these ideas, a single 800-point window was selected from the logistic map (a = 3.9074 to 3.9873) as shown in Figure 7. This window represents a purely chaotic window with no regular periodicities present. The signal looks stochastic, but because it is derived by iterating a nonlinear equation it is actually fully deterministic. This does not mean that the DETERM is 100%, however, because the various orbits are not equally close to one another. Increasing the radius, includes more and more orbits within the recurrence zone. low, the number of recurrent points will be too sparse; if radius is too high, the number of recurrence points will be too dense. Indeed, when the radius equals or exceeds the maximum distance in the distance matrix, the recurrence plot will completely saturate (RECUR = near 100% and DETER = near 100% for LINE = 2). Taking advantage of this principle, it is possible to increment the radius and see at which point the entropy values for the triangular RQA and boxed RQA diverge. This was accomplished by using the new boxed RQA program RQSB. Secondly, these divergent points must be a function of the embedding dimension. To quantify these ideas, a single 800-point window was selected from the logistic map (a = 3.9074 to 3.9873) as shown in Figure 7. This window represents a purely chaotic window with no regular periodicities present. The signal looks stochastic, but because it is derived by iterating a nonlinear equation it is actually fully deterministic. This does not mean that the DETERM is 100%, however, because the various orbits are not equally close to one another. Increasing the radius, includes more and more orbits within the recurrence zone.

**Figure 7.** Fully chaotic mode of the logistic equation within an 800-point window (a = 3.9074 to 3.9873). The logistic x variable fluctuates wildly within the bounds of 0.0 (minimum) to 1.0 (maxi-**Figure 7.** Fully chaotic mode of the logistic equation within an 800-point window (a = 3.9074 to 3.9873). The logistic x variable fluctuates wildly within the bounds of 0.0 (minimum) to 1.0 (maximum).

mum). The computational results are shown in Figure 8, where traditional entropy calculations are displayed by the blue lines and the boxed entropy calculations are displayed by the red lines (programs RQS and RQSB). As can be seen, both entropy values track very closely until the density of the recurrence plots becomes too great at very high radius values. At these turning points the two entropy values begin to diverge as identified by the black dots in each trace. As the embedding dimension increases the turning point shifts to the left. For the five increasing embedding dimensions, the five decreasing radius values at the turning points are, respectively: RADIUS = 94.94%, 85.44%, 80.95%, 79.04% and 77.86%. At the various turning points, all entropy values, triangular and boxed, are bunched within a rather tight or narrow range (ENT = 7.438 to 7.913 bits/bin) albeit at high

values. Additionally, as expected, when the radius value reaches 100% of the maximum distance, the two entropy values depart to their separate corners: regular entropies to their maxima of 9.64 bits/bin; boxed entropies to their minima of 0.00 bits/bin. What can be concluded from these results is that the triangular and boxed entropy values are so similar that the triangular RQA computations for entropy remain valid just as long as the systems under study are not fully periodic in the absence of noise. values. Additionally, as expected, when the radius value reaches 100% of the maximum distance, the two entropy values depart to their separate corners: regular entropies to their maxima of 9.64 bits/bin; boxed entropies to their minima of 0.00 bits/bin. What can be concluded from these results is that the triangular and boxed entropy values are so similar that the triangular RQA computations for entropy remain valid just as long as the systems under study are not fully periodic in the absence of noise.

The computational results are shown in Figure 8, where traditional entropy calculations are displayed by the blue lines and the boxed entropy calculations are displayed by the red lines (programs RQS and RQSB). As can be seen, both entropy values track very closely until the density of the recurrence plots becomes too great at very high radius values. At these turning points the two entropy values begin to diverge as identified by the black dots in each trace. As the embedding dimension increases the turning point shifts to the left. For the five increasing embedding dimensions, the five decreasing radius values at the turning points are, respectively: RADIUS = 94.94%, 85.44%, 80.95%, 79.04% and 77.86%. At the various turning points, all entropy values, triangular and boxed, are bunched within a rather tight or narrow range (ENT = 7.438 to 7.913 bits/bin) albeit at high

*Entropy* **2022**, *24*, x FOR PEER REVIEW 11 of 18

**Figure 8.** Line entropy values as a function of increasing radius values for five different embedding dimensions (E1, E2, E3, E4 and E5) for regular RQA (blue lines) and boxed RQA (red lines) computations. The black points indicate the radius value when the two entropy values diverge (turning points). The full range of radius values (top traces) are magnified for the upper quartile of radius values (lower traces) to show better discrimination of the parameter settings. Parameter settings: DELAY = 1; EMBED = 5; NORM = Euclid; WINDOW = 500; RESCALE = maximum distance; RADIUS = 10%; LINE = 2.

#### **8. Dow Jones Industrial Averages**

As shown above, the noise-free sinewave signal presented above represents an extreme example of divergent entropy computations for verses triangular recurrences verses boxed recurrences (ENT = 3.907 vs. 0.000 bits/bin). Adding noise to the sinewave converges the two entropy values (ENT = 3.000 vs. 2.807 bits/bin). However, what happens to all eight recurrence variables when the recurrence area is constrained (masked) to the tilted box as compared against the standard recurrence triangular recurrence area?

To answer this question, real-world data (noisy) were studied by both recurrence methods and the results were compared. The input time series consisted of downloaded scores of the Dow Jones Industrial Average for over a decade [17]. Recurrence quantifications were computed within a sliding window (programs RQE and RQEB) and the results are shown in Figure 9. Since the window size was 500 points, the length of the box side was 167 points (500/3) and %Box was 44.58% (see Equations (4)–(6)). Again the mean and standard deviation were computed for the middle third of the time series within each window, and the recurrence variables were left time shifted as explained in Figure 4.

**Figure 9.** The mean values (left) and standard deviation values (right) are plotted for both the triangular area (blue) and tilted-box area (red). The means line up nicely, but there is less agreement for the standard deviation values. The raw Dow Jones scores are superimposed in the graph of means (green). Second, nonlinear descriptors of the Dow Jones scores are shown in the remaining eight panels (Figure 9). Each recurrence variable will be reviewed one by one, remembering that blue curves refer to triangular recurrences while red curves refer to tilted-box recurrences. **Figure 9.** The mean values (left) and standard deviation values (right) are plotted for both the triangular area (blue) and tilted-box area (red). The means line up nicely, but there is less agreement for the standard deviation values. The raw Dow Jones scores are superimposed in the graph of means (green). Second, nonlinear descriptors of the Dow Jones scores are shown in the remaining eight panels (Figure 9). Each recurrence variable will be reviewed one by one, remembering that blue curves refer to triangular recurrences while red curves refer to tilted-box recurrences.

#### **9. Square Wave Pulse**

To better understand timing and amplitude relationships among all eight recurrence variables computed by the standard RQA versus the boxed RQA, the two methods were compared using a simple square wave input signal. This signal consisted of 7500 points in which the first 1500 points and last 4000 points were set to 0, but the middle 2000 points were set to 10. The first and second transitions were abrupt (pulse up: 0 to 10; pulse down: 10 to 0, respectively). Figure 10 presents the results of this square wave as processed by both the traditional RAQ and boxed RQA computations for all eight recurrence variables (programs RQE and RQEB). The sliding window approach was taken at the highest sensitivity (only 1 point between epochs) to allow for clear timing delineations. Again It should be noted that the boxed recurrences were aligned with the regular recurrences, plotting both time series from point 1. That is the boxed recurrences were plotted from point 1, not point 167 (or 500/3). Additionally, the double vertical green lines in each panel show the exact onset and offset of the square wave pulse.

The first observation is that width of the means and standard deviations are broader for the standard RQA than for the boxed RQA. This is because the tilted box spans fewer pulse points than the triangular space. Both means start from the same base (0.0) and rise to the same height (10.0), but the rise in the boxed mean begins later and has a greater slope. It can also be noted that the two standard deviations both reach identical maxima.

Second, prior to the onset and offset of the square pulse, the paired recurrence variables have peaks or nadirs that line up in time. Additionally, as is necessary, the interval between peaks or nadirs is exactly 2000 points which is the width of the square wave pulse.

Third, the contours and amplitudes of the paired recurrence variables are not the same. This is explained by the fact that the boxed recurrence area consists of only 55,611 pixels which is but 44.58% of the triangular recurrence area (124,750 pixels) (see Equations (4)–(6)).

Fourth, during the steady-state phases of the input signal (string of 0 s pre- and postpulse; string of 10 s during the pulse), the regular entropy is pegged at 8.963 bits/bin (incorrect) whereas the boxed entropy is minimized to 0.000 bits/bin (correct). Only during the two transition phases do the two entropy values approach each other, but never meet. With the exception of these two entropy values and possibly the LAMINAR and TRAPTIME variables, the other five RQA variables exhibit similar directional changes during the two transition phases.

**Figure 10.** Comparison of traditional and boxed recurrences on a single square wave of 7500 points for 2 linear and 8 nonlinear recurrence variables. Triangular RQA computations (blue traces) and boxed RQA computations (red traces) are performed within a sliding window of 500 points offset by a single point between epochs for a total of 6977 epochs. However, only the first 5,000 epochs are plotted. The onset of the rising square wave is indicated by the first vertical green line (pulse on); the end of the falling square wave is indicated by the second vertical green line (pulse off). Both time series were aligned to point 1 (see text). Parameter settings: DELAY = 1; EMBED = 5; NORM = Euclid; WINDOW = 500; SHIFT = 1; RESCALE = maximum distance; RADIUS = 10%; LINE = 2. **Figure 10.** Comparison of traditional and boxed recurrences on a single square wave of 7500 points for 2 linear and 8 nonlinear recurrence variables. Triangular RQA computations (blue traces) and boxed RQA computations (red traces) are performed within a sliding window of 500 points offset by a single point between epochs for a total of 6977 epochs. However, only the first 5000 epochs are plotted. The onset of the rising square wave is indicated by the first vertical green line (pulse on); the end of the falling square wave is indicated by the second vertical green line (pulse off). Both time series were aligned to point 1 (see text). Parameter settings: DELAY = 1; EMBED = 5; NORM = Euclid; WINDOW = 500; SHIFT = 1; RESCALE = maximum distance; RADIUS = 10%; LINE = 2.

#### **10. Discussion and Conclusions**

Back in 2015, Marwan and Webber [18] discussed in detail the mathematical and computational foundations of recurrence quantifications. However, no mention was made about masking the recurrence matrix to modify the triangular border truncations of diagonal lines. Later, both of these authors independently posited a tilted rectangle [9] or tilted box (present paper) masking of the recurrence matrix such that long diagonal lines are clipped to identical lengths thereby minimizing the line entropy computations. It should be emphasized that both entropy calculations, triangular entropies and boxed entropies, are accurate, but have different meanings. In the traditional case, line-entropy values really describe the actual probability distributions of truncated lines. However, in the tilted box case, line-entropy values better correspond to the underlying dynamic itself. That is, low entropy values are related to simple (or random) systems, whereas high entropy values are usually correlated with more complex systems, but not necessarily always. As mentioned previously, Censi et al. [11], Eroglu et al. [12], Corso et al. [13] and Leonardi [14] all have thoughts about revised entropy computations. The question is not which method is the best, but rather, based upon the assumptions of each approach, how do the various entropy computations perform on systems of all sorts. This becomes an invitation for performing comparative studies applying these differing methodologies to representative systems of interest (beyond the scope of this paper).

For example, Leonardi [14] has a very nice summary description of the meaning of entropy, especially the information entropy of Shannon [8]. Information entropy values approaching maximum entropy (such as for random systems) contain high levels of information, are complex in nature, and unpredictable in principle. However, information entropy values approaching minimum entropy (such as for highly periodic systems) contain low levels of information, are simple in construct, and can be very predictable. This last statement is qualified for white-noise processes which present with very sort line structures which translate to very low entropy values as well. In short, the higher the entropy value, the more interesting is the system. Conversely, the lower the entropy value, the less interesting, even boring, is the system (excepting white noise). So for the sinewave, which is definitely periodic, very linear, and very repetitive (if you have seen one, you have seen them all), the entropy minimum of 0.000 bits/bin best describes the diagonal-line entropy. The tentative conclusion is that boxed recurrence entropies may have an advantage over triangular recurrence entropies. However, confirmation of such a claim must be verified by other correction schemes [9].

Many observations in this paper suggest that entropy values, if not most other recurrence variables, computed from triangular and boxed recurrences may not necessarily be that different after all. First, the addition of subtle noise to a pure sine wave collapses both entropies to nearly the same values insofar that long diagonal lines are disallowed (Figures 1–3). Second, the logistic equation in its non-periodic and chaotic modes assigns very similar quantitative (overlapping) values to all recurrence variables (Figure 6). Third, both entropies track together for the chaotic state of the Logistic system until breaking points of separation occur only at very high radius values (Figure 8). Additionally, fourth, moving-window recurrence analyses of the Dow Jones Industrial Average data reveal correlations between triangular and boxed recurrence variables that have similar qualitative (REC, DET and LAM) or quantitative (ENT, TREND and TRAPTIME) profiles. With the exception of two poorly correlated variables (DMAX and VMAX), such patterns would lead to singular, not diverse, interpretations of the underlying dynamic (Figure 9).

As an aside, it will be appreciated that the two tilted maskings of the recurrence matrix are qualitatively similar for both Kramer and Marwan [9] and the present paper. What distinguishes the two approaches quantitatively, however, is that the area of the masked rectangular recurrence plot is a full 50% in the limit of the triangular matrix whereas the masking of the tilted box is only 44%. This 8% difference may constitute an advantage for the former insofar that it captures a greater proportion of recurrent points. This slight edge may be a fine tune difference between masked versus unmasked recurrence plots and quantifications. This hypothesis could be studied using various signals.

Finally, the moving-window recurrence analyses of the square wave pulse are provided basically for heuristic purposes alone (Figure 10). The input signal is very clean (noise-free) and two instantaneous jumps occur in the time series (pulse on and pulse off). One can visualize the unmasked triangular window and the masked boxed window both moving through the vertical lines. The responses in all eight recurrence variables reflect the two transitional state changes encountered. In addition, the effect of the left shift of the boxed recurrence variables can also be studied clearly. These plots can assist in learning how to read recurrence plot variables when using masked versus unmasked methodologies.

Perhaps the most surprising conclusion of this paper is that box masking of the recurrence matrix is not necessarily called for in most practical situations. Nevertheless, this study reveals that border truncation effects are real when dealing with long (even infinite) time series. If this observation can be verified for different systems it would mean that all previous studies relying on the traditional entropy calculations for systems that are even near periodic, the entropy values and other recurrence variables computed and reported are valid. A good example is one recent master's thesis on which the conclusions relied heavily on traditional triangular entropy values to come to consistent and meaningful conclusions [19].

**Funding:** This research received no external funding.

**Data Availability Statement:** The data reported in this paper can be requested directly from the author.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **Appendix A**

The new software programs used in this paper, all written by the author, are bundled within the suite of RQA 2021 programs. These new programs include RQDB (recurrence display) RQEB (recurrence epoch), and RQES (recurrence scale) which were run to compute the boxed recurrences in this paper. They correspond to traditional triangular recurrence programs RQD, RQE and RQS. All RQA programs are embedded within a single ZIP file [20] which is easily downloaded at no cost. C-source codes will be made available to anyone for the asking.

#### **References**


Notice that it is mandatory to have a method to compose power expansions of the manifold with the map itself. Sometimes, when, the map is explicit, one is able to find a suitable recurrence expression to compute the terms of this composition. If a recurrence is not available one can compute higher order terms by automatic differentiation1 . In the case we are interested, the map is not given explicitly but comes from a numerical integration of a differential equation. Here, the only reasonable strategy seems to use Jet Transport. This technique is based in the idea of, instead of integrating a single point, integrate a function given by its expansion in Taylor series. That is, one transports the jet, the set of derivatives of the function at a given point, up to a given order. It is straightforward to construct an integrator of jets from an integrator of numbers. It is only a matter of replacing the standard floating arithmetic by a polynomial arithmetic [Farrés et al., to appear].

There is another obstacle that can appear when dealing with Poincaré maps: if the orbit is very unstable, the hyperbolic direction may lead to a huge error propagation. This problem can be avoided by using parallel shooting. The idea behind parallel shooting is to enlarge the dimension of the system in order to decrease the time of integration. Let us denote by ϕ tf t0 (x) the solution of the differential equation with initial condition (t0, x) evaluated at time tf . Fix k ∈ N, the number of sections, and set h = T/k. For i = 0, . . . k, we define τi = ih. If m = nk, we define the function F : V 7→ R m

$$F: \begin{pmatrix} \varkappa\_1 \\ \varkappa\_2 \\ \vdots \\ \varkappa\_k \end{pmatrix} \mapsto \begin{pmatrix} f\_k(\varkappa\_k) \\ f\_1(\varkappa\_1) \\ \vdots \\ \vdots \\ f\_{k-1}(\varkappa\_{k-1}) \end{pmatrix}.$$

Here V = U k is an open set of R m, and, for x ∈ U, fj(x) = ϕ τj τj−1 (x). The differential map DF is given by

$$DF = \begin{pmatrix} 0 & \dots & \dots & Df\_k \\ Df\_1 & \dots & \dots & \vdots \\ \vdots & \dots & \dots & \vdots \\ 0 & \dots & Df\_{k-1} & 0 \end{pmatrix}.$$

For x¯ ∈ U, name y = (x1, . . . , xk ) where x1 = ¯x and xj = fj(xj−1) if 1 < j ≤ k. Then:

### 4. THE BICIRCULAR PROBLEM AND THE TRIANGULAR POINTS

The BCP is a perturbation of the RTBP. It is usual to take the units and the synodic coordinates of the Earth-Moon RTBP (see **Figure 2**). The BCP is not coherent, that is, the trajectories followed by the primaries do not obey Newton's laws. This is not an inconvenient since the model has been shown to be useful to describe the dynamics near the triangular points [51]. As a dynamical system, the BCP is a Hamiltonian system with three and a half degrees of freedom, i.e., a non-autonomous periodically time dependent with three degrees of freedom. The Hamiltonian function, written in the RTBP coordinates and units, is given by

$$H = \frac{1}{2} (p\_x^2 + p\_y^2 + p\_z^2) - \chi p\_\mathcal{Y} + \chi p\_\mathcal{x} - \frac{1-\mu}{r\_{\rm PE}} - \frac{\mu}{r\_{\rm PM}}$$

$$-\frac{m\_\mathcal{S}}{a\_\mathcal{S}^2} (\mathbf{y}\sin\theta - \mathbf{x}\cos\theta) - \frac{m\_\mathcal{S}}{r\_{\rm PS}}.\tag{3}$$

Here µ, rPE and rPM denote the same quantities as in (1). Moreover, mS denotes the mass of Sun, aS the averaged semimajor axis of Sun, θ = ωSt, ωS is the frequency of Sun in this system of reference, TS = 2π ωS is its period and finally, r 2 PS = (x−aS cos θ) 2 +(y−aS sin θ) 2 +z 2 . Notice that this Hamiltonian can be splitted in two parts:

$$H\_{BCP} = H\_{RTBP}(X) + H\_S(\theta, X).$$

1 In principle, one could also use numerical differentiation but it is a bad approach in terms of efficiency and accuracy.

The time dependent part contains two terms, the Coriolis effect due to the rotating frame of coordinates and Sun's gravitational potential. The Taylor expansion of the potential starts with

$$\frac{1}{a\_{\mathcal{S}}} \left( 1 + \frac{\varkappa \cos \theta - \wp \sin \theta}{a\_{\mathcal{S}}} \right)^2$$

.

Therefore, the Hamiltonian, if we truncate the Sun's potential at linear order is written as

$$H\_{\rm BCP}^{<2} = H\_{\rm RTBP} - \frac{m\_{\rm S}}{as}.$$

So, the Coriolis term and the truncated Sun's potential cancel out and the dynamics is the one of the RTBP. This is to say that the contribution due to Sun's potential starts at order two, that is, the BCP is a periodic time dependent perturbation with size O( mS a 3 S ) ≈ 0.0056. Anyhow, it is large enough to produce remarkable changes on the dynamics, especially near the triangular points. In **Figures 3** (left), **6** (left) we display continuations from the RTBP to the BCP. The vertical axis in these plots represent an artificial parameter ε which multiplies the mass of Sun. Therefore, when ε = 0 the model corresponds to the RTBP and, when ε = 1 the model corresponds to the BCP. We shall comment these figures in detail in the next sections.

### 4.1. Dynamical Equivalents of the Triangular Points

First of all let us mention that, due to a symmetry, the dynamics near L4 is the same as the dynamics near L5 (in fact, this symmetry maps orbits in the region y > 0 to orbits in the region y < 0). A feature of the BCP to be stressed is that the region around the geometrically defined triangular points is unstable. The influence of Sun's potential is enough to produce a bifurcation in the periodic orbit that replaces L4 (i.e., L5). It is well known [51] that each triangular point is replaced by three periodic orbits with the same period as Sun. One small and unstable (the actual replacement of L4) and two which are stable. We have named these orbits PO1, PO2, and PO3. See **Figure 3** (left) a continuation diagram from the RTBP to the BCP. The two additional periodic orbits are produced by an imperfect pitchfork bifurcation (i.e., a pitchfork bifurcation broken due to a loss of symmetry).

One may ask which is the model that displays the perfect bifurcation and which is the broken symmetry. To address this question we take a look at the order two of the Taylor expansion of Sun's gravitational potential. We have

$$H\_{\mathcal{S}}^2(\theta, \chi, \boldsymbol{\chi}, \boldsymbol{z}) = \frac{1}{a\_{\mathcal{S}}^3} (\frac{3}{2}T(\boldsymbol{\chi}, \boldsymbol{\chi}, \theta)^2 - \frac{1}{2}(\boldsymbol{\chi}^2 + \boldsymbol{\chi}^2)).$$

We have named T(x, y, θ) = −x cos θ + y sin θ. We would like to stress again that H2 S is the first contributing non-autonomous term in the model due to the cancellation produced by the Coriolis acceleration. This term is invariant under the symmetry (x, y, x, θ) 7→ (x, −y, x, −θ). The order three of the expansion is given by:

$$H\_{\mathbb{S}}^3 = \frac{1}{a\_{\mathbb{S}}} \left( \left( \frac{\rho}{a\_{\mathbb{S}}} \right)^3 \frac{\frac{5}{2}T^3 - \frac{3}{2}T}{a\_{\mathbb{S}}\rho} \right).$$

Here ρ 2 = x 2 + y 2 . The polynomial in T is no longer even. This breaks the symmetry and, hence, the pitchfork bifurcation. The non-autonomous model that displays the perfect bifurcation is:

$$H\_{\rm BCP}^{\leq 2} = H\_{\rm RTBP}(X) + H\_{\rm S}^2(X, \theta).$$

The perfect (non broken) pitchfork bifurcation in **Figure 3** (left, curve in red) shows the continuation diagram from the RTBP to this simplified version of the BCP. Due to the symmetry, periodic orbits PO2 and PO3 only differ on the phase on the orbit.

### 4.2. Phase Space of the Stroboscopic Map Near the Triangular Points

The three periodic orbits appear as fixed points of the stroboscopic map. We recall that their linear normal behavior is of type saddle×center×center for PO1 and totally elliptic for PO2 and PO3. There are several ways to justify that, from the elliptic directions of each fixed point, there is a family of invariant curves whose frequency tends to the normal modes of the fixed points [44, 52].

Therefore, we have a family of invariant curves for each elliptic direction, that is, two for PO1 (HF1 in the horizonal plane and VF1 in the vertical direction), three for PO2 (HF2F1 and HF2F2 are horizontal, and VF2 is vertical) and three for PO3 (HF3F1 and HF3F2 horizontal, and VF3 vertical). The remaining eigendirection of PO1 is hyperbolic. There exist stable and unstable one-dimensional invariant manifolds associated to these hyperbolic directions. Initial conditions near the triangular points shadow the unstable manifold which wonder some time around the periodic orbits PO2 and PO3 and finally abandon the vicinity of the triangular points. These manifolds are of special interest if one plans to put or take out objects near L4 and L5. The stable and unstable manifolds related to PO1 can be computed up to high order directly on the stroboscopic map (section 3.2 and [Farrés et al., to appear]). In **Figure 4** we observe a projection of the phase portrait of the map. The three points displayed with crosses correspond to PO1 (in the middle), PO2 and PO3. It is displayed as well, semianalytical approximations of the stable and unstable invariant manifolds. We have used an approximation of order 64. The width curve are the pieces given by the parameterization. The thin curve correspond to some iterations of these pieces. It can be observed, also, some invariant curves growing from PO2 and PO3. These invariant curves are totally elliptic near the fixed points.

### 4.3. Regions of Effective Stability Near the Triangular Points

As we have observed, the triangular points are replaced by three periodic orbits, one of them unstable. This is the reason why the BCP is an interesting model [12]. Indeed, the unstable manifold of the triangular periodic orbit takes initial conditions away from the vicinity of the triangular points. However, we can pursue on the seek for regions of (effective) stability out of the plane of motion of the primaries. The mechanism that suggest the existence of regions of effective stability is the stickiness of normally elliptic low dimensional invariant tori. See [41, 52]. As we discussed before, there are families of invariant tori emanating from the periodic orbits PO2 and PO3. These families are elliptic close enough to the fixed points. This results on two small regions of effective stability in the plane of motion of the primaries related to the totally elliptic orbits.

We put our attention on the vertical families (one for each orbit) of invariant tori. In **Figure 3**, Right, we display how these families vary when they grow out of the plane of motion of the primaries. We observe that the three families display a broken pitchfork symmetry, analogous to the one of the periodic orbits. The linear normal behavior of the tori is the same as the periodic orbit near the plane. As a consequence of the pitchfork bifurcation, at some some distance of the plane, the surviving family is totally elliptic. Therefore, the tori are sticky and regions of effective stability are to be expected. We label these families VF1, VF2 and VF3 after the corresponding periodic orbits. The families VF1 and VF2 are connected, as it happens with PO1 and PO2. On the other hand VF3 reaches high amplitudes in the (z, pz) plane. It is known that, skipping resonances, the three families have the same stability as the corresponding periodic orbits. Therefore, the tori of VF1 have hyperbolic directions while the ones of VF2 are normally elliptic (except for small intervals of instability produced by resonances involving internal and normal frequencies). Recall that both families are connected and the change of stability takes place at a turning point. The tori of VF3 are normally elliptic up to very high values (again, except for resonances).

Normally elliptic lower dimensional tori induce regions of effective stability. Numerical estimations of the shape and the size of these regions show that, in the case of VF2, the regions are small and narrow while in the case of VF3 large regions exist for sufficiently high values of the vertical amplitude. In **Figure 5** we show two stability regions out of the horizontal plane. These regions seem to persist in the real model for time spans of 1,000 years [14, 53]. The effect of Solar Radiation Pressure on the effective stability regions of the BCP is discussed in Jorba-Cuscó et al.[54].

Let us explain how **Figure 5** is obtained. Each of the vertical family of invariant tori can be identified by the value of the coordinate pz when z = 0 and pz > 0. Denote by a(pz) ∈ R 6 the coordinates of the point that identifies a torus. We have to select a set of initial conditions near a(pz) and integrate them for a long time span. Let us be more precise on how to select the initial conditions. We use a two dimensional grid, the coordinates z, px, py and pz will be fixed by the corresponding values of a(pz). To adapt to the shape of the regions, we use a polar-like grid, centered at Earth:

$$\begin{cases} \varkappa\_{i\bar{j}} = r\_i \cos \alpha\_{\bar{j}} + \mu, & r\_{\bar{i}} = 1 + ih\_r, \\ \nu\_{i\bar{j}} = r\_i \sin \alpha\_{\bar{j}}, & \alpha\_{\bar{j}} = 2\pi j h\_\alpha, \end{cases}$$

where hr and hα are used to control the density of the grid. The computation goes as follows. Take a point of the grid and integrate the vector field 15000 Moon revolutions. At each integration step, we test if there is a collision with Earth or Moon. If there is a collision, or the coordinate y becomes negative, we stop the integration (Recall that we are interested in the points that remain close to L4). We have used hr = 0.001 and hα = 0.0002. The difference on the sizes of these small quantities is aimed to produce a nearly squared grid.

### 4.4. A Weakness of the BCP

The translunar point is one of the most interesting locations at the Earth-Moon system. The reason is that L2 seems suitable to observe the far side at Moon. Taking into account that, a

natural criticism to the BCP is that it does not have a dynamical replacement of L2. In **Figure 6** (left) it is displayed a continuation of L2 as a periodic orbit from the RTBP to the BCP. The periodic orbits are identified by their coordinates as fixed points of the Stroboscopic map. These orbits have been computed by means of the parallel shooting method. Again, the vertical axis is an additional parameter, ε, multiplying the mass of Sun. The point L2 is the middle crossing of the characteristic curve with the homotopy level corresponding to the RTBP, at the bottom. The other two points of the RTBP correspond to the same 1:2 resonant planar Lyapunov orbit but with different initial times. We observe that the continuation of L2 reaches a turning point and it never reaches the homotopy level of the BCP. The result is that the translunar dynamical structure is lost in the BCP. This suggest that a more complex model needs to be used to analyze the natural behavior near the translunar point. The resonant Lyapunov orbit can be continued to the BCP. The result is a large orbit that remains away from the translunar point, see **Figure 6** (right).

### 5. THE QUASI-BICIRCULAR PROBLEM AND THE COLLINEAR POINTS

The quasi-bicircular solution of the Earth-Moon-Sun system is planar i.e., the three bodies move in the same plane. After the quasi-bicircular solution is computed one can write the equations of motion of the test particle, prescribing the quasi-bicircular solution as motion for the primaries. It is usual to compute the quasi-bicircular solution in the Jacobi frame, however, if one has the purpose of describing the dynamics in the Earth-Moon vicinity, it is suitable to use the frame of coordinates corresponding to the Earth-Moon RTBP. To do so, one has to perform three different transformations. First, one has to use a translation to move the origin from the global barycenter to Earth's and Moon's center of masses. Second, one has to use a rotating (synodic) frame to keep Earth and Moon fixed on the horizontal axis. Third, the unit of length is scaled so the distance between Earth and Moon is equal to one. The units of mass and time which are usually selected in the Earth-Moon RTBP can be imposed already in the Jacobi formulation of the Three Body Problem.

The resulting model is a Hamiltonian system with three and a half degrees of freedom. The Hamiltonian function can be written as

$$H = \frac{1}{2}\alpha\_1(p\_x^2 + p\_y^2 + p\_z^2) + \alpha\_2(p\_x x + p\_y y + p\_z z) + \alpha\_3(p\_x y - p\_y x)$$

$$+ \alpha\_4 x + \alpha\_5 y - \alpha\_6 \left(\frac{1 - \mu}{r\_{pe}} + \frac{\mu}{r\_{pm}} + \frac{m\_S}{r\_{ps}}\right),\tag{4}$$

TABLE 2 | Continuation of the low order resonant orbits from the RTBP to the QBCP.

The first column contains the label of the orbits corresponding to the RTBP. The second column contains the order of the resonance. The third columns contains number of bifurcating orbits. The fourth column contains the label of the orbits corresponding to the QBCP. See [15] for more details. See text for the color code.

We only put three for each orbit. The rest are given by their inverses due to the symplectic character of the stroboscopic map.

where, r 2 pe = (x − µ) 2 + y 2 + z 2 , r 2 pm = (x − µ + 1)2 + y 2 + z 2 , r 2 ps = (x − α7) 2 + (y − α8) 2 + z 2 , and for i = 1, . . . , 8 αi : T 7→ R are periodic functions. That is,

$$\alpha\_i(\theta) = a\_0^i + \sum\_{k \ge 0} a\_k^i \cos k\theta + \sum\_{k \ge 0} b\_k^i \sin k\theta. \tag{5}$$

Here, θ = ωSt and ωS is the frequency of Sun. Moreover, αi is odd for i = 1, 3, 4, 6, 7 and even for i = 2, 5, 8. Obviously one can only have a numerical approximation of these functions. In this case, we take advantage on the computations done in [15] and take the same values for the Fourier coefficients of the periodic functions αi 's. To end, and taking into account the properties of the functions αi 's, it is easy to see that the Hamiltonian function (4) has the symmetry (θ, x, y, z, x˙, y˙, z˙) 7→ (−θ, x, −y, z, −˙x, y˙, −˙z), x˙ = px + y, y˙ = py − x, z˙ = pz.

The meaning of these periodic functions is the following:

1. (α7, α8, 0) is the position of Sun in the plane of motion of the primaries.

### 5.1. Dynamical Equivalents of the Collinear Points

In this section we give some words about the minimal periodic orbits that replace the collinear points in the QBCP. In **Figure 7** we display the dynamical equivalents, from left to right, of L1, L2 and L3. We observe that the orbits replacing L1 and L2 are small, their maximal distance to the corresponding equilibrium point is of order O(10−6 ). As the original equilibrium points, the linear normal behavior of these orbits is of type saddle×center×center. In **Table 1** we display the eigenvalues of each orbit. We notice that the unstable direction of L1 (of order 108 ) and the unstable direction of L2 (of order 106 ) are large and this implies huge propagation of error near these orbits. On the other hand, the dynamical equivalent of

L3 has a very weak unstable direction, at least compared to the other two.

### 5.2. Resonant Orbits of Low Order

As the QBCP is a TS-periodic system, the simplest invariant objects are TS-periodic orbits. We already have mentioned that the equilibrium points are replaced by these periodic orbits of minimal period. Periodic orbits of the RTBP whose frequency is resonant with the one of Sun also persist as TS-periodic orbits in the QBCP. The Lyapunov and Halo families of periodic orbits related to the equilibrium points L1 and L2, are a source for these kind of resonant orbits. In contrast with the families related with L3, the families of the two first libration points are nourished with low order resonant orbits. A relation of low order resonant periodic orbits of the RTBP can be found in [15]. In [29] the authors show the ranges for the admissible periods for each family. The families related to L3 are of relatively large period and there are not many periodic orbits whose frequency are in low order rational relation with the frequency of Sun. There is, however, a 1:1 resonant periodic orbit near the end of the vertical family. This orbit is enormous in size and cannot be considered in the vicinity of L3. In **Table 2** details of the continuations of low order resonant orbits from the RTBP to the QBCP are given: The first column corresponds to resonant periodic orbits of the RTBP. The label in this first column consist in three numbers that encode each orbit. The first is a zero and indicates that the orbit belongs to the RTBP (this is intended to distinguish them from the orbits

#### TABLE 4 | Coefficients of the functions αj , j = 1, . . . , 8, in (5).

in the last column corresponding to the QBCP). The second number refers to the libration point related to each orbit (all of them belong to Lyapunov and Halo families related to L1 and L2). The third number is just an enumeration. The second column indicates the order of the resonance. We stress that the influence of Sun is relevant enough to produce bifurcating orbits in each of the continuations. The third column shows how many orbits bifurcate from the original ones when they are continued to the QBCP. Finally the last column contains the labels of the resulting orbits in the QBCP. **Table 2** can be found originally in Andreu [15]. We have added the order of the resonance and the color code to indicate the linear normal behavior of each orbit. Labels in blue stand for orbits of type saddle×center×center. Labels in green denote linear character of the kind saddle×saddle×center. Names in cyan denote totally hyperbolic orbits. The color yellow denotes totally elliptic orbits. The continuation for the orbits in red do not reach the homotopy level of the QBCP and, therefore, are not considered.

### 5.3. High Order Approximation of the Unstable Manifolds of the Collinear Periodic Orbits

This section is devoted to the results of implementing the algorithm explained in section 3.2 to the dynamical equivalents of the collinear points. **Figure 8** shows pieces of the stable (dashed) and unstable (solid) manifolds related to the three

Due to the symmetries of the model, each αj only contains either sin or cos terms, so we only list either the ak or bk coefficients.

collinear periodic orbits (the other branches can be obtained by symmetry). From left to right: the one related to L1, the one related to L2 and the one related to L3. We would like to remark that these pieces are obtained directly from the evaluation of the approximation (of order 64) of the manifolds. The error is controlled by checking that the contribution of the last term of the approximating polynomial is small. It is also checked that all the points in the pieces verify the invariance equation with high accuracy. We can observe that these approximations already give large excursions far away from the collinear points. Especially in the case of L3, where the piece of the manifold passes very close to the triangular points. The axes of **Figure 8** show the x and y values. These pieces can be mapped through the stroboscopic map to obtain larger pieces of the manifolds if it is necessary. The point of giving high order approximations of the manifold is that, just a fewer number of iterates are necessary. For the computation of the manifold related to L3, a simple shooting method has been used. Indeed, the instability associated to this libration orbit is very weak. For the computation of the manifold related to L2, multiple shooting is required. We have used two sections. For the computation of the manifold related to L1, the most unstable one, we have used a single shooting strategy but with an extended precision arithmetic of 128 bits. This last approach makes the program far slower but very simple to code. In **Figure 9** (left) we show the resonant orbit 2G of **Table 2**. We display also (right) the stable (dashed) and unstable (manifolds).

### 6. CONCLUSIONS

We have presented two alternatives to the RTBP for the study of the motion of a test particle in the Earth-Moon system. Both models, the BCP and the QBCP, depend periodically on the time. We use the so-called stroboscopic map to study the minimal periodic orbits of the systems and the invariant manifolds related to them.

The BCP is as useful model for the study of the triangular points. The simplicity of the vectorfield is a strong point, especially in problems related to effective stability where massive integrations are mandatory. We have also stressed its weakness: it is not suitable to understand the dynamics around the collinear points. The BCP is useless to describe the vicinity of the translunar point.

### REFERENCES

We have used the parametrization method to obtain high order approximations of the unstable manifolds related to the minimal periodic orbits that replace the collinear points in the QBCP. This is helpful to design long excursions between the two primaries and the collinear points. The main novelty is that we have computed the manifolds directly on the stroboscopic map. The QBCP is a complicated model with a numerically computed vectorfield. This makes it a bad candidate (in front of the BCP) to be the model used to face the problems involving massive simulations related to the triangular points.

We would like to stress that the BCP should be used to face problems related to the triangular points. Especially if this problems involve large time integrations to seek for regions of practical stability. The QBCP should be used when dealing with problems involving the collinear points.

### 7. TECHNICAL DETAILS

All the computations appearing in the Figures of this paper, also the ones which appear in the literature, have been performed by the authors. The integrations for the RTBP, the BCP and the QBCP have used a Taylor method with variable order and stepsize. The demanded accuracy for the standard double precision has been 10−16. The computations in multiple accuracy have been done using the library mpfr. The LAPACK library has also been used for some computations related to linear algebra. The rest of the programs have been written by the authors in C and C++ languages from the scratch. **Table 3** and **Table 4** contains the values of the parametres used for the computations.

### AUTHOR CONTRIBUTIONS

AF, ÀJ, and MJ-C split, more or less, equally the tasks regarding the computation and writing of the paper.

### ACKNOWLEDGMENTS

The authors would like to thank J. Gimeno and N. Miguel for their help. This work has been supported by the Spanish grants MTM2015-67724-P (MINECO/FEDER), MDM-2014-0445 and the Catalan grant 2017 SGR 1374. The project leading to this application has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 734557.

International Sysmposium on Spacecraft Flight Dynamics, ESTEC. Noordwijk: ESA Publications Division (1991). p. 35–41.

With Three or More Degrees of Freedom, NATO Advance Science Series C: Mathematical and Physical Science, Held in S'Agaró, Kluwer Academic Publishers, Dordrecht (1999). p. 430–434.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Jorba-Cuscó, Farrés and Jorba. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.