A Piecewise Linear Approach for Implementing Fractional-Order Multi-Scroll Chaotic Systems on ARMs and FPGAs

: This manuscript introduces a piecewise linear decomposition method devoted to a class of fractional-order dynamical systems composed of piecewise linear (PWL) functions. Inspired by the Adomian decomposition method, the proposed technique computes an approximated solution of fractional-order PWL systems using only linear operators and specific constants vectors for each sub-domain of the PWL functions, with no need for the Adomian polynomials. The proposed decomposition method can be applied to fractional-order PWL systems composed of n th PWL functions, where each PWL function may have any number of affine segments. In particular, we demonstrate various examples of how to solve fractional-order systems with 1D 2-scroll, 4-scroll, and 4 × 4-grid scroll chaotic attractors by applying the proposed approach. From the theoretical and implementation results, we found the proposed approach eliminates the unneeded terms, has a low computational cost, and permits a straightforward physical implementation of multi-scroll chaotic attractors on ARMs and FPGAs digital platforms.


Introduction
Fractional calculus is the branch of mathematics that studies the non-integer order integrals and derivatives, called fractional derivatives and fractional integrals [1].Nowadays, we can found several applications of fractional calculus in diverse fields of science [2][3][4][5] and engineering [6][7][8][9].Among them, fractional-order chaotic systems are of great interest because they represent the intersection of chaos theory with fractional calculus [10,11], forming a dynamical system with non-local characteristics and high complexity [12][13][14][15].
The interest in fractional-order systems is driven by their advantages over traditional integer-order systems.While integer-order systems have been fundamental in understanding the chaotic dynamics, as demonstrated by foundational works such as those by Pecora and Carroll on synchronization of chaotic systems [16] and Ott, Grebogi, and Yorke on controlling chaos [17], they have limitations in describing dynamical systems with memory effects [18].On the other hand, Fractional-order systems, possess memory properties and hereditary characteristics, which are absent in integer-order systems, and hence, they provide a more comprehensive framework for modeling and analyzing the non-local interactions of dynamical systems.These properties also improve the effectiveness for chaos-based applications as follows, • Extra-degree of fredom.The most significant advantage of fractional-order systems is the extra degree of freedom provided by the non-integer order derivatives.This additional parameter allows for greater flexibility in modeling and control.In applications such as chaos synchronization, it is possible to tune and optimize the system behavior by varying the fractional order, thereby achieving precise synchronization schemes and control outcomes [18].• Memory Properties and Hereditary Characteristics.Fractional-order systems possess memory properties and hereditary characteristics, which are absent in integer-order systems.These properties enhance the modeling of systems with memory effects and improve the effectiveness of chaos-based applications.For instance, the memory effects can introduce additional layers of complexity and unpredictability to the chaotic behavior, which can be used to create robust security applications [19,20].

