On the Solutions of a Class of Discrete PWC Systems Modeled with Caputo-Type Delta Fractional Difference Equations

Abstract

In this paper, it is shown that a class of discrete Piece Wise Continuous (PWC) systems with Caputo-type delta fractional difference may not have solutions. To overcome this obstacle, the discontinuous problem is restarted as a continuous fractional problem. First, the single-valued PWC problem is transformed into a set-valued one via Filippov’s theory, after which Cellina’s theorem allows the restart of the problem into a single-valued continuous one. A numerical example is proposed and analyzed.

1. Introduction

PWC real-valued functions $f:D\subseteq \mathbb{R}\times {\mathbb{R}}^n\to {\mathbb{R}}^n$ are time-continuous but discontinuous with respect to the state variable, x, are defined in a finite domain D of an $(n+1)$-dimensional $(t,x)$ space, where D consists of a finite number of domains, $D_i$, $i=1,2,\dots ,k$, in each of which f is continuous up to the boundary of the domains. We denote by $\mathcal{M}$ the discontinuity set containing the boundary points. The considered discontinuity is of the jump type when in the points of $\mathcal{M}$, the function jumps (switches) and left-hand and right-hand limits exist and are different. Inside the domains, f is continuous [1].

Throughout the paper, discontinuity is considered only with respect to the state variable.

Dynamical systems modeled by this kind of PWC functions appear in many different branches of engineering and applied sciences, such as dry friction, impacting machines, systems oscillating due to earthquakes, impacts in mechanical devices, power circuits, forced vibrations, elasto-plasticity, switching in electronic circuits, uncertain systems and many others (see, e.g., [2,3,4,5,6] and their references). The vast majority of such systems are defined by time-continuous Initial Value Problems (IVPs), modeled by ODEs of integer order. The numerical integration of such IVPs is a difficult task for which only special difference methods can be used (see, e.g., [4]). While the standard methods for continuous systems rely heavily on linearization, in the case of PWC systems modeled by ODEs, they do not require linearization in general. Another difficulty is that the underlying IVP might not even admit solutions (see, e.g., [7]). Further, the PWC systems could have trajectories colliding with the discontinuity surfaces, thereby generating a new kind of bifurcation [8,9,10]. To overcome the problem of having no solutions, tools of differential inclusions of integer order can provide a possible resolution (see, e.g., the method proposed in [11,12,13,14]), where the discontinuous single-valued IVP is transformed to a set-valued IVP. Next, to obtain a numerical solution, either special numerical schemes for differential inclusions [15,16,17,18] can be utilized or, via the selection theory, continuous or even smooth approximations in the neighborhood of discontinuity can be adopted [19,20,21].

On the other side, the main existing definitions of fractional order derivatives are based on the formulae presented by Caputo, Riemann–Liouville and Grünwald–Letnikov [22,23,24,25]. However, the number of proposed definitions based on these derivatives became huge such that it represents an obstacle to the diffusion of fractional calculus in Science and Engineering [26]. For practical applications, due to the considerable advantage of allowing the coupling of differential equations with classical initial conditions as for differential equations of integer order, compared to other non-integer derivatives, the Caputo derivative is one of the most commonly used derivatives for solving fractional differential equations.

Equally interestingly, PWC systems modeled by Fractional Differential Equations (FDEs) can be approached numerically via fractional differential inclusions [27]. Once the set-valued IVP is transformed into a single-valued IVP, it can be numerically integrated using one of the existing schemes for FDEs, such as the Adams–Bashforth–Moulton method [28] (for fractional differential inclusions, see [29,30,31,32,33] and references therein).

