How to empower Grünwald–Letnikov fractional difference equations with available initial condition?

. In this paper, the initial condition independence property of Grünwald–Letnikov fractional difference is revealed for the ﬁrst time. For example, the solution x ( k ) of equation G a ∇ αk x ( k ) = f ( x ( k )) , k (cid:62) a + 1 , cannot be calculated with initial condition x ( a ) . First, the initial condition independence property is carefully investigated in both time domain and frequency domain. Afterwards, some possible schemes are formulated to make the considered system connect to initial condition. Armed with this information, the concerned property is examined on three modiﬁed Grünwald–Letnikov deﬁnitions. Finally, results from illustrative examples demonstrate that the developed schemes are sharp.


Introduction
Fractional calculus, born at the end of seventeenth century, has provided many hot topics of research in many fields of science and engineering [4,19]. The reason of the success of the fractional calculus is its nonlocality. One of the possibility is the consideration of discrete time systems, which have the advantages of simple calculation and small singularity [24,25]. For the details of the most recent advances, one can refer to some excellent monographs [5,7,13] and the references cited therein as well.
The Grünwald-Letnikov definition was originally proposed for the fractional derivative (time-domain case). To avoid the inconvenience of using its limit form, the infinitesimal sampling period h was usually set as a finite small value in numerical computation process [10,12,16]. Besides, with the discretization approximation of Grünwald-Letnikov derivative, fractional-order systems can be implemented by FPGA (Field Programmable Gate Array) [15,18]. Along this way, the interpretation of Grünwald-Letnikov derivative was studied [6,17]. The initial value problem was evaluated by [1,3]. Actually, when h is finite, especially when h = 1, it does denote the fractional difference [13]. For the discrete case, [14] investigated the system order identification, [9] addressed the system pseudo state estimation, [11,20] discussed the stability analysis issue, and [8] considered the perfect control problem. Due to their efforts, a series of wonderful and meaningful results have been achieved [2].
Despite the wide use of Grünwald-Letnikov fractional difference, it has a major drawback (initial condition independence), which is the restriction on the dynamic behavior or even the stability. For example, by using the definition given in [20], the state response of G a ∇ α k x(k) = λx(k) + bu(k) is meaningless, which brings many difficulties for analysis and synthesis. To the best of the authors' knowledge, no results on this property have been reported, while it is indeed an essential property to be explored before Grünwald-Letnikov fractional difference can be used as a viable tool. Bearing this in mind, this paper plans to investigate the initial condition independence problem. It should be pointed out that this work will deepen our understanding on Grünwald-Letnikov fractional difference and enlarge its applicability.
The outline of the rest paper is organized as follows. Section 2 introduces some basic knowledge of nabla discrete fractional calculus. Section 3 reveals the initial conditions independent property of Grünwald-Letnikov fractional difference and provides three schemes to empower Grünwald-Letnikov fractional difference equations with available initial condition. Section 4 discusses several relevant definitions briefly and presents some detailed numerical examples to confirm the correctness and effectiveness of the obtained results. Finally, the paper ends with a conclusion in Section 5.
It can be found that (1) and (2) have a similar form. The distinction lies in that the nabla difference has local memory, and the nabla sum has nonlocal memory. However, to avoid the singularity, i.e., Γ(−n) = ∞, n ∈ N, the generalized binomial coefficient can be calculated as alternative expression For the special case of n = 0, let us make an assumption Extending the order n in Definition 1 from the positive integers to the positive real number, the following definition can be derived.
When the order . When the order α = −n ∈ Z − , G a ∇ −n k x(k) in (4) reduces to the classical nabla sum a ∇ −n k x(k) in (2). When the order α ∈ R − \ Z − , G a ∇ α k x(k) becomes the Grünwald-Letnikov fractional sum, which is also known as the nabla Riemann-Liouville fractional sum [7]. When the order α = n ∈ Z + , G a ∇ n k x(k) is not identical to ∇ n x(k) in (1), especially for those k − a − 1 n, since the upper limit of summation in formula (4) is a time-varying variable k − a − 1 instead of the constant n. If the number of x(k − i) in the summation is called as the memory length, it is worth noting that G a ∇ n k x(k) has a time-varying memory length, while ∇ n x(k) has a fixed memory length, e.g., , k a+2. When the order α ∈ R + \ Z + , G a ∇ α k x(k) becomes the Grünwald-Letnikov fractional difference. Due to the perfect match in format of fractional difference and fractional sum, the following relationship follows [23]: where α ∈ (n − 1, n), n ∈ Z + , k ∈ N a+1 and a ∈ R. By combining the rising function p q := Γ(p+q)/Γ(p) with the convolution operation, the Grünwald-Letnikov fractional difference/sum can be rewritten as where α ∈ R \ N, k ∈ N a+1 and u(k) * v(k) := k−a−1 i=0 Given a function x : N a+1 → R, a ∈ R, if there exists some r > 0 such that X(s) = N a {x(k)} := +∞ k=1 (1 − s) k−1 x(k + a) converges for |s − 1| < r, then the following equality holds [22]: where |s−1| < min{r, 1} and α ∈ R\Z + . Notably, the sufficient and necessary condition for lim k→+∞ x(k) = 0 is that the multiplicity of pole s = 0 for X(s) is less than 1 and all other principal poles of X(s) satisfy |s − 1| > 1.
It can be observed from (5) and (7) that This property will make it impossible to get x(k) from G a ∇ α k x(k) = λx(k), which means that the state or output feedback control for the system G a ∇ α k x(k) = f (x(k), u(k)), y(k) = g(x(k), u(k)) is meaningless. For this reason, the objective of this work is to analyze this default property and determine a recipe for eliminating this independence.