•
Larger Key-space.Fractional-order chaotic systems offer a significantly larger keyspace compared to integer-order systems due to the additional degrees of freedom provided by the fractional derivatives [21].
Building on these advantages, fractional-order multi-scroll chaotic systems have become an exciting research topic because their unique attributes can be used to enhance the performance of practical applications, as demonstrated by recent works in areas such as image and data encryption [22][23][24], audio encryption [25], Internet of Things (IoT) [26,27], random number generators [28], artificial neural networks (ANN) [29], to mention a few.As well known, multi-scroll chaotic attractors present numerous scrolls in phase space, which can be propagated on a 1D line, 2D lattice, and 3D grid [30].This kind of attractor is commonly generated by staircase functions, sawtooth functions, switching functions, saturated series functions, etc., which are generalized as piecewise linear (PWL) functions in many dynamical systems.As a consequence of those features, a multi-scroll chaotic attractor may have much more abundant dynamic behaviors and better unpredictability since its trajectory can randomly jump among multiple scrolls [31][32][33], which has contributed to boosting the security level in data encryption.Additionally, fractional derivatives increase complexity in both time and frequency domains [34][35][36], making them an excellent option for expanding the encryption keys space for secure communications.
However, when fractional-order multi-scroll chaotic systems are used in engineering, we face two challenges.The first one is that we are invariably required to generate chaotic time series.This task is not trivial in real-world scenarios because standard numerical methods can be computationally expensive and require significant storage space [37][38][39].The other issue is related to the fact that the approach to computing a solution for the underlying system should be viable and flexible to be implemented with distinct digital devices [40][41][42].
Therefore, we need to have a method that can efficiently compute an approximate solution of fractional-order PWL systems.The proposed approach should also preserve the computational cost while implementing multi-scroll chaotic attractors on advanced RISC machines (ARM) processors and field-programmable gate arrays (FPGAs).This manuscript goes in the same direction and presents a positive answer to these challenges.
In the literature, we found predominantly three solving methods to implement fractional-order PWL chaotic systems in digital platforms: (i) frequency-domain approximations (FDA) [43], (ii) numerical algorithms such as Grünwald-Letnikov (GL) [44] and Adams-Bashforth-Moulton (ABM) [45], and (iii) semi-analytical methods like the Adomian decomposition method (ADM) [46,47].Although FDA is easy to implement with electronic circuits, its results are questionable due to equilibrium points inconsistencies between the original and approximated fractional system [48,49].While ABM is highly reliable, it requires extensive simulation times and memory blocks to store plenty of solution data, making circuit implementation impractical in constrained digital platforms [44,45].Some efforts to reduce the computing time of numerical algorithms have been introduced in Refs.[50][51][52] by exploiting the short-memory concept.Since its introduction by Deng [53], the short memory principle has garnered attention in numerical computations to reduce space memory and computational cost by applying a memory window, i.e., a time interval where solution data are stored and used to compute the next iteration.For instance, the short-memory GL numerical algorithm [41] was employed for designing fractional-order chaotic systems on ARMs [54], FPGAs [40,55], and FPGA designs including coordinate rotation digital computer (CORDIC) algorithms s [23].Even so, a relatively large iteration number is still necessary since the size of the memory window has significant consequences to obtain a valid chaotic behavior [56,57].Also, the GL algorithm may not be reliable because of the sensitivity to determine the first few elements of chaotic solutions [58].
On the other hand, ADM is a suitable alternative for solving fractional differential equations because it preserves the memory effects and reduces computer time with respect to traditional numerical algorithms [46,47].The idea behind ADM consists of reformulating the system nonlinearities into the so-called Adomian polynomials to provide a solution in the form of a series expansion [59].According to Cafagna and Grassi [60], only five terms of the series are sufficient to obtain a solution with proper precision.In fact, ADM is one of the most common approaches for implementing fractional-order chaotic systems on digital hardware [40,[61][62][63], due to its time complexity O(n) and space complexity O (1).
Related works to the digital implementation of fractional order chaotic systems employing ADM are as follows.Gu et al. introduced a four-dimensional fractional chaotic system composed of quadratic and quadratic cross-product nonlinearities.Based on ADM, double-scroll chaotic attractors with offset boosting were evaluated on a Digital Signal Processing (DSP) board [64].Yang et al. presented a mathematical model of a fractional-order Hopfield neural network chaotic system that includes various hyperbolic tangent functions.Using ADM, one-scroll chaotic attractors were observed using a DSP platform [65].Liu et al. proposed a multi-attribute chaotic system with an exponential function.A numerical solution was computed using ADM to implement a double-scroll chaotic attractor on a DSP chip [66].Rajagopal et al. analyzed a fracmemristor-type chaotic system, whose memconductance function is a fourth-order polynomial.By applying ADM, two-and fourscroll chaotic attractors were verified experimentally on an FPGA [67].Xu et al. examined a fractional-order system with two exponential nonlinearities applying ADM [68].Fractionalorder Lorenz and Lorenz-Stenflo systems with a double-scroll chaotic attractor have been solved and implemented in DSP hardware using ADM [63,69,70], and so on [71,72] Based on the aforementioned literature review, we can distinguish that there are some aspects of the standard ADM formulation that have not been addressed yet, such as:

•
Piecewise linear (PWL) functions.The reported digital implementations are predominantly based on fractional-order systems with quadratic functions, cubic functions, or quadratic cross-product functions.This could be caused because determining Adomian polynomials of complex functions is challenging [73].Even when there are approaches to efficiently compute the Adomian polynomials of trigonometric, exponential, and logarithmic functions [59,74]; and improved ADM versions [75][76][77][78]; they still cannot represent PWL functions.

•
Multi-scroll chaotic attractors.The informed physical implementations of chaotic systems using ADM are mostly confined to one-and two-scroll attractors.Indeed, the classical ADM does not cope with PWL functions, which are a fundamental mechanism to get multi-scroll chaotic attractors [79][80][81][82][83]. Thus, a relevant open problem is how to incorporate fractional-order dynamical systems based on PWL functions in the sense of ADM.
In this manner, this manuscript presents a generalized piecewise linear decomposition method (PWL-DM) for implementing fractional-order multi-scroll chaotic attractors on digital hardware.Based on the Cartesian product of the finite partitions of phase space, an iterative numerical solution of fractional-order chaotic systems with PWL-type nonlinearities can be attained while preserving the advantages of standard ADM formulation in terms of computational complexity.As a consequence, the proposed PWL-DM offers a significant simplification, reducing computational complexity while preserving the complex dynamics of these systems, and hence, they can be efficiently applied for implementing fractionalorder PWL chaotic systems on ARMs and FPGAs no matter the chosen arithmetic formats, i.e., floating-and fixed-point.Thus, the PWL-DM can be considered as a multi-platform numerical algorithm.It is noteworthy to note that by the nature of PWL functions, this new decomposition technique does not require the computation of Adomian polynomials.So, it eliminates the unneeded terms, avoiding cumbersome digital designs.Particularly, the resulting implementations with this approach can be useful on real-time applications of secure communications [32,33], cryptography algorithms [27,84] and authentication schemes [85,86].Moreover, it is important to note that several studies have successfully used PWL functions to approximate nonlinear components, and despite the local linearization inherent in PWL approximations, the overall system retains its nonlinear characteristics such as chaos beahaviour [87][88][89].Building on this established approach, the proposed work can benefit the modeling and analysis of complex systems, allowing to simplify the handling of what would typically be complex nonlinear terms if they were expressed in their original forms.
To the best knowledge of the authors, this work provides, for the first time, a generalized approach to solving and implementing fractional-order multi-scroll chaotic systems in the sense of ADM.
The article is organized as follows.Section 2 introduces the mathematical foundations related to fractional calculus.Section 3 presents the proposed piecewise linear decomposition method, while Section 4 gives three specific examples.More particularly, the numerical solution of fractional-order 1D 2-scroll, 1D 4-scroll, and 2D 4 × 4-scroll chaotic systems is detailed.Section 5 demonstrates the usefulness of the proposed approach by implementing fractional-order chaotic systems using ARM and FPGA digital platforms.Finally, Sections 6 and 7 discuss the results and conclude the work, respectively.

Background
Fractional calculus generalizes the integer-order differentiation and n-fold integration.From the Riemann-Liouville theory, the fractional integral operator is given as follows.
Definition 1 ([90]).Let Γ(q) be the Gamma function with q ∈ R + .The operator J q t 0 , defined on Lebesgue space L 1 [t 0 , t 1 ] by is called the Riemann-Liouville fractional integral operator of order q with J 0 t 0 = I, the indentity operator.
For γ = 0 in the last equation of (2), we obtain the On the other hand, the classical fractional derivatives are given by Riemann-Liouville, Caputo, and Grünwald-Letnikov (GL) definitions that present global properties [1] and satisfy the Ross criterion [91].In particular, the Caputo operator clarifies the physical meaning of initial conditions for initial value problems [90].
Additionally, the Caputo derivative is the left inverse of the Riemann-Liouville integral but not its right inverse [39], as shown by: where the superscript + denotes causality.It means f (t) being a causal function of t, that is, f (t) = 0 for t < 0.

Proposed Piecewise Linear Decomposition Method (PWL-DM Approach)
Let us consider a class of fractional-order PWL systems defined by , is the constant for the autonomous system.n ∈ N + stands for the number of system variables, and U : R n → R n represents a piecewise vector which commutes as follows: where Φ = {Φ 1 , Φ 2 , . . ., Φ p 1 } is a finite partition of the phase space, which satisfies R n = p 1 ς=1 Φ ς .The vector U changes depending on which domain Φ ς ⊂ R n the trajectory is located.For that purpose, the vectors U ς with ς = 1, 2, . . ., p 1 , can be defined as U ς = u ς1 , . . ., u ςn T ∈ R n , which entries are given by: Herein, f 1 , f 2 , . . ., f n denote piecewise linear functions (PWL) whose commutation law obeys: where m i j , c i j ∈ R for i = 1, 2, . . ., φ, j = 1, 2, . . ., k i , represent the slope and intercept of every straight line in the PWL functions.φ ∈ N + indicates the number of PWL functions considered in Equation ( 8), k i ∈ N + denotes the total number of affine segments forming each function f i (x i ), and the subscript j symbolizes the element (slope or intercept) belonging to a specific PWL function given by superscript i.Hence, each function f i (x i ) generates sub-domains . ., H φ , where the surfaces S i in terms of s i j with 1 ≤ j ≤ k i − 1 specify the boundaries between two consecutive sub-domains H i j .Note that the intersections among those sub-domains compose the finite partition of the phase space Φ ς for ς = 1, 2, . . ., p 1 with p 1 = ∏ φ i=1 k i .In particular, the role of matrix B in Equation ( 6) is to define which and how many PWL functions are active at a given variable.For the general case, it is assumed that the number of PWL functions equates to the system variables, i.e., φ = n.Now, we consider the scenario with a single PWL function (φ = 1) composed by k 1 affine segments and X ∈ R 3 , then vectors U ς are represented by u ς1 , u ς2 , u ς3 T = f 1 (x 1 ), 0, 0 T .As a result, the affine linear system that governs the dynamics in the domain Φ ς , with ς = j In this manner, with some abuse of notation, the component U(X) in Equation ( 6) can be split into a piecewise vector U L (X) and piecewise constant vector U C according to the slopes m i j and intercepts c i j of PWL functions, respectively.Then, the fractional-order PWL system (6) becomes to: By applying the fractional integral operator J q t 0 and property (5) to both sides of Equation (10) with initial state we obtain: with ω = 1.Based on the Adomian decomposition method, the numerical solution of Equation ( 11) is: where ℓ = 0, 1, . . ., ∞.
Remark 1. Equation ( 12) is the exact solution of Equation ( 6), but it is important to note that in practical implementations, it is impossible to calculate infinite items of X ℓ .Nevertheless, because the Adomian decomposition method converges very fast [60,92], we only need to take the first items for the approximate solution of Equation (6).To discretize this equation, we divide the interval [t 0 , t] into subintervals [t n , t n+1 ].In this manner, the ρ−term approximate solution can be expressed as . Then we can obtain the discrete iterative form X(t n+1 ) = F X(t n ) , which is denoted as X(n + 1) = F(X(n)) with ρ = 5 for general cases.
It is important to remark that the Adomian polynomials computation is not required by the proposed PWL-DM because there are no explicit nonlinear terms in Equation (12).Now, the next step consists of deriving the piecewise-linear decomposition series from (12).The first step is to determine the finite partition Φ ς of phase space.Since the PWL functions (φ) may not necessarily equate to the system variables since φ ≤ n, and knowing that each PWL function can have an arbitrary number of affine segments k i , i.e., a different length of elements (slopes and intercepts), we establish Theorem 1.
Theorem 1.For a class of fractional-order PWL systems given by Equation (6) containing various PWL functions, i.e., φ ≥ 1, with each function having an arbitrary length k i , the finite partition Φ ς for ς = 1, 2, . . ., p 1 of the phase space with H = H 1 , H 2 , . . ., H φ = H i and H . ., φ, j = 1, 2, . . ., k i , is determined by: being H 1 × H 2 × . . .× H φ , with some abuse of notation, the Cartesian product defined as the set of all possible ordered φ-tuples.∩(•) stands for the intersection among each one of the resulting φ-tuples.Thus, R n is partitioned in p 1 sub-domains with p 1 = ∏ φ i=1 k i .
Proof.Without loss of generality, let us assume three PWL functions (φ = 3), having lengths k 1 = 2, k 2 = 3, and k 3 = 4, respectively.Thus, the number of sub-domains is p 1 = ∏ 3 i=1 k i = 24 as the cardinality of a Cartesian product of three sets is equal to the product of the cardinalities of the sets: , and 3 , H 3 4 , respectively.As a result of the p 1 combinations, the finite partition Φ of phase space contains the sub-domains: where the possible combinations can be recast as the Cartesian product of H 1 × H 2 × H 3 from all possible ordered 3-tuples.This result holds for an arbitrary number of elements k i composing the φ-PWL functions.Let H 1 , . . ., H φ be φ-non-empty sets.The Cartesian product H 1 × . . .× H φ is given by the set of all possible ordered φ-tuples (H 1 1 , ..., H i j ), where H i j ∈ H i , i = 1, . . ., φ Besides, the cardinality of the Cartesian product of φ-sets is obtained by , demonstrating that the total number of sub-domains is p 1 = ∏ k+1 i=1 k i , completing the proof.
Once the sub-domains in the finite partition of phase space have been determined, the next task consists of incorporating the slopes m i j and intercepts c i j of all PWL functions in the proposed approach.We thus propose a generalized matrix-based algorithm in Lemma 1 for a systematic determination of the decomposition series with computational purposes.Lemma 1.For the class of fractional-order PWL systems in Equation (6) satisfying Theorem 1, the components U L (X) and U C in the piecewise linear decomposition method (11) can be rearranged as a pair of matrices M and C, composed of the slopes and intercepts of PWL functions: where M i , C i , i = 1, 2, . . ., φ, are row vectors whose elements are given by: The brace symbol means that those elements will be repeated several times as the product of the k i 's indicated below them.M and C are n × p 1 matrices.
Proof.Suppose that M 1 and M 2 are non-empty sets containing the slopes m 1 j and m 2 j , with i = 1, 2, . . ., φ, j = 1, 2, . . ., k i , of two PWL functions f 1 (x 1 ) and f 2 (x 2 ), respectively.Here, k i means the number of segments in ith PWL function and, therefore, the number of slopes while φ is the number of PWL functions.Also, let us consider k 1 = 3 and k 2 = 2. Based on Theorem 1, the Cartesian product of the elements of vectors 2 )}.For the sake of simplicity, we have omitted the variables x i .To construct matrix M, the vector M i is repeated a certain number of times for each i.In this case: Thus, we have: where M 1 is rewritten by the repetition (p 3 times) of its elements: and M 2 is repeated p 2 times, with p 2 = ∏ φ−1 i=1 k i , and p 3 = ∏ φ i=2 k i .. In this manner, the final form of matrix M is: From a computational point of view, we can now straightforwardly access the ςth column of matrix M, denoted by M ( * ,ς) , ς = 1, 2, . . ., p 1 with p 1 = ∏ φ i=1 k i , which ensure that each slope and intercept match with its respective sub-domain of the finite partition of phase space.
Since a similar procedure can be developed for the intercepts C i to get the matrix C, it will not be explained here.In conclusion, we show how the components U L (X) and U C in the PWL-DM ( 11) can be rearranged into the matrices M and C by considering the slopes and intercepts for PWL functions with M i and C i with i = 1, 2, • • • φ, and an arbitrary number of segments k i .This ends the proof.
Based on Theorem 1 and Lemma 1, we can now recast the numerical solution of Equation ( 12) into a piecewise-linear decomposition series, named herein as PWL-DM, as follows: where M ( * ,ς) and C ( * ,ς) , ς = 1, 2, . . ., p 1 represent the ςth column of matrices M and C in (15).In particular, the notation ( * , ς) refers to the indexing expression to obtain all row elements in the ςth column.Additionally, Φ ς denotes the finite partition of the phase space given by the p 1 intersections in Equation ( 13), i.e., R n = ς=1 Φ ς = ∅, whereas the superscript ℓ symbolizes the element of the decomposition series.Thus, the decomposition series X 0 , X 1 , . . ., X ℓ in each sub-domain is obtained from: where the ρ−term approximate vector solution X ρ is given by X ρ = ∑ ρ−1 ℓ=0 X ℓ .An exact solution of X can be obtained when X = lim ρ→∞ X ρ .However, Cafagna and Grassi proved that a practical value ρ = 5 [60] is sufficient to get a suitable precision [60,92].