In recent years, discrete fractional calculus has gained considerable interest, and now the study of ordinary difference equations is widespread. However, the theory of fractional difference equations, a very new area for scientists, is still evolving [34]. The left and right Caputo fractional sums and differences, as well as their properties with relation to Riemann–Liouville differences, are studied in [34] (an early paper on the theory of fractional finite difference equations), initial value problems in discrete fractional calculus are analyzed in [34,35,36,37,38,39], with existence results for nonlinear fractional difference equations presented in [34,36,38,40,41,42]. For further reading of qualitative properties of fractional difference equations, see [35,43,44,45], and for applications of discrete fractional calculus, see [46,47].

Compared to PWC systems modeled by FDEs, which represent the subject of several works, there are no results on discrete PWC systems modeled by fractional differences. Therefore, we are motivated to propose a new class of fractional discrete PWC systems modeled by Caputo delta differences. The existence of the solutions of the underlying IVPs is also studied. Moreover, because the continuity is considered as a required property for the existence of the solutions of fractional difference equations (see, e.g., [36,48]), in this paper, the continuous approximation of the PWC function is proposed. One of the significant advantages of the considered Caputo fractional difference operator over the other fractional difference operators is that it includes traditional initial and boundary conditions in formulating the problem. In addition, another advantage is the fact that the Caputo fractional difference for a constant is zero. Because the systems modeled by the considered class of discrete PWC problems might have no solutions, continuous approximation is proposed.

This paper is structured as follows: Section 2 presents the class of fractional discrete PWC systems, modeled by Caputo-type delta fractional difference, and the existence of the solutions. Section 3 deals with the approximation of the PWC function. While in Section 4, some representative numerical simulations underline the theoretical results. In the end, we give our conclusions.

2. PWC Systems Modeled with Caputo-Type Delta Fractional Difference Equations

In this paper, the considered systems are time-independent, i.e., autonomous systems.

Let us first consider a PWC system modeled by the following FDE [11]