Main results
This section includes two parts, i.e., showing the independence of initial conditions and solving the initial condition independence problem.

The independent case
Consider the following system with Grünwald-Letnikov definition: where α ∈ (0, 1), k ∈ N a+1 , x(k) ∈ R n , is the pseudo state, u(k) ∈ R p is the input, y(k) ∈ R q is the output, f (·, ·) and g(·, ·) are linear or nonlinear functions.
Defining v(k) = f (x(k), u(k)) and using the property in (5), the first equation of (8) can be rewritten as , which implies that when u(k) = 0, k ∈ N a+1 , is considered, x(k), k ∈ N a+1 , cannot be calculated with initial condition x(a). Inspired by [21], a discrete time frequency distributed model is introduced here: : where z(ω, a) = 0, ω ∈ [0, +∞). It can be checked that the relation from u(k) to x(k) is equivalent to that of (8). Along this way, different equivalent representation of 1/s α could give different equivalent models of (8). For simplicity, it will not be expanded here. To get the equivalent relationship, zero initial condition should be configured for the equivalent model like (9), which is different from the Riemann-Liouville case and the Caputo case [21].
In system (8), the fractional difference is not calculated directly. To continue, consider the fractional difference operation directly where α ∈ (0, 1), v(k) is known, x(k) is to be calculated, and k ∈ N a+1 . Since s α = +∞ 0 can be calculated by the following systems: Note that x(k) can be accurately calculated without the dependance on initial condition v(a).

Remark 1.
It can be observed that no matter the system output of (8), and the difference output of (10) can be accurately calculated without considering the initial condition. Even if the order is extended from 0 < α < 1 to n − 1 < α < n, the independence of initial conditions is always the same. This point is actually implied in Definition 2. More specially, the calculation of G a ∇ α k x(k), k ∈ N a+1 , only needs the value of x(k), k ∈ N a+1 .
Up till now, the initial conditions independence property of system (8) has been discovered. To further show its influence, let us consider the following linear system governed by with α > 0, k ∈ N a+1 , and b = 0. x(k) can be calculated as where It is worth mentioning that in both (12) and (13), only the input response can be found, and there is no state response. This result is quite coincident with the previous discussion.
Remark 2. One point worth emphasizing is that the exact calculation of x(k) in (8). Two methods will be discussed.
The first one is to solve fractional difference equation directly. By using (4) and (6), it follows where α, a, λ, b, u(k), k ∈ N a+1 , are known, and x(k), k ∈ N a+1 , is to be calculated. (14) can be rewritten as From this the value of x(k) can be calculated sequentially.
, which involves the gamma function. When calculating the value of gamma function using MATLAB code gamma(), one has Γ(171) = 7.257415615308056 · 10 306 and Γ(171) = ∞. To avoid the overflow problem, two available ways can be taken. The one is recursive computation by using the property The other is equivalent transform by using the relation where the code gammaln() is helpful. Some other treatment can be made before introducing logarithmic function when Γ(i − α) < 0.
The second one is to solve fractional sum equation indirectly. By using (5), (11) can be rewritten as , while the unexpected case Γ(i − α) < 0 will not appear any more.