Convergence Analysis
In order to determine the convergence order and accuracy of the PWL-DM, we solve the following linear initial value problem: The analytical solution of ( 22) is given by: where E q,1 (z) is the two-parameter Mittag-Leffler function.In particular, the solution when λ = 1 and b 0 = 0.01 is plotted in Figure 2 for various fractional orders q.After that, we determine the approximated solution of Equation ( 22) by using the PWL-DM, yielding: which is compared with the exact solution in Figure 2. In the literature, we can find various possible criteria for evaluating the performance of the new algorithms.On one side, when one disposes of an analytical expression to find an exact solution, the accuracy and efficiency can be evaluated using the Mean Absolute Error (MAE) defined as [93]: with an Absolute Error (AE) AE = |x j − x N j |, x j and x N j are the exact and the approximate responses, respectively, while N means the number of interior mesh points.In addition, we can also estimate the Expected Convergence Order (ECO) as a function of the step size h by using [93]: Table 1 provides a summary of the accuracy and efficiency of PWL-DM.The findings reveal that reducing the step size h in PWL-DM yields a reduction in the mean absolute error while maintaining the expected convergence order within acceptable ranges.Both Figure 2 and Table 1 illustrate that PWL-DM offers suitable ECO and MAE values for various fractional orders and mesh refinements.
Table 1.MAE and ECO for expression (24) obtained from the proposed PWL-ADM under distinct step sizes h and fractional orders q.