$$\left\{\begin{array}{c}{D}_{*}^{q}x=2-3sgn\left(x\right),\hfill \\ x\left(0\right)=x_0,\hfill \end{array}\right.$$

(1)

where $D_{*}^q$ is Caputo’s derivative with starting point 0, $q\in (0,1)$, and the right-hand side is a jumping PWC function.

Proposition 1. 

The IVP (1) has no classical solutions.

Proof. 

In this case, $D_1:=D_{-}=(-\infty ,0)$ and $D_2:=D_{+}=(0,\infty )$ and $M=\left\{0\right\}$. Since $D_{*}^q\left(0\right)=0\ne 2=2-3sgn\left(0\right)$, $x\left(t\right)=0$ is not a solution. Therefore, there are no solutions starting from $x\left(0\right)=0$. Further, if one chooses $x_0\in D_{+}$, there exists a solution $x\left(t\right)=x_0-t^q/\mathrm{\Gamma }(1+q)$, but this is defined only on $[0,T^{\prime })$ with $T^{\prime }={\left(x_0\mathrm{\Gamma }(1+q)\right)}^{1/q}$ and, because it tends to the line $x=0$, cannot be extended onto larger intervals larger than $[0,T^{\prime })$. Similarly, for $x_0\in D_{-}$, the solution $x\left(t\right)=x_0+5t^q/\mathrm{\Gamma }(q+1)$ exists but, again, only a finite interval $[0,T^{″})$, with $T^{″}={(x_0\mathrm{\Gamma }(1+q)/5)}^{1/q}$. In both cases, the obtained solutions tend to the line $x=0$ but, as seen above, they hold at points $A(T^{\prime },0)$ and $B(T^{″},0)$, respectively, and cannot extend along this line (see Figure 1a, where $q=0.8$ and $x_0^{\prime }=0.2$, $x_0^{″}=-0.3$). □

Note that running this equation for some numerical scheme (such as an Adams–Bashforth–Moulton scheme for FDEs [28]), one can obtain a numerical result, but due to Proposition 1, this does not represent the correct numerical solution. For example, for $q=0.8$, Figure 1b presents the “solution” for $x_0=0$. Due to the finite precision in which computers perform calculations (see Section 2.1), the utilized numerical method can pass through points A and B, i.e., at these points, $x\left(t\right)$ is not (exactly) zero and, therefore, one obtains a wrong solution.

To overcome this obstacle, the problem has to be restarted as a differential inclusion and next as a continuous problem, which admits solutions (see details in [11,12,13,14]).

To introduce the class of discrete PWC fractional systems, some basic notions necessary to introduce the class of fractional PWC systems are presented next.

Denote by ${\mathbb{N}}_c=\{c,c+1,c+2,\dots \}$ and ${\mathbb{N}}_c^d=\{c,c+1,c+2,\dots ,d\}$, for any real numbers c and d such that $d-c\in {\mathbb{N}}_1$.

Definition 1. 

The Euler gamma function is defined by

$$\mathrm{\Gamma }\left(z\right)={\int }_0^{\infty }{e}^{-t}{t}^{z-1}dt,\Re \left(z\right)>0.$$

Using its reduction formula, the Euler gamma function can also be extended to the half-plane $\Re \left(z\right)<0$ except for $z\in \{\cdots ,-2,-1,0\}$.

Definition 2. 

Assume $u:{\mathbb{N}}_a^b\to \mathbb{R}$ and $N\in {\mathbb{N}}_1$. The first-order forward difference of u is defined by

$$\mathrm{\Delta }u\left(t\right)=u(t+1)-u\left(t\right),t\in {\mathbb{N}}_a^{b-1},$$

and the $N^{th}$-order forward difference of u is defined recursively by

$${\mathrm{\Delta }}^{N}u\left(t\right)=\mathrm{\Delta }\left({\mathrm{\Delta }}^{N-1}u\left(t\right)\right),t\in {\mathbb{N}}_a^{b-N}.$$

Finally, ${\mathrm{\Delta }}^0$ denotes the identity operator.

Definition 3. 

Let $u:{\mathbb{N}}_a\to \mathbb{R}$ and $q>0$. The $q^{\mathrm{th}}$-order delta fractional sum of u based on a is given by

$${\mathrm{\Delta }}_a^{-q}u\left(t\right)=\frac{1}{\mathrm{\Gamma }\left(q\right)}\sum _{s=a}^{t-q}\frac{\mathrm{\Gamma }(t-s)}{\mathrm{\Gamma }(t-s-q+1)}u\left(s\right),t\in {\mathbb{N}}_{a+q}.$$

Definition 4. 

Let $u:{\mathbb{N}}_a\to \mathbb{R}$, $q>0$, and $q\notin {\mathbb{N}}_1$. The $q^{\mathrm{th}}$-order Caputo delta fractional difference of u based on a is given by

$${\mathrm{\Delta }}_{a*}^{q}u\left(t\right)={\mathrm{\Delta }}_a^{-(N-q)}\left({\mathrm{\Delta }}^{N}u\left(t\right)\right),t\in {\mathbb{N}}_{a+N-q},$$

where $N=\left[q\right]+1$. If $q=N\in {\mathbb{N}}_1$, then

$${\mathrm{\Delta }}_{a*}^{q}u\left(t\right)={\mathrm{\Delta }}^{N}u\left(t\right),t\in {\mathbb{N}}_a.$$

Consider now the class of fractional PWC systems with Caputo-type delta fractional difference, modeled by the following IVP

$$\left\{\begin{array}{c}{\mathrm{\Delta }}_{*}^{q}x\left(n\right)=f\left(x(n+q-1)\right),n\in {\mathbb{N}}_{1-q},\hfill \\ x\left(0\right)=x_0,\hfill \end{array}\right.$$

(2)

where $N_{1-q}=\{1-q,2-q,3-q,\cdots \}$, ${\mathrm{\Delta }}_{*}^q$ represents the qth fractional Caputo-like difference in the usual case of a zero starting point $a=0$, with $q\in (0,1)$ [34] and f is a jump discontinuous scalar function of the following form

$$f\left(x\right)=\{\begin{array}{c}{f}_1\left(x\right),x\in (-\infty ,a],\hfill \\ f_2\left(x\right),x\in (a,\infty ).\hfill \end{array}$$

(3)

Function $f_{1,2}$ is continuous in its domain, with $f_1\left(a\right)\ne f_2\left(a\right)$.

If the solution of the fractional IVP (2) exists, it can be found with the following integral [34,41] (see [38] for ${\nabla }^q$ difference equations)

$$x\left(n\right)=x\left(0\right)+\frac{1}{\mathrm{\Gamma }\left(q\right)}\sum _{r=1-q}^{n-q}\frac{\mathrm{\Gamma }(n-r)}{\mathrm{\Gamma }(n-r-q+1)}f\left(x(r+q-1)\right),n\in {\mathbb{N}}_0.$$

(4)

To obtain a convenable numerical form, consider in (4) the following substitution $r+q=s$. Then, a convenient iterative numerical form of the sum of Equation (4) is given by

$$x\left(n\right)=x\left(0\right)+\frac{1}{\mathrm{\Gamma }\left(q\right)}\sum _{s=1}^n\frac{\mathrm{\Gamma }(n-s+q)}{\mathrm{\Gamma }(n-s+1)}f\left(x(s-1)\right),n\in {\mathbb{N}}_0.$$

(5)

The example considered in this paper is a fractional order variant of the model presented in [49], with the right-hand side

$$f\left(x\right)=\{\begin{array}{c}m-p{x}^2,x\in (-\infty ,0],\hfill \\ 1-p{x}^2,x\in (0,\infty ),\hfill \end{array}$$

(6)

with $D_1:=D_{-}=(-\infty ,0]$ and $D_2:=D_{+}=(0,\infty )$, p is a real parameter, $m\in (0,1)$ and $\mathcal{M}=\left\{0\right\}$. For all considered m values, function f has a jump discontinuity at $x=0$, $f_1\left(0\right)=m\ne 1=f_2\left(0\right)$.

The PWC (2) becomes

$$\begin{array}{c}{\mathrm{\Delta }}_{*}^{q}x\left(n\right)=\left\{\begin{array}{c}m-px{(n+q-1)}^2,x(n+q-1)\in (-\infty ,0],\hfill \\ 1-px{(n+q-1)}^2,x(n+q-1)\in (0,\infty ),\hfill \end{array}\right.\hfill \\ x\left(0\right)=x_0,n\in {\mathbb{N}}_{1-q}.\hfill \end{array}$$

(7)

Like in the case of the time-continuous PWC system (1), if one considers $x_0=0\in D_{-}$, one can see that, in this case, Equation (2) is not verified because ${\mathrm{\Delta }}_{*}^q\left(0\right)=0\ne 1=1-p\times 0$. For other values of $x_0\ne 0$, it is possible that after some iterations, the solution reaches the line $x=0$, which cannot be the solution (see the case of system (1)).

The same situation can happen in the case of the general IVP (2) with f given by (3): it is possible that, for some $x_0$, the orbit crosses the line $x=a$, which does not verify the equation. Therefore, one can be deduced that the IVP (2) with f given by (3) might have no solutions.

Remark 1. 

It is possible that, for some set of parameters and $x_0$ and q, the orbit does not cross line $x=a$ and remains in the same domain of $x_0$ (either $D_{-}$ or $D_{+}$) when the IVP admits solutions (see the example in Figure 4d, Section 4).

2.1. Computational Approach

Theoretically, it has been shown that it is possible that the solution to IVP (2) with f given by (3) can reach line $x=0$, i.e., $x\left(n\right)$ becomes 0, where the problem has no solution. However, a numerical method of calculation is an approximation that can be stable (meaning that it tends to reduce rounding errors) or unstable (meaning that rounding errors are magnified); therefore, very often, there are both stable and unstable solutions for a problem [50]. Further, in computer hardware, a value is not necessarily exactly computed, and the loss in precision could sometimes be inevitable. Moreover, considering any numeric representation that is limited to finite precision, for example, operating a decimal at 100,000,000 digits, which will be able to distinguish between the values that differ in their hundred millionth decimal places, there would still be an infinite number of values that would not be able to be exactly represented. In other words, considering the Pigeon Hole Principle, or Dirichlet drawer principle, a simple yet powerful idea in mathematics, which says that if you have n items to put into m containers where $n>m$, then at least one container must have more than one item [51,52,53]. Hence, if you have an N-bit representation of numbers, then at most $2N$ different numbers can be represented, and other numbers cannot exactly exist within that system. Therefore, it is easy to understand that the numerical schemes, such as integral (5), will not precisely meet the zero value, where $x\left(n\right)=0$ cannot be a solution. Namely it will either exceed the discontinuity or return to a previous value.

Remark 2.

(i)

Similar situations can arise in integer-order discrete systems, such as the system in [49], defined by the following IVP

$$\left\{\begin{array}{c}x(n+1)=f\left(x\right(n\left)\right),n\in \mathbb{N},\hfill \\ x\left(0\right)=x_0,\hfill \end{array}\right.$$

where f is some PWC function defined by (3);

(ii)

There are PWC systems with jump discontinuity, for which f is not defined at $x=a$ (see, e.g., [54]). In these cases, after some number of iterations, $n=k$, in the internal representation, as shown above, it is possible for $x\left(k\right)$ to enter a sufficiently small neighborhood of a, where $x\left(k\right)$ cannot be determined, and the software considers an unpredictable value for $x\left(k\right)$.

Concluding, to overcome this inconvenience, the continuous (even smooth) approximation of the PWC problem is proposed so that the underlying problem admits a solution.

3. Continuous Approximation of Map $\mathit{f}$

In this section, it is briefly shown how the PWC function, f, defined by (3), can be continuously approximated using Filippov’s approach [1,19,20] (details on the approximation algorithm can be found in [11,14]).

The PWC function f is transformed into a set-valued map $F:\mathbb{R}\rightrightarrows \mathbb{R}$ via the Filippov regularization [1], which is a map from $\mathbb{R}$ to the set of subsets of $\mathbb{R}$. F can be defined in several ways. One of the simplest forms of F is defined as follows

$$F\left(x\right)=\bigcap _{\epsilon >0}\bigcap _{\mu \left(\mathcal{M}\right)=0}\overline{conv}\left(f\right(y\in \mathbb{R}\setminus \left\{0\right\}:|y-x|\le \epsilon \left)\right),$$

(8)

where $\epsilon$ is the radius of the ball centered on x. At the points where f is continuous, $F\left(x\right)$ consists of one single point, i.e., $F\left(x\right)=\left\{f\right(x\left)\right\}$, while at the points $x\in \mathcal{M}$, $F\left(x\right)$ is given by (8). The set-valued function F defined by (8) has values in the convex subsets of $\mathbb{R}$.

To justify the use of the Filippov regularization in physical systems, the value of $\epsilon$ must be chosen to be small enough so that the motion of the physical systems approaches a certain solution (ideally, it coincides with the solution if $\epsilon \to 0$).

In the sketch in Figure 2a, the graph of a set-valued function F is plotted, while in Figure 2b, the closure of the convex hull is plotted in blue. The values of $F\left(x\right)$ for $x=x_1$ and $x=x_3$ are segments, while at $x=x_2$, $F\left(x_2\right)$ is a single point (see Condition $\gamma$ in [1] p. 68).

For function f given by (6), with the discontinuity set $\mathcal{M}=\left\{0\right\}$ and for $m=0.6$ and $p=1.5$, the graph is presented in Figure 3a. Consider a $\epsilon$-neighborhood of $x=0$. For clarity of the graphical exposition, the ray of the neighborhood is considered as $\epsilon =0.1$ (Figure 3b). The set-valued map $F:\mathbb{R}\rightrightarrows \mathbb{R}$ defined with (8) is

$$F\left(x\right)=\left\{\begin{array}{c}{f}_1\left(x\right),x<0,\hfill \\ [A\prime ,B\prime ],x=0,\hfill \\ f_2\left(x\right),x>0,\hfill \end{array}\right.$$

(9)

where $A^{\prime }$ and $B^{\prime }$ are the endpoints of the vertical segment at $x=0$. Points A and B are the intersections of the graph of f with the lines $x=-\epsilon$ and $x=\epsilon$.

With the Filippov regularization, the fractional discrete difference (2) is restarted as a set-valued fractional discrete difference (inclusion)

$$\left\{\begin{array}{c}{\mathrm{\Delta }}_{*}^{q}x\left(n\right)\in F\left(x(n+q-1)\right),\mathrm{for}\mathrm{almost}\mathrm{all}n\in {\mathbb{N}}_{1-q},\hfill \\ x\left(0\right)=x_0,\hfill \end{array}\right.$$

(10)

which is identical to (2) for those values of x for which $F\left(x\right)=\left\{f\right(x\left)\right\}$. For $x=0$, points A and B (Figure 3b) become $A^{\prime }(0,m)$ and $B^{\prime }(0,1)$, respectively, and, for $x=0$, system (7) transforms into the fractional differential inclusion, ${\mathrm{\Delta }}_{*}^q\left(0\right)\in [m,1]$, i.e., ${\mathrm{\Delta }}_{*}^q\left(0\right)$ could take every value within the line $[m,1]$.

Solutions to the set-valued IVP (10) (absolutely continuous functions satisfying (10) for almost all $n\in {\mathbb{N}}_{1-q}$) are not considered here (see, e.g., [1]).

Definition 5. 

A single-valued function $h:\mathbb{R}\to \mathbb{R}$ is called the approximation (selection) of the set-valued function F if $h\left(x\right)\in F\left(x\right)$, for all $x\in \mathbb{R}$.

Definition 6. 

As set-valued function $F:\mathbb{R}\rightrightarrows \mathbb{R}$ is upper semicontinuous at $x_0\in \mathbb{R}$, if for any open set B containing $F\left(x_0\right)$, there exists a neighborhood A of $x_0$ such that $F\left(A\right)\in B$.

It is said that F is upper semicontinuous if it is so at every $x_0\in \mathbb{R}$.

Remark 3. 

A set-valued function satisfies a property if and only if its graph satisfies it (i.e., symmetric interpretation of a set-valued function as a graph [19]). Therefore, a set-valued function is said to be closed if and only if its graph is closed. Further, a set-valued function $F:\mathbb{R}\rightrightarrows \mathbb{R}$ whose graph is closed is upper semicontinuous [19] p. 42.

Finding the approximations, which are locally Lipschitz, is allowed by the Approximate Selection Theorem (Cellina’s Theorem), whose proof presents an explicit way (see [19] p. 84 and [20] p. 358) to construct the approximation.

It is easy to see that the function defined by (9) has a closed graph and, therefore, it is upper semicontinuous and admits a continuous (even smooth) selection g (Figure 2b, red plot).

Theorem 1 

([21] (see also [19])). Let $F:X\rightrightarrows Y$ be an upper semicontinuous function from a compact metric space X to a Banach space Y. If the values of F are nonempty and convex, then for every $\epsilon >0$ there exists a locally Lipschitz single-valued map $g:X\to Y$ such that

$$Graph\left(g\right)\subset B\left(Graph\right(F),\epsilon ),$$

and for every $x\in X$, $g\left(x\right)$ belongs to the convex hull of the image of F.

Next, the main result can be presented

Theorem 2. 

The set-valued map $F:\mathbb{R}\rightrightarrows \mathbb{R}$, defined by (9), admits a locally Lipschitz selection $g:[-\epsilon ,\epsilon ]\to \mathbb{R}$.

Proof. 

From Remark 3 it follows that F defined by (9) is upper semicontinuous, nonempty and convex. Therefore, Theorem 2 applies. □

One of the simplest continuous approximations of the discontinuous function f (6) is the cubic polynomial (There exists an infinity of smooth functions to approximate the PWC function f).

$$g\left(x\right)=c_1{x}^3+c_2{x}^2+c_{3}x+c_4,c_i,\in \mathbb{R},i=1,2,3,4.$$

The approximated function (2) becomes

$$\tilde{f}\left(x\right)=\left\{\begin{array}{c}m-p{x}^2,x<-\epsilon ,\hfill \\ g\left(x\right),x\in [-\epsilon ,\epsilon ],\hfill \\ 1-p{x}^2,x>\epsilon .\hfill \end{array}\right.$$

Since $f_1$ and $f_2$ in (6) are smooth, to define g as connecting points A and B, the following “gluing” conditions are to be set

$$\begin{array}{c}\tilde{f}(-\epsilon )=g(-\epsilon ),\hfill \\ \tilde{f}\left(\epsilon \right)=g\left(\epsilon \right),\hfill \\ {\tilde{f}}^{{}^{\prime }}(-\epsilon -0)=g^{{}^{\prime }}(-\epsilon +0),\hfill \\ {\tilde{f}}^{{}^{\prime }}(\epsilon +0)=g^{{}^{\prime }}(\epsilon +0),\hfill \end{array}$$

which represents a system with unknown $c_i$, $i=1,2,3,4$. ${\tilde{f}}^{{}^{\prime }}(\pm \epsilon \pm 0)$ and $g^{{}^{\prime }}(\pm \epsilon \pm 0)$ are lateral limits of the derivatives $\tilde{f}$ and g at $\pm \epsilon$. Note that the last two equations represent the smoothness conditions on points A and B.

Solving the system, one obtains

$$\begin{array}{cc}\hfill c_1& =\frac{m-1}{4{\epsilon }^3},\hfill \\ \hfill c_2& =-p,\hfill \\ \hfill c_3& =-\frac{3m-3}{4\epsilon },\hfill \\ \hfill c_4& =\frac{m+1}{2}.\hfill \end{array}$$

Finally, the obtained smooth discrete fractional system is

$$\begin{array}{c}{\mathrm{\Delta }}_{*}^{q}x\left(n\right)=\tilde{f}\left(x\left(n\right)\right):=\left\{\begin{array}{c}m-px{(n+q-1)}^2,x\left(n\right)<-\epsilon ,\hfill \\ g\left(x\right(n+q-1\left)\right),x\left(n\right)\in [-\epsilon ,\epsilon ],\hfill \\ 1-px{(n+q-1)}^2,x\left(n\right)>\epsilon ,\hfill \end{array}\right.\hfill \\ x\left(0\right)=x_0,n\in {\mathbb{N}}_{1-q}.\hfill \end{array}$$

(11)

The existence of the following numerical integral

$$x\left(n\right)=x\left(0\right)+\frac{1}{\mathrm{\Gamma }\left(q\right)}\sum _{s=1}^n\frac{\mathrm{\Gamma }(n-s+q)}{\mathrm{\Gamma }(n-s+1)}\tilde{f}\left(x(s-1)\right),n\in {\mathbb{N}}_0,$$

(12)

is ensured by the smoothness of the right-hand side of (11) [41,42,55].

To computationally implement integral (12), the entire orbit history (main characteristic of fractional systems) must be taken into account. Therefore, a modality is inside the cycle that calculates the sum, every step $x(s-1)$ is tested for which the domain belongs.

4. Dynamics of the Approximated Fractional System (11)

Before studying the dynamics of fractional system (11), recall the following important result regarding continuous and discrete fractional systems [56,57].

Theorem 3. 

Autonomous, continuous-time and discrete fractional systems cannot admit non-constant exact periodic solutions.

Proof. 

This result regarding continuous fractional systems modeled by fractional order differential equations is proven in [56], while for discrete fractional systems, it is proven in [57]. □

Remark 4. 

Due to Theorem 3, periodicity cannot be considered in continuous or discrete fractional systems. Therefore, notions of stable cycles, bifurcation and even chaos (where unstable periodic orbits form the skeleton of chaotic dynamics) represent a delicate problem. Thus, following the definition given by, e.g., Wiggins in [58]: a non-constant solution $x\left(t\right)$ of a system is periodic if there exists $T>0$ such that $x\left(t\right)=x(t+T)$, for all $t\in \mathbb{R}$, it follows that even using some asymptotic approach, one cannot obtain periodic orbits in fractional systems. Instead, one can consider numerically periodic orbits in the sense that the trajectory, from the numerical point of view, up to some small error, can be considered in the state phase as a closed orbit. However, there are particular cases when one can talk about periodicity in the case of continuous fractional systems when the lower terminal of the fractional derivative is $\pm \infty$ (see, e.g., [59]). Further, in the case of discrete fractional systems, there could exist S-asymptotically periodic orbits [57].

To obtain the numerical results in this section, a Matlab code has been written.

Consider first the not approximated system (7) with parameters $m=0.92$, $p=1.556$ and $q=0.6$. Applying the integral, one obtains the chaotic orbit presented in Figure 4a. The value of the orbits close to 0 is plotted in red. As can be seen, integral (5) gives a numerical result (orbit), and $x\left(n\right)$ is close to 0 (see Section 2.1) at the $n_0\approx 1200$th iteration. However, as shown in Section 2, the result is not correct.

If one considers the approximated fractional system (11), with $\epsilon ={10}^{-3}$ and the same parameters $m=0.92$, $p=1.556$ and $q=0.6$, the obtained correct chaotic orbit is presented in Figure 4b. Note that the approximation is performed only in a relatively large neighborhood of the discontinuity. As can be seen in the images, as expected, the differences between the non-approximated and approximated cases appear only after the intersection with the line $x=0$ and within the neighborhood of the ray $\epsilon$, respectively (see the vertical dotted red line at $n_0$ in Figure 4a,b). Because, in the approximated case where the code also considers the function g, once the orbit enters the neighborhood after $n>n_0$, the orbits are different.

A numerically periodic orbit (Remark 4) can be obtained for $q=0.6$, $p=1.2$ and $m=0.92$, with $\epsilon ={10}^{-3}$ (Figure 4c).

An example of when the orbit remains within one of the domains $D_{-}$ or $D_{+}$, is presented in Figure 4d, where the numerically periodic orbit, obtained for $q=0.6$, $p=0.9$, $m=0.92$, $\epsilon ={10}^{-3}$ and $x_0\in D_{+}$, remains in $D_{+}$, while a numerically periodic orbit that visits both $D_{-}$ and $D_{+}$, obtained for $\epsilon =65\times {10}^{-4}$ and $q=0.71$, is presented in Figure 4e. While the orbit in Figure 4d is not related to discontinuity $x=0$ in the case presented in Figure 4e, although the orbit does not seem to depend on the discontinuity, the transient somehow meets the neighborhood of the discontinuity (red point).

Intensive numerical tests show that the smallest neighborhood size where the orbits could be identified is of order $\epsilon ={10}^{-3}$.

5. Conclusions

In this paper, it is shown that the fractional PWC systems (7) might have no solutions. Even if the use of the numerical integral (5) could offer a numerical solution, this could be incorrect. This characteristic is also explained computationally. A possible solution is to use the Cellina theorem, which allows the restarting of the PWC problem as a continuous one, where integral (5) can be applied.