The dependent case
The Scheme 1 Consider the following system: where α > 0, k ∈ N a+1 , λ ∈ R and m ∈ Z + . Using the calculating approach in Remark 2 yields It can be observed that when calculating x(k) for k ∈ N a+1 , the initial conditions x(a + 1 − m), x(a + 2 − m), . . . , x(a) are all adopted. This scheme just hints how to solve the initial conditions independence issue when applying the Grünwald-Letnikov fractional difference.
Defining X(s) := N α {x(k)}, the time-delay property of the nabla Laplace transform follows [22] from which the state response in frequency domain becomes Note that (18) coincides with (16) and (17). System (15) will be asymptotically stable if λ is designed such that the multiplicity of root s = 0 is less than 1 or all other the principal roots lie in the region |s − 1| > 1.
Scheme 2 Consider the following system: with α > 0, k ∈ N a+1 , λ ∈ R, and m ∈ Z + . Using the calculating approach in Remark 2 yields As such, the calculation of x(k) for k ∈ N a+m+1 strongly depends on the initial conditions x(a + 1), x(a + 2), . . . , x(a + m). By applying the following time advance property [22] N X(s) := N α {x(k)} can be expressed as To sum up, in face of the independence of initial conditions for Grünwald-Letnikov fractional difference, the autonomous system G a ∇ α k x(k) = λx(k) cannot work, while system (15) or (19) can be used as a substitute. The stable condition of the system (19) is the same with the system (15).

Remark 3.
Combing the two methods, a strange but available system follows G a ∇ α k x(k) = λx(k + m), where α > 0, k ∈ N a+1 , λ ∈ R, λ = 0, and m ∈ Z + . In this regard, one has the relationships x(k + m) = (1/λ) x(a + i), which means that the initial conditions x(a + 1), x(a + 2), . . . , x(a + m) are needed to calculate x(k) for k ∈ N a+m+1 . Different from the methods in (15) and (19), this method is not directly applicable to the nonlinear case, such as G a ∇ α k x(k) = f (x(k + m)).

Scheme 3
Afterwards, let us continue to consider the system G a ∇ α k x(k) = λx(k). If we assume k ∈ N a+1+m , m ∈ Z + and take x(a + 1), x(a + 2), . . . , x(a + m) as the initial condition, the system becomes with α > 0, k ∈ N a+1 , λ ∈ R, and m ∈ Z + . From this the initial condition has been exerted successfully.
Using the calculating approach in Remark 2 yields
It can be found that the initial conditions x(a+1), x(a+2), . . . , x(a+m) are implied in X(s). The corresponding stability condition becomes that the root of s α − λ = 0 lies in the region |s − 1| > 1. Compared with systems (15) and (19), it is more convenient to use frequency tool in system (21). However, the inconvenience is that this method is not applicable to the nonlinear case.
Till now, four available schemes have been presented already. Note that the systems in (15), (19), (21), and (22) only need the initial condition about x(k) not the difference or sum of x(k), which is helpful in practical applications. The nature of them is shifting. Along this way, there are various of methods to introduce the initial condition. The basic principle might be to connect the future signal with the past signal and calculate the unknown value with the known value.
Remark 4. The Grünwald-Letnikov fractional difference/sum adopted in [8,9,11,13,14] is which is a special case of (22) with m = 1. In that case, the initial conditions independence can also be solved. When k = a, system G a−1 ∇ α k x(k) = λx(k), k ∈ N a , becomes x(a) = λx(a), which is unexpected. To avoid this confusing phenomenon, the discrete time k ∈ N a+1 should be applied.

Further discussion
In this section, the developed results will be discussed for some variants of discrete Grünwald-Letnikov fractional calculus. Afterwards, three detailed numerical examples are provided to verify the correctness of the proposed theory.

Extension
In the previous discussion, the basic is assumed. Hereafter, three other definitions will be considered.
By using the nabla Laplace transform, it follows Although a new parameter λ is introduced, there is no item regarding to the initial condition in (23).
The elegant formula N a { G a ∇ α(k) k x(k)} = s α(k) X(s) does not hold, while this definition can be rewritten as G a ∇ α(k) k From this, when G a ∇ α(k) k x(k), k ∈ N a+1 , is calculated, no more information of x(k), k < a + 1 is needed. To be more precise, all the fractional differences in Definitions 2, 3, and 4 are independent on the initial condition. Fortunately, by using the schemes in Section 3.2, the initial condition dependence can be attained.
Notably, the equality G k−N ∇ N k x(k) = ∇ N x(k) strongly implies the possible connection to the initial condition. Another equality G k−N ∇ α k x(k) = G a ∇ α k x(k) − G a ∇ α k−N x(k) also shows that the calculation of G k−N ∇ α k x(k) for k ∈ N a+1 depends on the history information x(k), k < a + 1. Different from the previous definitions, here the memory length is a constant value N + 1 instead of a time varying value k − a. Along this way, both G a ∇ α k+N x(k) and G k ∇ α k+N x(k) might be the alternative schemes to generate the initial conditions.