Fractional Order
Step Size MAE ECO  24) and the numerical approximation with PWL-DM (25) using only two terms for distinct step sizes h and fractional orders q.Magnitude of the absolute errors (AE).

Illustrative Example
In this section, the proposed decomposition method (PWL-DM) is applied for obtaining the numerical solution of a fractional-order PWL system with known dynamic features, which has been reported in [79,80].It is worthy to remark that the PWL-DM algorithm can be applied to the general class of fractional-order PWL dynamical systems (6).

Example 1:
Step-by-Step Demonstration for a 2-scroll PWL Chaotic System Consider a fractional-order dynamical system given by Equation (6) in R 3 with linear operators: where α 31 , α 32 , α 33 , and β 31 are positive real constants.Additionally, the system possesses a single PWL function f 1 , expressed in the form of Equation ( 9) as: Therefore, the vector U in Equation ( 7) commutes among three sub-domains: while the vectors U ς in Equation ( 8) with ς = 1, 2, 3 have entries defined by: In this case, the vectors U ς can be rewritten as: This system presents a 2-scroll chaotic attractor when α 31 = 1.4,

Using PWL-DM
In this scenario, the fractional system possesses a single PWL function, i.e., φ = 1.Theorem 1 and Lemma 1 are used for getting the finite partition Φ of phase space as well as matrices M and C. We observe that f 1 (x 1 ) contains three straight lines, i.e., k 1 = 3, so the space R 3 is partitioned in three sub-domains: , and Φ 3 = H 1 3 with boundaries s 1 1 and s 1 2 .Besides, the slopes and intercepts of each straight line are expressed as . Thus, we have: A clear consequence of the above result is that the column vectors are given by: Thus, for this example, the PWL-DM (20) can be rewritten by: with p 1 = k 1 = 3, and considering that for this case M ( * ,1) = 0, M ( * ,3) = 0, C ( * ,2) = 0, G = 0, and ω = 1.Now, we are going to compute the ρ−term approximate vector solution with ρ = 5, i.e., the first five elements of the decomposition series X 0 , X 1 , X 2 , X 3 , and X 4 in each sub-domain.By exploiting the recursive relationship (21) and the PWL-DM (34) with initial conditions T , the first item X 0 in the sub-domain H 1 2 is obtained by: By defining the coefficients ∆ 0 1 = x 0 1 , ∆ 0 2 = x 0 2 , ∆ 0 3 = x 0 3 , and using the the fundamentals properties (2) of operator J q t 0 , the second item X 1 in the sub-domain H 1  2 can be obtained applying Equations ( 21) and ( 34): X 1 = J q t 0 AX 0 + J q t 0 BM ( * ,2) (X 0 ), with matrices A, B from ( 27) and vector M ( * ,2) ̸ = 0 in (32).Then, we have: with ).Thus, we obtain: being Based in the PWL-DM (34), and applying sequentially Equation ( 21), X 3 = J q t 0 AX 2 + J q t 0 BM ( * ,2) (X 2 ), X 4 = J q t 0 AX 3 + J q t 0 BM ( * ,2) (X 3 ), the fourth and fifth sets in sub-domain H 1  2 are given respectively as follows: X 4 with Once the terms X 0 , X 1 , X 2 , X 3 , X 4 have been determined, a solution for sub-domain H 1  2 can be stated.The next step consists of determining the solutions for the other partitions, i.e., H 1 1 and H 1  3 .Analogous to the procedure for the sub-domain H 1 2 , the first item of iterative Formula (21) for Equation (34) under sub-domains H 1 1 and H 1 3 is obtained applying X 0 = J q t 0 BC ( * ,1) + Υ 0 , and X 0 = J q t 0 BC ( * ,3) + Υ 0 , respectively.
X 0 We note that elements c 1 1 and c 1 3 have symmetry.Thus, we only present the results for sub-domain H 1  1 due to the partition H 1 3 will have the same form.In this manner, by defining , the second item is obtained by eq.X 1 = J q t 0 AX 0 , as follows: Then, we have the following coefficients: The third set of coefficients is computed with X 2 = J q t 0 AX 1 , resulting: Now, we obtain the coefficients: Similarly, the fourth and fifth sets X 3 = J q t 0 AX 2 , and X 4 = J q t 0 AX 3 are given by: and respectively, whereas the delta coefficients are defined as: and Finally, we can recast the solution of the fractional-order PWL dynamical system defined by Equations ( 27)-( 30), using the proposed piecewise linear decomposition method (20) as follows: According to Ref. [80], the fractional-order PWL chaotic system ( 27)-( 30) can generate a 2-scroll chaotic attractor when α 31 = 1.4, and s 1 2 = 1. Figure 3 shows the numerical simulation results obtained from the proposed PWL-DM (20) by setting step-size h = 0.01, initial conditions 0.1, 0.1, 0.1 T , and fractional-order q = 0.92.
In particular, the 2-scroll chaotic attractor in Figure 3 agrees well with the reported analytical results in Refs.[79,80], demonstrating the usefulness of the proposed decomposition method.Let us explain this statement.Tavazoei and Haeri [94] stated that a necessary condition for 2-scroll attractor existence in fractional-order systems is given by | arg(eig(A))| > qπ/2.For the fractional-order system ( 27)- (30), this condition is fulfilled for q = 0.92.Additionally, in chaotic systems, it was proved that if the system has a 2-scroll attractor, one of its equilibrium points is the saddle point of index 1, and the others are saddle points of index 2 [94].The scrolls are generated only around the saddle points of index 2, while saddle points of index 1 are responsible only for connecting scrolls.For the underlying system ( 27)-( 30) the equilibrium points are E 1 = (0, 0, 0), E 2 = (β 31 c 1 1 )/α 31 , 0, 0), 3 )/α 31 , 0, 0), which satisfies previous conditions [80].In this manner, we observe that the numerical solution obtained by the proposed PWL-DM (Figure 4) leads to the expected 2-scroll chaotic attractor where the scrolls are formed around the saddle points of index 2 (green dots) connected through the saddle point of index 1 (black point).the Newton-Gregory formula (NGF), and the Fractional Backward Differentiation Formula (FBDF), respectively.All solutions were computed with h = 0.01 and q = 0.98.The green dots denote saddle points of index 2, whereas the black point indicates a saddle point of index 1.
We also compare the numerical results obtained with the PWL-DM against the solutions given by well-known numerical algorithms, such as the Fractional Backward Differentiation Formula (FBDF), the Newton-Gregory formula (NGF), the generalization of the classical Trapezoidal rule (TR), and the predictor-corrector method of Adams-Bashforth-Moulton (PECE ABM) [95].As one can observe in Figure 4, the PWL-DM coincides with all other numerical results since all solutions converge to a similar chaotic 2-scroll attractor.It demonstrates the effectiveness of the proposed PWL-DM.Because chaotic systems are extremely sensitive to small variations in initial conditions and numerical schemes, the error cannot be determined straightforwardly.Therefore, in the previous section, we analyzed the convergence and accuracy of the PWL-DM using a linear fractional differential equation.
Finally, the trade-off between the discretization step h and the number of terms used to compute the decomposition series in (21) are numerically investigated in Table 2.We found that the proposed PWL-DM is not significantly influenced by those parameters since the expected 2-scroll chaotic attractor is obtained from three terms with diverse discretization steps.

Example 2: 1D Multiscroll PWL Chaotic System
In this example, we consider a fractional-order PWL system in the form of Equation (6) that generates a 1D 4-scroll chaotic attractor.Let us define the same linear operators as the ones given in Equations ( 27)-( 29), but the PWL function f 1 is now defined by: Now, the vector U in Equation ( 7) changes among seven sub-domains: whereas the vectors U ς in Equation ( 8) with ς = 1, 2, . . ., 7, can be defined as: By setting suitable parameters α 31 = 0.86, α 32 = 0.72, α 33 = 0.72, -scroll chaotic attractor can be obtained.

Using PWL-DM
Similar to the previous example, the fractional system possesses only a single PWL function.Based on Theorem 1 and Lemma 1, we find matrices M and C: Due to k 1 = 7, the finite partition of space is given as In this manner, the column vectors required by the decomposition (20) are: Thus, the piecewise linear decomposition method (20) is recast as: and ω = 1.We observe that the iterative formula for sub-domains H 1 4 , and can be computed as was derived in Equations ( 35)-( 39) and ( 40)-( 49), respectively.As a consequence of these findings, only iterative items for the remaining sub-domains H 1 2 and H 1 6 are needed.Because both domains are symmetrical, we only focus on the H 1 2 partition.Then, the first iterative item is attained with X 0 = J q t 0 BC ( * ,2) + Υ 0 , with B and C ( * ,2) from Equations ( 27) and ( 56) respectively: The second item X 1 is obtained applying X 1 = J q t 0 AX 0 + J q t 0 BM ( * ,2) (X 0 ), with matrices A, B from ( 27) and vector M ( * ,2) ̸ = 0 in (55).Then, we have: (59) Recursively, we obtain X 2 , X 3 , X 4 applying X 2 = J q t 0 AX 1 + J q t 0 BM ( * ,2) (X 1 ), X 3 = J q t 0 AX 2 + J q t 0 BM ( * ,2) (X 2 ), and X 4 = J q t 0 AX 3 + J q t 0 BM ( * ,2) (X 3 ), respectively.