Examples
Example 1. Consider the following five systems: It can be observed that with the shifting principle, the initial conditions are available to calculate x(k), k ∈ N 4 , in cases 1-5. Moreover, for different σ, x(k) can converge to 0 as k increases. In this connection, the results in Section 3.2 have been validated.

Example 2.
To evaluate the practicability of the proposed result, the following gradient methods are considered: .

(a)-3(f).
These figures clearly show that all the presented cases could find the exact extreme point x m for α < 2. Since the stability domain of (15) or (19) is difficult to obtain and therefore cases 2 and 4 are not included in the case of α > 2. For cases 1-6, the overshoot appears when α > 1. Notably, the initial value x(4) in case 1 is different from others (see yellow circle and red circle) since the calculation is x 1 (4) = 2(ρη/(2ρη + 1))x m , not the setting value 1.
In what follows, the quartic function J(x) = 3x 4 − 4x 3 − 12x 2 shown in Fig. 4 is considered. Because the gradient dJ(x)/dx = 12x(x − 2)(x + 1) is nonlinear, the exact solutions of cases 1, 3, 5, and 6 cannot be calculated directly. Therefore, only cases 2 and 4 are discussed hereinafter. With the following parameters α = 1.2, ρ = 0.016, and σ = −2.5, −2.4, . . . , 2.5, the simulation results can be obtained as Figs. 5(a)-5(b).  Since dJ(x)/dx = 0, the exact extreme points are x = −1, x = 0, and x = 2. More specially, x = −1 is a local minima point, x = 0 is maximum point, and x = 2 is a global minima point. From Fig. 5(a) it can be found that all the curves converge to the exact extreme points. Likewise, cases 2 and 4 give the similar results. The relationship between the convergent point and the initial condition is clearly shown in Fig. 5(b), which means that the initial conditions independence disappears with the special design.
Example 3. Construct the following fractional-order neural network:  With the given conditions, the phase diagram in the three dimensional space is depicted in Fig. 6(a), and the trajectories of x(k) over time is shown in Fig. 6(b). It could see that the proposed initial condition scheme is effective.
Before ending this section, the main contributions of this paper will be summarized to make it more readability.
(i) The independence on the initial condition is shown for Grünwald-Letnikov fractional difference equation. (ii) The equivalent frequency distributed models for a fractional-order system with Grünwald-Letnikov definition are derived. (iii) The frequency distributed model is also developed for Grünwald-Letnikov fractional difference calculation. (iv) The unavailability of state response for a linear fractional-order system is analyzed. (v) The dependent case is considered, which connects the system response with the initial condition. (vi) The independence problem is also discussed for some variants of Grünwald-Letnikov definition.

Conclusions
In this paper, the initial conditions independent property has been studied for the first time. Inspired by the meaningless system G a ∇ α k x(k) = λx(k), k ∈ N a+1 , the initial value problem on Grünwald-Letnikov fractional difference is discussed from the definition, the distributed frequency model, and the time-domain response. To improve the practicability, four schemes are developed by connecting the initial conditions to the considered system. Apart from this, three illustrative examples are provided to confirm our findings.
It is believed that the proposed principles could greatly enrich the comprehension of Grünwald-Letnikov fractional difference and facilitate its applications. Future research efforts will be directed towards the following topics.
(i) Discuss the stability condition of systems G a ∇ α k x(k) = λx(k − m), G a ∇ α k x(k) = λx(k + m), G a ∇ α k x(k + m) = λx(k) further and analyze the dynamic behaviour of its time-domain system response.
(ii) Extend the related results to other Grünwald-Letnikov based definitions, i.e., the right-hands, the nonsingular kernel case, time-scales case, etc. (iii) Apply the initial conditions dependent Grünwald-Letnikov difference to practical applications, such as modeling, controlling, filtering, optimization, etc.