Example 3: 2D Grid Multi-Scroll PWL Chaotic System
Consider a fractional-order dynamical system given by ( 6) in R 3 with linear operators: (70)

Physical Implementation of the PWL-DM Based on FPGA and ARM Digital Boards
As well known, fractional-order chaotic systems may enhance chaos-based engineering applications.Since cutting-edge technologies (IoT, robotics, cloud, machine learning, AI, and so on) are dominated by a digital world, we must explore the physical realization of chaotic systems on digital platforms.For fractional-order chaotic systems, there exist two standard implementation procedures [96].The first approach uses FPGA platforms configured with hardware description languages such as Verilog or VHSIC Hardware Description Language (VHDL), to describe the configuration of logic cells.The FPGA-based implementation provides a fast execution speed appropriate for real-time processing since several instructions can be executed in a single clock cycle.The second approach involves programming in C or Python for sequential computing in Digital Signal Processors (DSPs) and Advanced RISC Machines (ARMs) [96].Due to its serial mode operation, this approach is relatively slow and fails to meet real-time requirements.Nevertheless, data transmission is straightforward through built-in Ethernet or WiFi interfaces.The proposed PWL-DM can be easily adapted to various hardware.In the next subsections, we demonstrate the implementation of the fractional-order PWL chaotic systems in example-1 (Equations ( 34) and ( 50)), example-2 (Equations ( 57) and (63)), and example-3 (Equations ( 70) and ( 71)), using FPGAs and ARMs.

Non-Embedded Design: PWL-DM Implementation Based on FPGA
The FPGA-based implementation approach is realized on the Programmable Logic (PL) of the Xilinx Zynq-7000 by employing VHDL for the hardware configuration.Figures 7 and 8 illustrate the overall schematics of the digital design for the 2-scroll, 4-scroll, and 2D 4 × 4-grid scrolls chaotic attractors, respectively.For FPGA design, the arithmetic operations are performed in fixed-point Qa.b format (8 bits for integer representation plus one bit for sign and 23 bits for fractional representation).The proposed block design is divided as follows:

1.
Chaotic Oscillator Unit: The elements of the FPGA implementation for the ×2, ×4 and 4 × 4 scrolls are marked in color red.The highest level of the FPGA design for the ×2 and ×4-scrolls (see Figure 7) comprises two blocks, whereas the 4 × 4-scrolls (see Figure 8) contains three blocks, each one performs calculations separately.The PWL blocks that describe the PWL function (Equations ( 28), ( 51) and ( 65)) are shown in Figure 9. Finally, the ADM_description block executes the proposed PWL-DM (20).2.
Digital/Analog conversion: For visualization of the chaotic attractors on an oscilloscope, the Post_processing and dac_hdr blocks (blue square in Figures 7 and 8), control the 12-bit dual-channel D/A converter MCP4922 at 20 MHz.It is worth remembering that both blocks may be deleted in practical applications, i.e., they should not be considered as part of the final implementation since they are just used for visualization purposes.
A more detailed discussion on the generation mechanism of multi-scroll chaotic attractors is given as follows.

1D 2-scrolls generation mechanism
In Figure 10 is shown the 1D 2-scroll attractor generation, the process begins with the chaotic function generator f_x1 which feeds the signal into the Comparators block.This block contains several comparison units that determine whether the signal x 1 belongs to specific regions H 1 1 , H 1 2 and H 1 3 .1.
Signal Comparison.The signal x 1 is compared against the predefined regions H 1 1 , H 1 2 and H 1 3 .

2.
Multiplexer (MUX).Based on the comparison results, the MUX selects appropriate values for the parameters m j and b j .3.
PWL-DM Description.The parameters m 1 j and b 1 j are fed into the ADM_description block, which uses a Finite State Machine (FSM) that transitions from state s0 to state s5.Each state processes parallel calculations corresponding to a set of coefficients ∆ ℓ i .In state s5, the calculation results of the coefficients ∆ ℓ i are used to calculate the solution of the system and hence to generate a 2-scroll chaotic attractor.

1D 4-scrolls generation mechanism
In Figure 11 is shown the 1D 4-scroll attractor, the process is similar but involves more comparison units and states to accommodate the additional scrolls.

1.
Extended Signal Comparison.The signal x 1 is compared against more regions, from H 1 8 to H 1 7 , to determine which parameters to use.

2.
Expanded MUX.The MUX block is expanded to handle the additional regions and selects appropriate parameters m 1 j , b 1 j for generating a 4-scroll chaotic attractor.

3.
PWL-DM Description: Similar to the previous implementation, the FSM within the ADM Description block, transitions from the state s0 to the state s5, ensuring that the correct ∆ ℓ i are properly computed in each state, resulting in the generation of a 4-scrolls chaotic attractor.

2D 4 × 4-grid scrolls chaotic attractor generation mechanism
Figure 12 shows the generation of 2D 4 × 4-grid scrolls chaotic attractors, which involves a more complex mechanism as it uses two PWL functions, f_x1 and f_x2.

2.
Separate MUX Units: Each signal has its own MUX unit, which selects the parameters m 1 j , b 1 j for x 1 and m 2 j , b 2 j for x 2 .

3.
Integrated PWL-DM Description: The outputs of the two MUX units are fed into a single ADM_description block, which combines the parameters and processes them using an FSM.The FSM inside this block handles the transition from state s0 to state s5.
It ensures that the combination of x 1 and x 2 produces the desired 2D 4-scroll chaotic attractor.Consequently, the calculations in this implementation are relatively more complex compared to the previous two.
The generation of all cases of chaotic attractors is quite similar; the key difference lies in the number of regions to be evaluated, the number of PWL functions, and within the ADM_description block, the configuration and calculation of different sets of coefficients depending on each state.This configuration results in the generation of different multiscroll chaotic attractors.

Embedded Design: PWL-DM Implementation Realized on an ARM Platform
In this subsection, we present the methodology for the ARM implementation of the ×2, ×4, and 4 × 4 scroll chaotic attractors.The design process works as follows. 1.
Device configuration.The ARM-based design illustrated in Figure 13 is elaborated in the Xilinx Vivado Design Suite, version 2016.4.This block contains the ZYNQ7 Processing System IP, which is the software interface around the Processing System (PS) of the Xilinx Zynq-7000 board, currently working at 667 MHz.For including the data converter, a functional interface for communicating the PS with the PL of the Xilinx Zynq-7000 board, operating at 100 MHz, was designed (Figure 14).This interface comprises a Processor System Reset IP, an AXI Interconnect IP, and two AXI GPIO IPs.Later, the dac_hdr_v1_0 is a custom IP for controlling the MCP4922 D/A converter at 20 MHz. 2.
Calculation Preparation.From now on, the Xilinx Software Development Kit (XSDK), version 2016.4,will be used for the configuration of the PS of the Xilinx Zynq-7000 board in C programming language by using a single-precision (32-bits) floating-point format.The Algorithm 1 shows the pseudo-code for implementing the PWL-DM on an ARM digital platform.3.
Numerical solution .In this step, the PWL-DM in Equations ( 34), ( 57) and ( 70) is computed recursively using Algorithm 1 in the ARM.Finally, the chaotic time series are collected for further visualization on the Tektronix TDS3054C oscilloscope.

Discussion
For evaluation purposes, both implementations were carried out using the same platform, Zynq-7000 XC7Z020, which is a System on Chip (SoC) board constituted of a dual-core ARM Cortex-A9 Processing System (PS) and a Xilinx 7-series Programmable Logic (PL) equivalent to an Artix-7 FPGA. Figure 15a shows the experimental setup for the FPGA implementation.
Figures 15b and 16a,b present the experimental results for the FPGA realization.In particular, the phase portraits for ×2, ×4, and 4 × 4-grid scroll chaotic attractors are obtained from the proposed PWL decomposition method in Equations ( 34), ( 57) and (70), respectively.On the other hand, the experimental outcomes of PWL-DM ( 50), ( 63), (71) based on the ARM platform are presented in Figure 17a,b and Figure 18a,b, respectively.
In both physical implementations, the chaotic attractors observed experimentally are consistent with those obtained by numerical simulations in Figures 3, 5 and 6, and also with the underlying theory in Refs.[79,80].We also found that the resulting FPGA implementations are uncomplicated because the PWL-DM is only based on piecewise linear descriptions.In general, the FPGA realization of chaotic systems requires the construction of mathematical models demanding several multipliers, LUTs, and DSPs blocks to perform arithmetic operations.On the contrary, the proposed PWL-DM saves hardware by a reusable FPGA design in the form of simple blocks to perform basic operations (addition, subtraction, and two's complement multiplication).Thus, the "ADM description" block in Figures 7 and 8 computes the delta terms of the ADM-PWL efficiently.Also, the PWL functions can be implemented in parallel (see Figure 9) to optimize computation time.
Γ((k+1)q+1) .D/A conversion Output: x 1 , x 2 , . . ., x n For the case of the ARM-based implementation, the Algorithm 1 receives systems parameters, Euler Gamma values, and initial conditions.Then, the proposed PWL-ADM is evaluated recursively for each sub-domain based on the delta coefficients to derive a solution.Next, any pair of chaotic time series can be visualized in an oscilloscope through the data converter MCP4922.We recall that the pseudo-code (1) can be used for fractionalorder PWL systems with any number of variables by updating the delta terms.This is possible since the PWL-DM (20) possesses a matrix-based structure used as a numerical template, bringing flexibility to digital hardware designers.Therefore, this sub-section demonstrated that PWL-DM is a cross-platform approach for implementing a class of fractional-order PWL systems by straightforward procedures.
Besides, Table 3 summarizes the hardware resources utilization for both scenarios.We found that the three FPGA designs demand very few hardware elements, being the 4 × 4-grid scrolls chaotic attractor the more expensive, which is expected considering the complexity of the attractor.In the ARM scenario, we obtain hardware needs similar to those of the FPGA design.Since the ARM design runs the pseudo-code (1) on the PS of Xilinx Zynq-7000 board (Figure 13), the three implementations consume the same hardware.No matter the chosen digital platform, we have demonstrated that the proposed PWL-DM permits a feasible physical implementation.
Next, the throughput of the FPGA implementations has also been determined.Table 4 shows the bits per second rate for the 2-scrolls, 4-scrolls, and 4 × 4-grid scrolls chaotic attractors.As one can observe, the clock frequency and throughput are slightly less when the number of scrolls increases.Table 5 presents the running computation comparison among GL, short memory GL, and PWL-DM for calculating the numerical solution of the PWL multi-scroll chaotic attractor in a dual-core ARM Cortex-A9 Processing System.According to Ref. [40], the dominant approach for a digital implementation of fractional-order multi-scroll chaotic systems is the GL algorithm and its optimized version known as short-memory GL due to standard ADM cannot describe PWL functions.Therefore, the proposed PWL-DM has been compared against classical numerical algorithms such as ABM and short-memory GL algorithms.
In particular, we have implemented the 4-scrolls chaotic system of Section 4.2 on an ARM platform.For ABM and short-memory GL algorithms, we based on the methods presented in [45,54,97], p. 19, respectively, while we use the design guidelines in Section 5.2 for the PWL-DM approach.We found the computation time is improved significantly by PWL-DM.First, the ABM algorithm consumes a lot of computation time, which is impractical for digital implementation.
On the other hand, the short-memory GL requires 96 µs per iteration, whereas the PWL-DM needs just 6.6 µs per iteration.Hence, we obtain around a 15× speed increase.It means that short-memory GL required 12.91 s to execute 200,000 iterations by only 0.88 s of the PWL-DM.Moreover, the short-memory GL algorithm may suffer from a complicated determination of the memory window [54,57].Additionally, we have also compared our approach against the short-memory GL algorithm, which is one of the most dominant approaches to implementing fractionalorder chaotic systems in FPGAs [40] First, it is important to remark that for FPGA-based implementations, to the best knowledge of the authors, there are no works showing the implementation of fractional-order PWL chaotic systems using the standard ADM, which is expected since classical ADM cannot express PWL functions.In fact, that is the goal of our manuscript.Thus, we present the comparison in two tables since the hardware resource utilization depends on the FPGA chip and digital synthesis methodology employed for implementation.Table 6 shows the comparison of hardware resource utilization between the proposed PWL-DM and short-memory GL implementation of the PWL chaotic system with 2-scrolls shown in Example-1, Equations ( 27) and (28).In this manner, Table 6 presents a direct comparison since both FPGA-based designs are implemented on the same FGPA i.e., an Artix-7 FPGA.Table 7 compares the hardware resource utilization among proposed PWL-DM and other short-memory GL-based implementations for the case of fractionalorder PWL chaotic systems exhibiting 1D n-scrolls chaotic attractors.In particular, we select only those previous works related to fractional-order 1D multi-scroll chaotic systems based on PWL functions.From those results, we can observe that the proposed PWL-DM is a suitable approach for solving fractional-order PWL systems generating multi-scroll chaotic attractors with feasible implementations that perform better than classical numerical algorithms, yielding a flexible and cross-platform algorithm.

Conclusions
A generalized piecewise-linear decomposition method (PWL-DM) has been presented for computing the numerical solution of a class of fractional-order PWL dynamical systems exhibiting multi-scroll attractors.We demonstrated that the proposed decomposition method can be applied to fractional-order PWL systems composed of nth PWL functions, where each PWL function may have any number of affine segments.It has also been shown that the proposed technique depends only on linear operators and specific constant vectors for each sub-domain of the PWL function.In addition, we found that the finite partitions of phase space obey a combinatoric law given by their Cartesian product, providing a structured matrix-based algorithm.As a consequence of those findings, the PWL-DM approach can be easily implemented in ARM processors and FPGA technology since it is an algorithm flexible, cross-platform, and with a relatively low computational cost.Several numerical results and physical implementations of fractional-order multi-scroll chaotic systems demonstrated the usefulness of the proposed approach.
Therefore, the PWL-DM method can be viewed as an essential tool for broadening the engineering applications of PWL systems and should be explored further by those interested in this field.For instance, the proposed decomposition method may be the basis for developing lightweight data encryption schemes and Hopfield neural networks based on fractional-order multi-scroll chaotic systems.

Figure 1 .
Figure 1.Cartesian product of the sets M 1 and M 2 and its representation as matrix M.

Figure 2 .
Figure 2. Comparison of the exact expression(24) and the numerical approximation with PWL-DM (25) using only two terms for distinct step sizes h and fractional orders q.Magnitude of the absolute errors (AE).

Figure 4 .
Figure 4. Numerical results of the proposed PWL-DM versus the predictor-corrector method of Adams-Bashforth-Moulton (PECE ABM), the generalization of the classical Trapezoidal rule (TR), the Newton-Gregory formula (NGF), and the Fractional Backward Differentiation Formula (FBDF), respectively.All solutions were computed with h = 0.01 and q = 0.98.The green dots denote saddle points of index 2, whereas the black point indicates a saddle point of index 1.

Figure 8 .
Figure 8. RTL schematics of the FPGA-based implementation of the proposed PWL-DM for example-3 (70).

Figure 14 .
Figure 14.ARM design for communicating the Zynq PS with the PL required by the PWL-DM implementation.

Figure 15 .
Figure 15.(a) Experimental setup of the FPGA-based design given in Figures 7-9 using a Zynq-7000 XC7Z020 SoC digital platform.(b) Experimental result of the FPGA-based implementation (Figure15a) for a 2-scroll chaotic attractor with fractional-order q = 0.92 applying the PWL-DM(50).

Table 2 .
(21)rical solutions of the PWL-DM for fractional-order system (27)−(30) under various discretization steps h and a different number of terms in the decomposition series(21).

Table 3 .
Hardware utilization for FPGA-and ARM-based implementations of proposed PWL-DM.

Table 4 .
Clock frequency and throughput for the FPGA-based implementation of the proposed PWL-DM.

Table 5 .
Running time for computing the solution of fractional-order multi-scroll PWL chaotic system of Section 4.2 in a dual-core ARM Cortex-A9 Processing System using Adams-Bashforth-Moulton (ABM) algorithm, short-memory Grünwald-Letnikov (GL), and the proposed PWL-DM approach.

Table 6 .
(28)arison of the hardware resource utilization between proposed PWL-DM and short-memory GL implementation of the PWL chaotic system with 2-scrolls shown in Example-1, Equations (27) and(28).Both FPGA-based designs are implemented on an Artix-7 FPGA.

Table 7 .
Comparison of the hardware resource utilization among proposed PWL-DM and other implementations for the case of fractional-order PWL chaotic systems exhibiting 1D n-scrolls chaotic attractors.