# A Plea for the Integration of Fractional Differential Systems: The Initial Value Problem

Nezha Maamri, Jean-Claude Trigeassou

► To cite this version:

Nezha Maamri, Jean-Claude Trigeassou. A Plea for the Integration of Fractional Differential Systems: The Initial Value Problem. Fractal and Fractional, 2022, 6 (10), pp.550. 10.3390/fractalfract6100550. hal-04028402

HAL Id: hal-04028402

https://hal.science/hal-04028402

Submitted on 14 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Article

# A Plea for the Integration of Fractional Differential Systems:
The Initial Value Problem

Nezha Maamri <fnref n="1" />,* and Jean-Claude Trigeassou <fnref n="2" />

**Abstract:** The usual approach to the integration of fractional order initial value problems is based on the Caputo derivative, whose initial conditions are used to formulate the classical integral equation. Thanks to an elementary counter example, we demonstrate that this technique leads to wrong free-response transients. The solution of this fundamental problem is to use the frequency-distributed model of the fractional integrator and its distributed initial conditions. Using this model, we solve the previous counter example and propose a methodology which is the generalization of the integer order approach. Finally, this technique is applied to the modeling of Fractional Differential Systems (FDS) and the formulation of their transients in the linear case. Two expressions are derived, one using the Mittag-Leffler function and a new one based on the definition of a distributed exponential function.

**Keywords:** initial value problem; fractional differential systems; fractional integral equation; infinite-state approach; Riemann-Liouville integral; frequency distributed exponential function

<figure><img src="image_3.png" /></figure>

check for
updates

**Citation:** Maamri, N.; Trigeassou, J.-C. A Plea for the Integration of Fractional Differential Systems: The Initial Value Problem. *Fractal Fract.* **2022**, *6*, 550. https://doi.org/10.3390/fractalfract6100550

Academic Editors: Thach Ngoc Dinh, Shyam Kamal, Rajesh Kumar Pandey and Norbert Herencsar

Received: 14 July 2022

Accepted: 19 September 2022

Published: 28 September 2022

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

<figure><img src="image_4.png" /></figure>

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

## 1. Introduction

The integration of Fractional Differential Equations (FDE) and Systems (FDS) is considered to be a well-founded and approved topic for most fractional calculus researchers. Therefore, the title of the paper appears as an ingenuous and unrealistic objective to revisit an established mathematical result. Nevertheless, our purpose is to provide an objective analysis of this fundamental problem and to formulate a satisfactory solution to fractional-order initial-value problems.

In fact, the initial-value problem, or Cauchy problem, is obviously trivial in the integer order case [1,2]. On the other hand, the solution of the fractional-order case appears as a generalization of the integer-order one. However, due to the multiplicity of fractional-order derivative definitions, researchers have considered it necessary to adapt the classical approach by referring to a particular derivative and its corresponding initial conditions [3–5]. Practically, most of the time, the Caputo derivative [3,6] is used because its “initial conditions” can be physically interpreted. Many critics have already addressed this choice, based on initialization considerations [7–15]. In those papers, the authors emphasize the inability of the Caputo derivative technique to solve the initialization problem, but, contrary to the history function technique [9,10,16–21] and the infinite-state approach [13,22–24], they do not provide a solution to this problem. Recently, some solutions based on new fractional derivatives (see, for example, [25–28]), which are in fact local derivatives [29,30], have been proposed. Practically, the direct consequence of these multiple choices is that different theoretical free responses are possible for the same FDE/FDS problem, which is of course physically inconsistent.

Our objective in this paper is to prove (in fact to recall) theoretically, using an elementary initial value problem, that the solution predicted by the Caputo derivative approach

<footnote n="1">1 LIAS Laboratory, Poitiers University, 2 rue Pierre Brousse, CEDEX 9, 86073 Poitiers, France</footnote>

<footnote n="2">2 IMS Laboratory, Bordeaux University, 351 Cours de la libération, CEDEX, 33403 Talence, France</footnote>

<footnote>* Correspondence: nezha.maamri@univ-poitiers.fr</footnote>
leads to a false free response. Then, we treat the same example with the frequency-distributed model of the fractional integrator [31–34]. We demonstrate that using a distributed initial condition, in fact that of the fractional integrator, provides the good solution to the considered problem. The conclusion of this analysis is that any fractional-order initial-value problem has to be treated such as in the integer-order case, using essentially the fractional integrator and its distributed initial state or its initialization function. Then, this technique is applied to the modeling of FDE/FDS and the formulation of their transients in the linear case. Two expressions are derived, one using the classic Mittag–Leffler function and a new one based on the definition of a distributed exponential function.

The theory developed in the paper is not a new one, since the first paper [35] related to the fractional integrator was published in 1999. Since that original publication, this research has been applied to the modeling and identification of real-world diffusive processes: electrochemical [36], thermal [37], and rotor skin effect, see chapter 5, volume 1 of [34]. The modeling of fractional systems based on the fractional integrator, known as the infinite-state approach, has been presented in several articles, see, for example, [22,33,38], with a particular focus on system initialization [23,24]. Moreover, it has been applied to the stability analysis of linear and nonlinear systems with a distributed formulation of the Lyapunov function [39]. The theory of the infinite state approach and its applications to various domains of control theory are presented in a two-volume monograph [34]. However, in spite of its contributions to initialization and Lyapunov system stability, this theory is ignored or considered as an exotic contribution to fractional calculus, although it has been adopted by researchers for initialization purposes [21,40–42] and Lyapunov stability analysis [43–48]. Moreover, although the pseudo initial conditions of the Caputo derivative are frequently criticized [8,9,16,20], mainly for their use in system initialization [7,14,15,21], they are still used because they provide apparently simple solutions. Consequently, there is an important challenge to provide a general and satisfactory solution to the initialization problem, using the same approach as in the integer-order case, where the initial conditions are those of the fractional integrator.

Thus, this paper intends to treat the FDE initial-value problem with a new and theoretical presentation of the infinite state approach, demonstrating that we do not have to refer to any fractional derivative and, on the contrary, can focus on the Riemann–Liouville integral and its distributed initial conditions. It is important to note that the authors have privileged a theoretical formalism contrary to their previous publications, where numerical simulations were abundantly used. So, the reader can refer to these previous papers to find numerical illustrations related to initialization. A restricted version of the paper has already been published in a recent conference [49].

The paper is composed of six sections and a conclusion. Section 1 is the introduction. Sections 2–4 present the materials and methods related to initial-value problems and the infinite state approach. In Section 5, an elementary counter-example permits us to invalidate the usual Caputo derivative initial value approach. In Section 6, the frequency-distributed integrator model is used to solve the previous counter-example and to formulate a new approach to the FDE initial-value problem. This methodology is used in Section 6 to express the dynamics of the general FDS initial-value problem.

## 2. Materials and Methods

### 2.1. The ODE Initial-Value Problem (or Cauchy Problem [1])

Let us consider the following Ordinary Differential System:

$$\frac{dx(t)}{dt} = f(x(t), u(t)) \qquad (1)$$

where $x(t) = x(0)$ at $t = 0$ is the initial value.
The Picard-Lindelöf theorem [50] guarantees the existence of a unique solution to (1). The principle of this theorem consists in reformulating the problem as an equivalent integral equation:

$$x(t) = I_t^1 \left[ \frac{dx(t)}{dt} \right] + x(0) = \int_0^t f(x(\tau), u(\tau))d\tau + x(0) \quad (2)$$

and to construct a sequence of functions

$$\phi_{k+1}(t) = \int_0^t f(\phi_k(\tau), u(\tau))d\tau + x(0) \quad \text{with} \quad \phi_1(t) = x(0), \quad (3)$$

which converges to the solution of (1) and thus to the solution of the initial-value problem. Such a construction is called Picard's method [51] or the method of successive approximations. In the linear and multidimensional case, Equation (1) can be expressed as:

$$\frac{d\underline{x}(t)}{dt} = A\underline{x}(t) + Bu(t) \quad x(t) = x(0) \quad \text{at } t = 0, \quad (4)$$

where $x(t) \in R^N$, and $A$ and $B$ are matrices of appropriate dimensions. It is well known that its solution, based on the exponential matrix function or transition matrix

$$\Phi(t) = e^{At} \quad \text{with} \quad \Phi(t) \in R^{N \times N}$$

is given by [52]:

$$x(t) = \Phi(t)x(0) + \int_0^t \Phi(t-\tau)Bu(\tau)d\tau. \quad (5)$$

## 2.2. The FDE/FDS Initial Value Problem

Let us consider the elementary FDE

$$D_t^n(x(t)) = f(x(t), u(t)) \quad 0 < n < 1, \quad (6)$$

where $n$ is the fractional order and $x(t) = x(0)$ at $t = 0$.

Contrary to the integer-order case, several approaches are derived from the fractional derivative definitions of $D_t^n(x(t))$. The main popular ones are the Caputo and Riemann-Liouville derivatives [3].

Practically, Equation (6) is integrated with the Caputo derivative definition, since its initial condition is considered equal to $x(0)$.

Then, in order to prove the existence and the uniqueness of the solution $x(t)$ of (6), Picard's method [3-5] is frequently used.

In the linear multidimensional commensurate order case, Equation (6) becomes:

$$\left\{ \begin{array}{l} D_t^n(\underline{x}(t)) = A\underline{x}(t) + Bu(t) \\ y(t) = C\underline{x}(t) \end{array} \right. \quad 0 < n < 1, \quad (7)$$

where $x(t) \in R^N$, and $A$ and $B$ are matrices of appropriate dimensions. The general solution of (7), expressed in terms of the Mittag-Leffler matrix function [53]

$$\Phi(t) = E_{n,1}(At^n), \quad \Phi(t) \in R^{N \times N}$$

is

$$\underline{x}(t) = \Phi(t)\underline{x}(0) + \int_0^t \Phi(t-\tau)B\tilde{u}(\tau)d\tau \quad \text{with} \quad \tilde{u}(\tau) = D^{1-n}(u(\tau)). \quad (8)$$
As mentioned in the introduction, the main objective of this paper is to revisit the integration of the FDE/FDS initial-value problem, using the infinite-state approach, which is directly related to the integer order ODE case and does not need any derivative definition, as it is exhibited in Section 4.

### 3. Integration of FDE/FDS Based on Derivative Definitions

#### 3.1. Riemann–Liouville Integral

The fractional integral of a function $v(t)$, also called the Riemann–Liouville integral is defined by

$$x(t) = {}_0I_t^n(f(t)) = \int_0^t \frac{(t-\tau)^{n-1}}{\Gamma(n)} v(\tau)d\tau \quad 0 < n < 1, \quad (9)$$

where $(n)$ is the gamma function.

The fractional integral is in fact a convolution integral, characterized by the impulse response or Kernel, $h_n(t)$, such that:

$$h_n(t) = \frac{t^{n-1}}{\Gamma(n)} \quad \text{and} \quad x(t) = h_n(t) * v(t). \quad (10)$$

Using the Laplace transform, we obtain

$$L\{h_n(t)\} = \frac{1}{s^n}, \quad (11)$$

where $\frac{1}{s^n}$ corresponds to the fractional order integration operator.

#### 3.2. Fractional Derivatives Definitions

Contrary to the fractional integral, the fractional derivative is not uniquely defined. Usually, two main derivatives are considered; since they are used for the integration of FDE/FDS, we focus on the case $0 < n < 1$ [34].

##### 3.2.1. Caputo Derivative Definition

This definition corresponds to first differentiate $x(t)$ and then calculates a fractional integral with order $(1-n)$. Since $0 < n < 1$, then $0 < 1-n < 1$.

$${}^cD_t^n(x(t)) = I_t^{1-n}\left(\frac{dx}{dt}\right) = h_{1-n}(t) * \frac{dx}{dt} \quad (12)$$

Definition (12) clearly shows that the Caputo derivative corresponds to the Riemann–Liouville integral of the derivative of $x(t)$.

In the Laplace domain, the Caputo derivative definition leads to

$$L\{({}^cD_t^n(x))\} = \frac{1}{s^{1-n}}L\left\{\frac{dx(t)}{dt}\right\} = \frac{1}{s^{1-n}}[sX(s) - x(0)]. \quad (13)$$

$$\text{So, } X(s) = \frac{1}{s^n} L\{ {}^cD_t^n(x) \} + \frac{x(0)}{s}, \quad (14)$$

Or in the time domain,

$$x(t) = I_t^n({}^cD_t^n(x)) + x(0). \quad (15)$$

Thus, the solution of FDE/FDS (6) according to the Caputo approach is

$$x(t) = {}_0I_t^n(f(x(t), u(t))) + x(0) \quad \text{for } t \ge 0 \ 0 < n < 1, \quad (16)$$
where $x(0)$ is interpreted as the initial condition of the FDE/FDS and also as the initial value of the Riemann–Liouville integral.

This simple integral Equation (16), apparently equivalent to the integer order case (2), has made the success of the Caputo derivative approach.

### 3.2.2. Riemann–Liouville Derivative Definition

This definition shows that the Riemann–Liouville derivative corresponds to the integer-order derivative of the Riemann–Liouville integral of $x(t)$.

$$^{RL}D_t^n(x(t)) = \frac{d}{dt}[I_t^{1-n}(x)] = \frac{d}{dt}[h_{1-n}(t) * x(t)] \quad (17)$$

Using the Laplace transform, we obtain

$$L\{^{RL}D_t^n(x)\} = s(L\{I_t^{1-n}(x)\}) - g(0) = s(\frac{1}{s^{1-n}}X(s)) - g(0) = s^nX(s) - g(0)$$
with $g(0) = \{\lim_{s \to 0} I_t^{1-n}(s)\}_{t=0}$.

Since $g(0)$ does not have a physical and direct interpretation, the Riemann–Liouville derivative is generally not used to integrate system (6) or (7).

It is important to note that the two derivative definitions only require integer-order differentiation ($\frac{d}{dt}$) and fractional-order integration ($I^{1-n}$). So, contrary to a common belief, the basic operation of fractional calculus is not fractional differentiation but fractional integration.

### 3.2.3. The Grünwald–Letnikov Derivative

Instead of the two previous fractional derivatives, it is possible to use the Grünwald–Letnikov (G.L.) derivative, with appropriate initial conditions. In fact, it is preferable to consider the G.L. integrator that corresponds to the discretization of the Riemann–Liouville integral.

The Nth integer order Euler derivative of $x(t)$ is defined as

$$(D^N(x(t)))_{t=kT_e} = \lim_{T_e \to 0} \frac{(1-q^{-1})^N}{T_e^N} x_k, \quad (18)$$

where $T_e$ is the sample time, $x_k = x(kT_e)$, and $q^{-1}$ is the delay operator.

The generalization to the fractional order case provides the Grünwald–Letnikov derivative

$$(^{GL}D^n(x(t)))_{t=kT_e} = \lim_{T_e \to 0} \frac{(1-q^{-1})^n}{T_e^n} x_k, \quad 0 < n < 1. \quad (19)$$

Since $L\{q^{-1}\} = e^{-T_e s}$, we obtain

$$L\{^{GL}D^n(x(t))\} = \lim_{T_e \to 0} \frac{(1 - e^{-T_e s})^n}{T_e^n} L\{x(t)\} = s^n X(s). \quad (20)$$

Notice that

$$\frac{(1 - q^{-1})^n}{T_e^n} = \frac{1}{T_e^n} \left[ 1 + \sum_{i=0}^{\infty} \alpha_{i, GL} q^{-i} \right], \quad (21)$$

with $\alpha_{i, GL} = (-1)^{i} \frac{n}{2} \frac{n-1}{3} \frac{n-2}{3} \dots \frac{n-(i+1)}{i}$,

$$\text{so } (^{GL}D^n(x(t)))_{t=kT_e} = \frac{1 + \sum_{i=0}^{\infty} \alpha_{i, GL} q^{-i}}{T_e^n} x_k, \quad (22)$$

which is the Moving Average formulation of the Grünwald–Letnikov derivative.
Reciprocally, we can define the Grünwald–Letnikov integral operator [34] as

$$^{GL}I^n(f_k) = \frac{T_e^n q^{-1}}{1 + \sum_{i=0}^{\infty} \alpha_{i, GL} q^{-i}} f_k, \qquad (23)$$

which is the Auto-Regressive formulation of the Grünwald–Letnikov integrator.

Notice that

$$L\{^{GL}I^n(f(t))\} = \lim_{T_e \to 0} \frac{T_e^n e^{-T_e s}}{(1 - e^{-T_e s})^n} L\{f(t)\} = \frac{1}{s^n} F(s), \qquad (24)$$

which means that the Grünwald–Letnikov integrator is the time discretization of the Riemann–Liouville integral.

Consider now the elementary FDE initial value problem (6):

$$D_t^n(x(t)) = f(x(t), u(t)) \quad 0 < n < 1.$$

Using the Grünwald–Letnikov integrator, we can express {$x(k)$} as:

$$\{x(k+1)\} = {}^{GL}I^n(f(\{xk\}, \{uk\})) + g\{x_{init}\}, \qquad (25)$$

where {$x_{init}$} = {$x(0), x(-1), \dots, x(-i), \dots, x(-\infty)$}.

This means that the initial conditions are composed of all the past values of $x(-i)$, since $k = -\infty$.

Practically, this technique is used for the numerical simulation of the FDE/FDS problem.

The interested reader can refer to chapter 3 volume 1 of [34], where different initializations of the G.L. integral and the short memory principle [3] are analyzed.

# 4. The Infinite State Approach

## 4.1. Introduction

The infinite-state approach (do not confuse it with diffusive representation, see chapter 7 of [34]) is a modeling technique based on the fractional-integration operator, which is at the heart of any modeling and simulation system, either integer or fractional order, linear or non-linear [22,52].

## 4.2. The Frequency-Distributed Model of the Fractional Integrator

As shown in Section 3, the Riemann–Liouville integral of a function $v(t)$ is defined as the convolution of $v(t)$ with the impulse response $h_n(t)$ of the fractional integrator.

Another expression of $h_n(t)$ can be derived from the inverse Laplace transform of $\frac{1}{s^n}$ [31,32,34], for $0 < n < 1$, i.e., $h_n(t) = L^{-1}\{\frac{1}{s^n}\}$.

Using a Bromwich contour, we can write (see [34] and the references therein):

$$h_n(t) = \begin{cases} \frac{\gamma + j\omega}{2\pi j} \int_{\gamma - j\omega}^{\gamma + j\omega} \frac{1}{s^n} e^{st} ds & \text{for } t > 0 \\ 0 & \text{for } t < 0 \end{cases}.$$

Thus, we obtain

$$h_n(t) = \int_0^\infty \mu_n(\omega) e^{-\omega t} d\omega = \frac{t^{n-1}}{\Gamma(n)} \quad \text{for } 0 < n < 1 \qquad (26)$$
with $\mu_n(\omega) = \frac{\sin(n\pi)}{\pi} \omega^{-n}$
Note that, in the particular case where $v(t) = (t)$ is an impulse function, the output $x(t)$ corresponds to the impulse response $h_n(t)$ and is provided by the following distributed integer-order differential system

$$ \left\{ \begin{array}{l} \frac{\partial z(\omega, t)}{\partial t} = -\omega z(\omega, t) + \delta(t) \quad \omega \in [0, +\infty) \\ h_n(t) = \int_0^\infty \mu(\omega) z(\omega, t) d\omega \end{array} \right. \qquad (27) $$

Thus $z(\omega, t) = e^{-\omega t}$, which leads to $h_n(t) = \int_0^\infty \mu(\omega) e^{-\omega t} d\omega$.

More generally, for any input $v(t)$, the corresponding output $x(t)$ of the fractional integrator is provided by the following distributed frequency system:

$$ \left\{ \begin{array}{l} \frac{\partial z(\omega, t)}{\partial t} = -\omega z(\omega, t) + v(t) \quad \omega \in [0, \infty) \\ x(t) = \int_0^\infty \mu_n(\omega) z(\omega, t) d\omega \\ \text{with } \mu_n(\omega) = \frac{\sin(n\pi)}{\pi} \omega^{-n} \end{array} \right. \qquad (28) $$

It is fundamental to notice that the original model of the fractional integrator has been transformed into an infinite-dimension integer-order differential system (28), where integer-order differentiation $\frac{\partial}{\partial t}$ has been substituted to fractional-order differentiation, and where the fractional order $n$ appears in the weighting function $\mu_n(\omega)$. This means that the fractional integrator $\frac{1}{s^n}$ is an infinite-dimension linear system [54]. These models are the two “faces” of the fractional integrator.

In fact, according to linear system theory [52], the fractional integrator has two types of models, as does any other linear system:

- the Equation $X(s) = \frac{1}{s^n}V(s)$ is the input/output representation of the fractional integrator, characterized by its impulse response $h_n(t)$ and its frequency response $\frac{1}{(j\omega)^n}$.

- the distributed differential system (28) is the infinite-dimension state-space model of the integrator, where the internal state $z(\omega, t)$ permits a complete representation of system dynamics and particularly its free response from an initial condition $z(\omega, 0)$.

**Remark 1:** let us consider the Laplace transform of (27):

$$ L\{h_n(t)\} = L\left\{\int_0^\infty \mu(\omega)e^{-\omega t}d\omega\right\} = \int_0^\infty \mu_n(\omega)L\{e^{-\omega t}\}d\omega = \int_0^\infty \mu_n(\omega)\frac{1}{s+\omega}d\omega $$

The previous equality and Equation (11) lead to

$$ \frac{1}{s^n} = \int_0^\infty \mu_n(\omega) \frac{1}{s+\omega} d\omega \quad 0 < n < 1. \qquad (29) $$

This relation exhibits that the fractional integrator is composed of an infinity of modes $\omega$, ranging from 0 to $+\infty$, whereas the integer-order integrator corresponds to only one mode situated at $\omega = 0$. Figure 1 displays the graphic representation of Equation (28). Note that the distributed differential Equations of (28) correspond to the first-order systems displayed in Figure 1. Due to the distributed nature of the differential system, the graph of Figure 1 is composed of an infinity of first-order systems.
<figure><img src="image_5.png" /><figcaption>Figure 1. The frequency-distributed model of the fractional integrator.</figcaption></figure>

**Remark 2:** The classical Laplace transform of the integer order derivative is known as
$L\left\{\frac{dx(t)}{dt}\right\} = sX(s) - x(0)$, which corresponds in fact to the relation $X(s) = \frac{1}{s}L\left\{\frac{dx(t)}{dt}\right\} + \frac{x(0)}{s}$

i.e., in the time domain, to $x(t) = \int_0^t \frac{dx(\tau)}{d\tau} d\tau + x(0)$, which means that $x(0)$ is not the initial condition of the derivative but is, in fact, the initial condition of the integrator which has memorizing capability.

4.3. Transients of the Fractional Integrator

Consider the Laplace transform of (28):

$$
\left\{
\begin{array}{l}
sZ(\omega, s) - z(\omega, 0) = -\omega Z(\omega, s) + V(s) \quad \omega \in [0, \infty) \\
X(s) = \int_{0}^{\infty} \mu_{n}(\omega) Z(\omega, s) d\omega
\end{array}
\right.
\qquad (30)
$$

where $z(\omega, 0)$ is the initial value of $z(\omega, t)$ at $t = 0$.

$$
\text{This means that } Z(\omega, s) = \frac{z(\omega, 0)}{s+\omega} + \frac{V(s)}{s+\omega}
$$

$$
\text{and } X(s) = \int_{0}^{\infty} \mu_{n}(\omega) \frac{z(\omega, 0)}{s+\omega} d\omega + \int_{0}^{\infty} \mu_{n}(\omega) \frac{V(s)}{s+\omega} d\omega \quad (31)
$$

$$
\text{Since } \frac{1}{s^n} = \int_{0}^{\infty} \mu_{n}(\omega) \frac{1}{s+\omega} d\omega \quad 0 < n < 1
$$

we can write in the time domain:

$$
x(t) = \int_{0}^{\infty} \mu_{n}(\omega) z(\omega, 0) e^{-\omega t} d\omega + {}_{0}I_{t}^{n}(v(t)) \quad (32)
$$
where:

- $x(t) = \int_{0}^{\infty} \mu_{n}(\omega)z(\omega, 0)e^{-\omega t}d\omega$ is the free response of the fractional integrator initialized by the distributed initial conditions $z(\omega, 0) \quad \forall \omega \in [0, \infty)$.

- $I_{t}^{n}(v(t))$ is the forced response of the fractional integrator caused by the input $v(t)$.

Previously, using the definition of the Caputo derivative “initial condition”, we wrote $(20) x(t) = I_{t}^{n}(v(t)) + x(0)$.

The conclusion is that this expression of the free response is wrong, since $\int_{0}^{\infty} \mu_{n}(\omega)z(\omega, 0)e^{-\omega t}d\omega$ is the initialization function of the integrator. In fact, $x(0) = \int_{0}^{\infty} \mu_{n}(\omega)z(\omega, 0)d\omega$, and Equation (15) is correct only at $t = 0$ and is wrong for $t > 0$.

The conclusion is that the fractional-integrator transients require to refer to its distributed model.

Basically, the initial condition of the differential system $D_{t}^{n}(x(t)) = f(x(t), u(t))$ is related to the initial condition of the fractional integrator $\frac{1}{s^{n}}$ used for the integration of the FDE/FDS, not to the pseudo-initial condition of any fractional derivative. Notice that if this FDE/FDS is related to a real system, its dynamics must not depend on the fractional derivative definition choice of the user.

## 5. A Counter Example

In previous papers related to the infinite-state representation, we have already demonstrated that the so-called initial conditions of the Caputo derivative are unable to correctly express the dynamics of FDE/FDS free responses. However, attracted by the apparent simplicity of the Caputo initial conditions, most fractional calculus researchers ignore the more complex (in fact not too complex) infinite-state approach.

Consequently, this paper intends to prove the fundamental errors of the usual Caputo derivative approach using an elementary theoretical counter example. Then, with the frequency-distributed integrator allowing the theoretical computation of the true free response, we prove the necessity to use frequency-distributed initial conditions to solve any FDE/FDS initial-condition problem.

### 5.1. Problem Formulation

Consider the simplest FDE initial value problem

$$ \begin{cases} D^n(x(t)) = u(t) & 0 < n < 1 \\ x(t) = x(0) & \text{at } t = 0 \end{cases} \tag{33} $$

Consider the special function $u(t)$, composed of two delayed Heaviside functions, $UH(t)$ and $-UH(t-T)$, with $u(t) = UH(t) - UH(t-T)$.

Consequently, see Figure 2

$$ u(t) = \begin{cases} U & \text{for } 0 \le t < T \\ 0 & \text{for } t \ge T \end{cases} \tag{34} $$

Moreover, assume that the system is at rest at $t = 0$, i.e., $x(0) = 0$.

**Remark 3:** The interest of this example is to create a realistic initial condition at $t = T$, where the free response can be calculated with two approaches, the first one from $t = 0$ with no ambiguity using usual fractional calculus theory and the second one from $t = T$, using Equation (15) at $t_0 = T$.
<figure><img src="image_6.png" /><figcaption>Figure 2. True free response and Caputo derivative initialization for n = 0.5.</figcaption></figure>

## 5.2. The Exact Solution

Since $x(0) = 0$, we obtain $x(t) = I_t^n(UH(t) - UH(t-T)) = x^+(t) + x^-(t)$.

So $x^+(t) = h_n(t) * UH(t) = \int_0^t \frac{(t-\tau)^{n-1}}{\Gamma(n)} UH(\tau)d\tau = \frac{t^n}{\Gamma(n+1)}UH(t)$ and $x^-(t) = -h_n(t) * UH(t-T) = -\frac{(t-T)^n}{\Gamma(n+1)}UH(t-T)$.

Consequently, see Figure 2 (for n = 0.5):

$$x(t) = \begin{cases} \frac{t^n U}{\Gamma(n+1)} & \text{for } 0 \le t \le T \\ \frac{U}{\Gamma(n+1)}[t^n - (t-T)^n] & \text{for } t \ge T \end{cases} \quad (35)$$

where $x(t)$ for $t \ge T$ represents the free response of (33) at $t_0 = T$.

## 5.3. Solution Derived from the Caputo Derivative Definition

This free response can also be expressed using Equation (15) at $t_0 = T$ with $x(t_0) = x(T)$, i.e., $x(t) = I_t^n(u(t)) + x(T)$ for $t \ge T$.

Thus (see Figure 2),

$$x(t) = x(T) \text{ for } t \ge T \quad (36)$$

This result is obviously in complete contradiction with Equation (35), i.e., $x(t) = \frac{U}{\Gamma(n+1)}[t^n - (t-T)^n]$ for $t \ge T$.

Notice that, for $n = 1$, we obtain:

$$x(T) = UT \text{ for } t = T \text{ and } x(t) = U[t - (t-T)] = UT \text{ for } t \ge T.$$

Thus, we verify that $x(t) = x(T)$ for $t \ge T$, i.e., Equation (36) is only correct in the integer-order case.

With this very simple example, we have demonstrated that Equation (15) is wrong in the fractional-order case. So, what is the reason of this basic error?
5.4. Solution Derived from the Distributed Frequency Model of the Fractional Integrator

Consider again the elementary example (33):

$$D^n(x(t)) = u(t) \quad 0 < n < 1$$

This system is supposed at rest at $t = 0$, i.e., $z(\omega, 0) \quad \forall \omega \in [0, \infty)$. So, using (30) we obtain

$$Z^{+}(\omega, s) = \frac{U}{s(s + \omega)} \quad \text{for } u(t) = UH(t)$$

$$\text{Thus, according to (29), } X(s) = \int_{0}^{\infty} \mu_{n}(\omega) \frac{1}{(s+\omega)^{\frac{1}{s}}} \frac{U}{s} d\omega = \frac{U}{s^{n+1}}$$

$$\text{and} \quad \begin{cases} x^{+}(t) = \frac{t^n U}{\Gamma(n+1)} H(t) \\ x^{-}(t) = -\frac{U(t-T)^n}{\Gamma(n+1)} H(t-T) \end{cases}$$

Obviously, we recover the same result as (35) using the distributed model. Moreover, this model allows us to express $z(\omega, t)$, i.e.,

$$z^{+}(\omega, t) = \frac{U}{\omega}(1 - e^{-\omega t})H(t)$$

$$\text{So} \quad \begin{cases} x^{+}(t) = UH(t) \int_{0}^{\infty} \frac{\mu_{n}(\omega)}{\omega} (1 - e^{-\omega t}) d\omega \\ x^{-}(t) = -UH(t-T) \int_{0}^{\infty} \frac{\mu_{n}(\omega)}{\omega} (1 - e^{-\omega(t-T)}) d\omega \end{cases} \quad (37)$$

Consequently, the free response is expressed as:

$$x(t) = U \int_{0}^{\infty} \frac{\mu(\omega)}{\omega} (e^{-\omega(t-T)} - e^{-\omega t}) d\omega \quad \text{for } t \ge T, \quad (38)$$

which is the distributed equivalent of Equation (35).

Moreover, we can verify that it is now possible to calculate the response of the integrator for $t \ge T$ using the expression (32).

So, consider the response initialized at $t = T$ with $u(t) = 0$ for $t \ge T$.

Since $z(\omega, T)$ is the initial condition at $t = T$, with

$$z(\omega, T) = z^{+}(\omega, T) = \frac{U}{\omega}(1 - e^{-\omega T}), \quad (39)$$

and since $I_t^n(0) = 0$, we obtain

$$x(t) = U \int_{0}^{\infty} \mu_{n}(\omega) z(\omega, T) e^{-\omega(t-T)} d\omega \quad \text{for } t \ge T$$

So

$$x(t) = U \int_{0}^{\infty} \frac{\mu_{n}(\omega)}{\omega} (1 - e^{-\omega T}) e^{-\omega(t-T)} d\omega \quad \text{for } t \ge T$$

And we obtain the same result as previously (38), i.e.,

$$x(t) = U \int_{0}^{\infty} \frac{\mu_{n}(\omega)}{\omega} (e^{-\omega(t-T)} - e^{-\omega t}) d\omega \quad \text{for } t \ge T$$
We can conclude that the distributed state-space model provides the exact expression of the free response using the usual tools of linear system theory. Consequently, this distributed model is the necessary tool to express transients of the fractional integrator.

Notice that numerical simulations corresponding to this counter example are available in [33].

## 5.5. Conclusions

Two main conclusions can be stated from this counter example:

- The integration of FDE/FDS based on the Caputo derivative definition (or on the Riemann-Liouville derivative) are wrong approaches leading to erroneous free responses.

- The frequency-distributed state-space model provides the exact expression of the free response using the usual tools of linear system theory. Consequently, this distributed model is the necessary tool to express transients of the fractional integrator and thus those of FDE/FDS.

Notice that the Grünwald-Letnikov approach based on relation (25) provides a correct solution to the integration of FDE/FDS. However, its initial conditions, composed of past values of $x(-i)$ since $k = -\infty$, are not easy to use, particularly for an initialization objective.

## 5.6. The Caputo Derivative Definition Revisited

We have demonstrated with the previous elementary counter example that the integration technique based on the Caputo derivative definition is unable to provide a correct expression of the free response of an elementary initial-value problem. Of course, we have pointed out the reason of this failure, i.e., $x(0)$ does not represent the initial condition of the fractional integrator. Basically, what is the origin of this error?

In order to understand why so many researchers have been misled by the so-called “initial condition” $x(0)$, we have to once again consider the definition (12) of the Caputo derivative:

$$ {}^{C}D_{t}^{n}(x(t)) = I_{t}^{1-n}\left(\frac{dx(t)}{dt}\right) $$

This derivative relies on a fractional integrator $\frac{1}{s^{1-n}}$, so we have to take into account its internal state variables $z_C(\omega, t)$ at $t = 0$ (notice that $z_C(\omega, t) \neq z(\omega, t)$) [34].

Thus, the distributed-frequency model of the Caputo derivative corresponds to that of the fractional integrator $I_t^{1-n}(\cdot)$, where, in this case, the input and the output are $\frac{dx(t)}{dt}$ and ${}^C D_t^n(x(t))$, respectively:

$$ \left\{ \begin{array}{l} \frac{\partial z_C(\omega, t)}{\partial t} = -\omega z_C(\omega, t) + \frac{dx(t)}{dt} \quad \omega[0, \infty) \\ {}^C D_t^n(x(t)) = \int_0^\infty \mu_{1-n}(\omega) z_C(\omega, t) d\omega \\ \mu_{1-n}(\omega) = \frac{\sin((1-n)\pi)}{\pi} \omega^{-(1-n)} \quad \text{and } 0 < n < 1 \end{array} \right. \quad (40) $$

with the initial condition $z_C(\omega, 0)$ and $\omega \in [0, \infty)$.

By using the Laplace transform, we can write

$$ \left\{ \begin{array}{l} Z_C(\omega, s) = \frac{z_C(\omega, 0)}{s + \omega} + \frac{L\left\{\frac{dx(t)}{dt}\right\}}{s + \omega} \quad \omega \in [0, \infty) \\ \text{with } L\left\{\frac{dx(t)}{dt}\right\} = sX(s) - x(0) \end{array} \right. $$

Thus, we obtain

$$ L\{ {}^C D_t^n(x(t)) \} = \int_0^\infty \mu_{1-n}(\omega) Z_C(\omega, s) d\omega = s^n X(s) - \frac{x(0)}{s^{1-n}} + \int_0^\infty \frac{\mu_{1-n}(\omega)}{s + \omega} z_C(\omega, 0) d\omega \quad (41) $$
We can conclude that the usual initial condition of the Caputo derivative is wrong because it does not take into account the transients of its associated integrator $\frac{1}{s^{1-n}}$, i.e., its distributed initial conditions $z_C(\omega, 0)$.

Notice that the same conclusions apply to the Riemann–Liouville derivative [34].

Consequently, the Caputo derivative approach to the integration of FDE/FDS must be rejected because it provides wrong solutions to fractional initial-value problems. This approach is wrong for two main reasons:

- The exact initial conditions of the Caputo derivative are $x(0)$ and the distributed state variable initial condition $z_C(\omega, 0)$.

- The technique based on the Caputo derivative is not natural because the true and physical initial conditions are those of the fractional integrator $\frac{1}{s^n}$ of Equation (30), i.e., $z(\omega, 0)$, such as in the integer order case.

Numerical simulations of the Caputo and Riemann–Liouville derivatives exhibiting the role of their initial conditions are available in [38] and in Volume 1 of [34].

# 6. Fractional Differential Systems Transients

## 6.1. Integration of a FDE

We demonstrated in the previous section that it is necessary to use the frequency-distributed model of the fractional integrator to take into account the transients of the free response. Thus, we have to apply the same approach for the integration of any FDE/FDS initial value problem (6,7). As we have shown, the solution is provided by the fractional integral equation, which is equivalent to the integer order case (2):

$$x(t) = I_t^n(f(x(t), u(t))) + x_0(t), \quad (42)$$

where $x_0(t)$ is the initialization function of the Riemann–Liouville integral

$$x_0(t) = \int_0^\infty \mu_n(\omega)z(\omega, 0)e^{-\omega t}d\omega. \quad (43)$$

Notice that (42) is a Volterra integral equation.

Fundamentally, Figure 3 displays the graphical representation of the integral Equation (42), which underlines the closed-loop behavior of the FDE/FDS based on the fractional integrator $\frac{1}{s^n}$ with the initialization function $x_0(t)$.

<figure><img src="image_7.png" /><figcaption>Figure 3. Closed-loop model of the FDE/FDS.</figcaption></figure>
However, we can also use the frequency-distributed model of the fractional integrator, where its input is $v(t) = D^n(x(t)) = f(x(t), u(t))$, which leads to the distributed representation of the FDE:

$$ egin{cases} 
\frac{\partial z(\omega, t)}{\partial t} = -\omega z(\omega, t) + f(x(t), u(t)) \\
x(t) = \int_{0}^{\infty} \mu_n(\omega) z(\omega, t) d\omega \\
\mu(\omega) = \frac{\sin(n\pi)}{\pi} \omega^{-n} \quad 0 < n < 1
\end{cases} 
\quad (44) $$

In this case, the solution $z(\omega, t)$ is provided by the distributed integer-order integral equation:

$$ egin{aligned} z(\omega, t) &= \int_{0}^{t} [-\omega z(\omega, t) + f(x(\tau), u(\tau))] d\tau + z(\omega, 0) \quad \forall \omega \in [0, \infty) \\ &= I_{t}^{1}(-\omega z(\omega, t) + f(x(t), u(t))) + z(\omega, 0) \end{aligned} 
\quad (45) $$

where $z(\omega, 0)$ is the initial condition of the integer order integral.

We can represent (see Figure 4) this frequency distributed system graphically, where frequency varies from $\omega = 0$ to $\omega = +\infty$, according to its Laplace transform. This graph corresponds to Figure 3, where the fractional integrator is replaced by its distributed graph of Figure 1.

<figure><img src="image_8.png" /><figcaption>Figure 4. Closed-loop model of the FDE/FDS based on the fractional integrator distributed-frequency model.</figcaption></figure>

Equations (42) and (43) and the graph of Figure 3 focus on the pseudo-state variable $x(t)$, whereas Equation (44) and the graph of Figure 4 focus on the internal state variables $z(\omega, t)$ of the integrator, which are in fact those of the FDE/FDS.

Notice that, for the isolated integrator (Figure 1), with input $v(t)$ and output $x(t)$, the state variables are decoupled and evolve independently. On the other hand, in system (44), i.e., in the graph of Figure 4, the state variables are coupled by the relation $v(t) = f(x(t), u(t))$. This means that the evolution of the state variable $z(\omega, t)$ (for the particular value) depends on all the other state variables $z(\xi, t)$ $\xi \in [0, \infty)$. Namely, the
original FDE/FDS (6) has been transformed into an infinite-dimension system of first-order differential Equations (44).

These are the two “faces” of the same problem:

- Equation (42) and Figure 3 correspond to the pseudo-state variable $x(t)$, directly taking into account the fractional order $n$.

- Equation (44) and Figure 4 correspond to the system of distributed state variables $z(\omega, t)$, whose solution is obtained through an integer-order approach. The pseudo-state variable $x(t)$ is provided by a weighted integral, where $\mu_n(\omega)$ is the link between the integer order and fractional order domains.

In Figures 3 and 4, there is no hypothesis about the nature of $f(x(t), u(t))$ which can be either linear or nonlinear. In the nonlinear case, the integral formulation of the FDE/FDS initial-value problem leads to (42) or (45). The solution of these integral equations can be obtained with Picard's method, which is currently used to solve nonlinear fractional order problems (see for instance [3]).

Of course, the nonlinear case is a very wide topic, and our objective is not to treat it in this paper. In fact, the linear case is also of fundamental interest, and we propose to revisit it, essentially to formulate free responses.

## 6.2. FDS Transients and the Mittag–Leffler Function

First, we consider the elementary FDE initial-value problem:

$$D^n(x(t)) = ax(t) \quad 0 < n < 1 \qquad (46)$$

with the distributed initial condition $z(\omega, 0) \quad \forall \omega \in [0, \infty)$.

System (46) can be transformed into the distributed frequency one:

$$\left\{ \begin{array}{l} \frac{\partial z(\omega, t)}{\partial t} = -\omega z(\omega, t) + ax(t) \\ x(t) = \int_{0}^{\infty} \mu_{n}(\omega) z(\omega, t) d\omega \end{array} \right. \qquad (47)$$

with the initial condition $z(\omega, 0)$.

Using the Laplace transform, we can write

$$sZ(\omega, s) - z(\omega, 0) = -\omega Z(\omega, s) + aX(s),$$

$$\text{i.e., } X(s) = \frac{s^n}{s^{n-a}} \int_{0}^{\infty} \mu_{n}(\omega) \frac{z(\omega, 0)}{s+\omega} d\omega$$

$$\text{or } X(s) = s \frac{s^{n-1}}{s^n - a} X_0(s), \qquad (48)$$

$$\text{with } X_0(s) = \int_{0}^{\infty} \mu_n(\omega) \frac{z(\omega, 0)}{s+\omega} d\omega$$

So,

$$x_0(t) = L^{-1} \left\{ \int_{0}^{\infty} \mu_n(\omega) \frac{z(\omega, 0)}{s+\omega} d\omega \right\} = \int_{0}^{\infty} \mu_n(\omega) z(\omega, 0) e^{-\omega t} d\omega \qquad (49)$$

Let us remember that

$$L\{E_{n,1}(at^n)\} = \frac{s^{n-1}}{s^n - a} \qquad (50)$$

where $E_{n,1}(at^n)$ is the Mittag–Leffler function:

$$E_{n,1}(at^n) = \sum_{k=0}^{\infty} \frac{(at^n)^k}{\Gamma(nk + 1)}, \qquad (51)$$
which is the generalization of the exponential function, since for $n = 1$ we obtain
$$E_{1,1}(at^1) = \sum_{k=0}^{\infty} \frac{(at^1)^k}{k!} = e^{at}.$$

Consequently, we obtain

$$x(t) = L^{-1}\{X(s)\} = \frac{d}{dt}\{E_{n,1}(at^n) * x_0(t)\}. \qquad (52)$$

We can easily generalize this result to any linear FDS initial value problem:

$$D^n(\underline{x}(t)) = A\underline{x}(t) + B\underline{u}(t) \quad 0 < n < 1, \qquad (53)$$

with the distributed initial condition $z(\omega, 0) \quad \forall \omega \in [0, \infty)$.

Defining the matrix Mittag-Leffler function:

$$E_{n,1}(At^n) = \sum_{k=0}^{\infty} \frac{(At^n)^k}{\Gamma(nk + 1)}, \qquad (54)$$

we can write

$$\underline{x}(t) = \frac{d}{dt}\{E_{n,1}(At^n) * \underline{x}_0(t)\} + \int_{0}^{t} E_{n,1}(A(t-\tau)^n) \underline{B}\tilde{u}(\tau)d\tau, \qquad (55)$$

where $\tilde{u}(\tau) = D^{1-n}(u(\tau))$ and $\underline{x}_0(t) = \int_{0}^{\infty} \mu_n(\omega)\underline{z}(\omega, 0)e^{-\omega t}d\omega$.

Obviously, the free response $\frac{d}{dt}\{E_{n,1}(At^n) * \underline{x}_0(t)\}$ of $\underline{x}(t)$ is more complex than the wrong usual one derived from (8), i.e., $\underline{x}(t) = E_{n,1}(At^n)\underline{x}(0)$.

Moreover, there is a major difficulty, i.e., the convolution between the Mittag-Leffler function and $\underline{x}_0(t)$, which is not a straightforward operation. Thus, the interest of this expression is essentially theoretical.

Notice also that the matrix Mittag-Leffler function does not verify the semi-group properties of the matrix exponential function [53,54].

**Remark 4:** The use of the Caputo derivative is based on the assumption $\underline{x}_0(t) = x(0) = cte \quad \forall t$. Is this requirement always wrong? In addition, if it is correct, what does it mean?

Consider the following system

$$D^n(x(t)) = ax(t) + bu(t) \quad 0 < n < 1$$

and suppose that it is at rest at $t = 0$ and that $u(t) = UH(t)$.

Since, in this case,

$$Z(\omega, s) = \frac{aX(s) + b\frac{U}{s}}{s + \omega} \quad \text{and} \quad X(s) = \int_{0}^{\infty} \mu_n(\omega)Z(\omega, s)d\omega, \text{ we obtain } Z(\omega, s) = b\frac{s^n}{s^n - \omega} \frac{1}{s + \omega} \frac{U}{s}.$$

Thus

$$z(\omega, \infty) = \lim_{s \to 0} Z(\omega, s) = 0 \quad \forall \omega \neq 0 \text{ and } z(0, \infty) = \lim_{s \to 0} \frac{bU}{s^{1-n}} = \infty \qquad (56)$$

On the other hand, we can also write $x(\infty) = -\frac{b}{a}U$ for a stable system, i.e., $a < 0$ in our case.

Assume that we have applied the same step input $UH(t)$, but, at $t = -\infty$, this means that at $t = 0$, we obtained $x(0) = -\frac{b}{a}U$ for a long time in the past ($t < 0$).

Consequently, the condition $x_0(t) = x(0) = cte \quad \forall t$ requires that the system has been at rest for a very long time.
Then

$$x_0(t) = \int_0^\infty \mu_n(\omega)z(\omega, 0)e^{-\omega t}d\omega = \int_0^\infty \mu_n(\omega)z(0, 0)e^{-0t}d\omega = x(0) = cte. \quad (57)$$

Since $z(\omega, \infty) = 0 \forall \omega \neq 0$.

Theoretically, the condition $x_0(t) = x(0) = cte \ \forall t$ can be achieved, but it is completely unrealistic. Moreover, it would require infinite energy [34].

**Remark 5:** Consider again the revisited definition of the Caputo derivative (see Section 5.6):

$$C_t^n D_t^n(x(t)) = I_t^{1-n} \left( \frac{dx(t)}{dt} \right) \Rightarrow \begin{cases} \frac{\partial z_C(\omega, t)}{\partial t} = -\omega z_C(\omega, t) + \frac{dx(t)}{dt}, & \omega \in [0, \infty) \\ C_t^n D_t^n(x(t)) = \int_0^\infty \mu_{1-n}(\omega) z_C(\omega, t) d\omega \\ \mu_{1-n}(\omega) = \frac{\sin((1-n)\pi)}{\pi} \omega^{-(1-n)} \quad \text{and} \quad 0 < n < 1 \end{cases}$$

with $z_C(\omega, t) = z_C(\omega, 0)$ at $t = 0$.

We have demonstrated (41) that:

$$L\{C_t^n D_t^n(x(t))\} = s^n X(s) - \frac{x(0)}{s^{1-n}} + \int_0^\infty \frac{\mu_{1-n}(\omega)}{s+\omega} z_C(\omega, 0) d\omega.$$

If $x(t) = cte$, for a long time in the past ($t < 0$), then $\frac{dx(t)}{dt} = 0$. This means that $z_C(\omega, t) = 0 \ \forall \omega$ for a long time in the past ($t < 0$), and consequently $z_C(\omega, 0) = 0 \ \forall \omega$.

Then, we can write $L\{C_t^n D_t^n(x(t))\} = s^n X(s) - \frac{x(0)}{s^{1-n}}$.

Obviously, the previous conditions are very restrictive, and the usual initial condition is wrong as soon as there is a variation of $x(t)$.

**Remark 6:** We have demonstrated (52) that:

$$x(t) = \frac{d}{dt}\{E_{n,1}(at^n) * x_0(t)\} = E_{n,1}(at^n) * \frac{d}{dt}x_0(t)$$

Since $x_0(t) = x(0) = cte \ \forall t$, we can write $x_0(t) = x(0)H(t)$, so $\frac{d}{dt}x_0(t) = x(0)\delta(t)$.

Then,

$$x(t) = E_{n,1}(at^n) * x(0)\delta(t) = E_{n,1}(at^n)x(0) \quad (58)$$

is the usual result corresponding to the Caputo derivative assumption.

The conclusion is that many results obtained with the Caputo derivative approach are not necessarily wrong, but they require the previous very restrictive assumptions related to $x(0)$.

### 6.3. FDS Transients Expressed with the Distributed Exponential Function

Previously, the problem of fractional transients has been focused on the pseudo-state variable dynamics, with no insight in the distributed state variable. So, we propose now to express the dynamics of $z(\omega, t)$. Let us start again with the elementary system (46):

$$D^n(x(t)) = ax(t) \quad 0 < n < 1$$

with the distributed initial condition $z(\omega, 0) \ \forall \omega \in [0, \infty)$.
The corresponding distributed frequency model (47) leads to:

$$ \frac{\partial z(\omega, t)}{\partial t} = -\omega z(\omega, t) + a \int_{0}^{\infty} \mu_n(\xi) z(\xi, t) d\xi, \quad \forall \omega \in [0, \infty). \quad (59) $$

Notice that we have to separate the current frequency $\omega$ from all the other ones $\xi \in [0, \infty)$. This means that the behavior of $z(\omega, t)$ depends on all the behaviors of the other state variables $z(\xi, t)$.

Let us define $\delta(\xi)$, which is the frequency Dirac impulse verifying $\int_{0}^{\infty} \delta(\xi) d\xi = 1$.

Then,

$$ \int_{0}^{\infty} \omega \delta(\xi - \omega) z(\xi, t) d\xi = \omega \int_{0}^{\infty} \delta(\xi - \omega) z(\xi, t) d\xi = \omega z(\omega, t). \quad (60) $$

Let us define

$$ \psi(\omega, \xi) = -\omega \delta(\xi - \omega) + a \mu_n(\xi) \quad (61) $$

Then, (59) can be expressed as:

$$ \frac{\partial z(\omega, t)}{\partial t} = -\omega \int_{0}^{\infty} \delta(\xi - \omega) z(\xi, t) d\xi + a \int_{0}^{\infty} \mu_{n}(\xi) z(\xi, t) d\xi = \int_{0}^{\infty} \psi(\omega, \xi) z(\xi, t) d\xi \quad (62) $$

The solution $z(\omega, t)$ requires the integration of the integer-order distributed system (59) with the initial condition $z(\omega, 0) \quad \forall \omega \in [0, \infty)$. Basically, the solution verifies the integral relation:

$$ \begin{aligned} z(\omega, t) &= \int_{0}^{t} \left[ \int_{0}^{\infty} \Psi(\omega, \xi) z(\xi, \tau) d\xi \right] d\tau + z(\omega, 0) \\ &= 0 I_{t}^{1} \left[ \int_{0}^{\infty} \Psi(\omega, \xi) z(\xi, \tau) d\xi \right] + z(\omega, 0) \quad \omega \in [0, \infty) \end{aligned} \quad (63) $$

This integration is performed with Picard's method, which is an iterative technique (3). At the first iteration, $z(\omega, t)$ is approximated by $z(\omega, 0)$.

Since $z(\omega, 0)$ and $\int_{0}^{\infty} \psi(\omega, \xi) d\xi$ are constants for $0I_{t}^{1}$, we obtain

$$ z_{1}(\omega, t) = z(\omega, 0) + I_{t}^{1} \left[ \int_{0}^{\infty} \psi(\omega, \xi) z(\omega, 0) d\xi \right] = z(\omega, 0) + z(\omega, 0) t \int_{0}^{\infty} \psi(\omega, \xi) d\xi. \quad (64) $$

At the second iteration, $z(\omega, t)$ is approximated by $z_1(\omega, t)$. So, we obtain

$$ z_2(\omega, t) = z(\omega, 0) + z(\omega, 0)t \int_{0}^{\infty} \psi(\omega, \xi) d\xi + z(\omega, 0) \frac{t^2}{2} \left[ \int_{0}^{\infty} \psi(\omega, \xi) d\xi \right]^2 \quad (65) $$

Additionally, at iteration $k$, we obtain

$$ z_k(\omega, t) = \sum_{j=0}^{k} \frac{t^j}{j!} \left[ \int_{0}^{\infty} \psi(\omega, \xi) d\xi \right]^j z(\omega, 0) \quad (66) $$
$$ \text{Thus, } z(\omega, t) = \lim_{k \to \infty} \sum_{j=0}^{k} \frac{t^j}{j!} \left[ \int_{0}^{\infty} \psi(\omega, \xi) d\xi \right]^j z(\omega, 0), \quad (67) $$

$$ \text{i.e., } z(\omega, t) = \left[ \sum_{k=0}^{\infty} \frac{t^k}{k!} \left[ \int_{0}^{\infty} \psi(\omega, \xi) d\xi \right]^k \right] z(\omega, 0) \quad (68) $$

Notice that, in the integer order case, $\frac{dx(t)}{dt} = ax(t)$, so $\psi(\omega, \xi) \equiv a$.

$$ \text{Then, } x(t) = \left[ \sum_{k=0}^{\infty} a^k \frac{t^k}{k!} \right] x(0) = e^{at} x(0). \quad (69) $$

Thus, in the distributed case, we can define the distributed exponential function:

$$ \phi(t) = \exp \left( t \int_{0}^{\infty} \psi(\omega, \xi) d\xi \right) = \sum_{k=0}^{\infty} \frac{t^k}{k!} \left[ \int_{0}^{\infty} \psi(\omega, \xi) d\xi \right]^k \quad (70) $$

and

$$ z(\omega, t) = \exp \left( t \int_{0}^{\infty} \psi(\omega, \xi) d\xi \right) z(\omega, 0) = \phi(t) z(\omega, 0) \quad \omega \in [0, \infty), \quad (71) $$

which is the distributed generalization of Equation (69).

Then, we can generalize the distributed exponential to the linear multidimensional non-commensurate order case:

$$ D^{\underline{n}}(\underline{x}(t)) = A\underline{x}(t) \dim(\underline{x}(t)) = N, \quad (72) $$

where $\underline{n}^T = [n_1 \dots n_i \dots n_N]$ $0 < n_i \le 1$.

Moreover, $0 < n_i \le 1$ means that the FDS system may include integer-order derivatives, since, with real systems, the dynamics are caused either by integer-order or fractional-order derivatives.

So, (72) corresponds to the integer-order frequency-distributed differential system:

$$ \left\{ \begin{array}{l} \frac{\partial \underline{z}(\omega, t)}{\partial t} = -\omega \underline{z}(\omega, t) + A \int_{0}^{\infty} [\mu_{\underline{n}}(\underline{\xi})] \underline{z}(\underline{\xi}, t) d\underline{\xi} \quad \dim(\underline{z}(t)) = N \\ \underline{x}(t) = \int_{0}^{\infty} [\mu_{\underline{n}}(\underline{\xi})] \underline{z}(\underline{\xi}, t) d\underline{\xi} \end{array} \right. \quad (73) $$

with the initial condition $\underline{z}(\omega, 0) \quad \forall \omega \in [0, \infty)$,

$$ \text{and } [\mu_{\underline{n}}(\underline{\xi})] = \begin{bmatrix} \mu_{n_1}(\underline{\xi}) & & 0 \\ & \mu_{n_i}(\underline{\xi}) & \\ 0 & & \mu_{n_N}(\underline{\xi}) \end{bmatrix} \quad (74) $$

Let us define the matrix

$$ \Psi(\omega, \xi) = -\omega\delta(\xi - \omega)I + A[\mu_{\underline{n}}(\underline{\xi})] \quad (75) $$

where *I* is the Identity matrix with appropriate dimension.

Then, the differential system (73) can be expressed as

$$ \frac{\partial \underline{z}(\omega, t)}{\partial t} = \int_{0}^{\infty} \Psi(\omega, \xi) \underline{z}(\underline{\xi}, t). \quad (76) $$
Therefore, the distributed exponential $\exp\left(t\int_{0}^{\infty}\Psi(\omega, \xi)d\xi\right)$ function is replaced by the distributed exponential matrix:

$$ \Phi(t) = \exp\left(t\int_{0}^{\infty}\Psi(\omega, \xi)d\xi\right) = \sum_{k=0}^{\infty} \frac{t^k}{k!} \left[\int_{0}^{\infty}\Psi(\omega, \xi)d\xi\right]^k, \quad (77) $$

and the solution of (76) is:

$$ \underline{z}(\omega, t) = \exp\left(t\int_{0}^{\infty}\Psi(\omega, \xi)d\xi\right) \underline{z}(\omega, 0) = \Phi(t)\underline{z}(\omega, 0) \quad \omega \in [0, \infty). \quad (78) $$

Notice that, contrary to the Mittag-Leffler matrix, the distributed exponential matrix shares the same semi-group properties (see chapter 2 of [54]) as its integer order analog, i.e.,

$$ \Phi(t, t_0) = \Phi(t, \tau)\Phi(\tau, t_0) \quad t_0 < \tau < t \quad (79) $$

Finally, the solution of the linear multidimensional non-commensurate order FDS (72) with the initial condition $\underline{z}(\omega, 0) \quad \forall \omega \in [0, \infty)$ is similar to the integer-order case (5), i.e.,

$$ \underline{z}(\omega, t) = \Phi(t)\underline{z}(\omega, 0) + \int_{0}^{t} \Phi(t - \tau)\underline{B}u(\tau)d\tau \quad \omega \in [0, \infty), \quad (80) $$

where $\Phi(t)$ is the distributed exponential transition matrix.

## 6.4. Computation of the Distributed Exponential Function

The interest of the distributed exponential function, either in scalar or matrix form, is essentially theoretical. Practically, in order to compute it, the distributed exponential function requires frequency discretization, i.e., a discretization of the continuous distributed integrator model:

$$ \left\{ \begin{aligned} \frac{\partial z(\omega, t)}{\partial t} &= -\omega z(\omega, t) + v(t) \quad \omega \in [0, \infty) \\ x(t) &= \int_{0}^{\infty} \mu(\omega)z(\omega, t)d\omega \end{aligned} \right. \quad (81) $$

Direct discretization is possible, but indirect discretization is more efficient [34,35]. It is based on the approximation of the frequency response of $\frac{1}{s^n} \quad 0 < n < 1$, on a frequency interval $[\omega_{min}, \omega_{max}]$ with J first-order cells, associated with an integer-order integrator at =0. So, the global model is composed of J + 1 cells (or modes) such as:

$$ \left\{ \begin{aligned} \frac{dz_j(t)}{dt} &= -\omega_j z_j(t) + v(t) \quad j = 0 \text{ to } J \\ x(t) &= \sum_{j=0}^{J} c_j z_j(t) \end{aligned} \right. \quad (82) $$

and the elementary system $D^n(x(t)) = ax(t)$ becomes:

$$ \left\{ \begin{aligned} \frac{dz_j(t)}{dt} &= -\omega_j z_j(t) + a \sum_{j=0}^{J} c_j z_j(t) \quad j = 0 \text{ to } J \\ x(t) &= \sum_{j=0}^{J} c_j z_j(t) \end{aligned} \right. \quad (83) $$

where $\omega_j$ varies from $\omega_0 = 0$ to $\omega_J$.
Let us define:

$$ 
\begin{gathered}
\underline{z}_I^T(t) = [z_0(t) \dots z_j(t) \dots z_J(t)] \quad \dim(\underline{z}_I) = J+1 \\
A_I = \begin{bmatrix} -\omega_0 & & \\ & -\omega_j & \\ & & -\omega_J \end{bmatrix} \quad \underline{B}_I^T = [1 \dots 1 \dots 1] \quad \underline{C}_I = [c_0 \dots c_j \dots c_J]
\end{gathered}
\quad (84) $$

Then, the previous system (83) is transformed into the following one

$$ \begin{cases} \frac{d\underline{z}_I(t)}{dt} = A_I \underline{z}_I(t) + a \underline{B}_I \underline{C}_I \underline{z}_I(t) = A_{\text{syst}} \underline{z}_I(t) & \text{with } A_{\text{syst}} = A_I + a \underline{B}_I \underline{C}_I \\ x(t) = \underline{C}_I \underline{z}_I(t) \end{cases} \quad (85) $$

Thus, the solution of system (85) with the initial condition $\underline{z}_I^T(0) = [z_0(0) \dots z_j(0) \dots z_J(0)]$ is expressed as

$$ \underline{Z}(t) = e^{A_{\text{syst}}t} \underline{Z}(0) = \phi_d(t) \underline{Z}(0), \quad (86) $$

where $\phi_d(t) = e^{A_{\text{syst}}t}$ is the matrix exponential corresponding to the frequency discretization of the distributed exponential function $\phi(t) = \exp\left(t \int_0^\infty \psi(\omega, \xi) d\xi\right)$.

Numerical examples of the distributed exponential function are available in chapter 9 volume 1 of [34].

## 7. Conclusions

In this paper, it has been proved, thanks to an elementary counter example, that the solution of the FDE initial value problem cannot be provided by the well-known Caputo derivative approach. The origin of this error relies on ignoring the dynamics of the fractional integrator, which has caused confusion between pseudo initial conditions of the Caputo derivative and the initialization function of the Riemann–Liouville integral.

Furthermore, it has been demonstrated that it is necessary to take into account two complementary models of the fractional integrator: the classical model used to express the pseudo state variable $x(t)$ is an input/output model, whereas the distributed model is an infinite-dimension state variable model, which is adapted to the formulation of free responses thanks to the initial values of its distributed state variable $z(\omega, t)$, as proved by the counter-example.

These two models of the fractional integrator generate two complementary approaches to the modeling of FDE and FDS one being focused on the pseudo-state variable, whereas the other permits to express internal transients linked to the distributed state variable. Two expressions of their free responses have been formulated in the linear case, one with the help of a convolution of the Mittag–Leffler function with the initialization function, and the other thanks to the definition of a distributed exponential function, which is a straightforward generalization of the integer-order case.

Beyond the misuse of fractional derivatives to solve the initial-value problem, the main conclusion of this paper is that any fractional-order system or one characterized by a long memory phenomenon is an infinite dimension system, whatever its input/output representation, linear or nonlinear.

Consequently, the use of an internal distributed representation is not an option but is necessary to correctly express initialization problems and dynamical transients.

**Author Contributions:** N.M. and J.-C.T. contributed equally to the paper. All authors have read and agreed to the published version of the manuscript.

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

**Conflicts of Interest:** The authors declare no conflict of interest.
## References

1. Coddington, E.A.; Levinson, N. *Theory of Ordinary Differential Equations*; Mc Graw Hill: New York, NY, USA, 1955.

2. Tenenbaum, M.; Pollard, H. *Ordinary Differential Equations*; Harper and Row: New York, NY, USA, 1963.

3. Podlubny, I. *Fractional Differential Equations*; Academic Press: San Diego, CA, USA, 1999.

4. Diethelm, K. *The Analysis of Fractional Differential Equations*; Lecture Notes in Mathematics; Springer: Berlin/Heidelberg, Germany, 2010.

5. Kilbas, A.A.; Srivastava, H.M.; Trujillo, J.J. *Theory and Applications of Fractional Differential Equations*; Elsevier: Amsterdam, The Netherlands, 2006.

6. Caputo, M. *Elasticita e Dissipazione*; Zanichelli: Bologna, Italy, 1969.

7. Fukunaga, M.; Shimizu, N. Role of pre-histories in the initial value problems of fractional viscoelastic equations. *Nonlinear Dyn.* **2004**, *38*, 207–220. [CrossRef]

8. Du, M.; Wang, Z. Initialized fractional differential equations with Riemann-Liouville fractional order derivative. In Proceedings of the ENOC 2011 Conference, Rome, Italy, 24–29 July 2011.

9. Du, M.; Wang, Z. Correcting the initialization of models with fractional derivatives via history dependent conditions. *Acta Mech. Sin.* **2016**, *32*, 320–325. [CrossRef]

10. Hartley, T.T.; Lorenzo, C.F. The error incurred in using the Caputo derivative Laplace transform. In Proceedings of the ASME IDET-CIE Conferences, San Diego, CA, USA, 30 August–2 September 2009.

11. Sabatier, J.; Merveillaut, M.; Malti, R.; Oustaloup, A. How to impose physically coherent initial conditions to a fractional system? *Commun. Non Linear Sci. Numer. Simul.* **2010**, *15*, 1318–1326. [CrossRef]

12. Ortigueira, M.D. *Fractional Calculus for Scientists and Engineers*; Springer Science: New York, NY, USA, 2011.

13. Trigeassou, J.C.; Maamri, N.; Oustaloup, A. The Caputo derivative and the infinite state approach. In Proceedings of the 6th Workshop on Fractional Differentiation and its Applications, Grenoble, France, 4–6 February 2013.

14. Sabatier, J.; Farges, C. Comments on the description and initialization of fractional partial differential equations using Riemann-Liouville's and Caputo's definitions. *J. Comput. Appl. Math.* **2018**, *339*, 30–39. [CrossRef]

15. Sabatier, J.; Farges, C. Initial value problems should not be associated to fractional model descriptions whatever the derivative definition used. *AIMS Math.* **2021**, *6*, 11318–11329. [CrossRef]

16. Lorenzo, C.F.; Hartley, T.T. Initialization in fractional order systems. In Proceedings of the European Control Conference (ECC’01), Porto, Portugal, 4–7 September 2001; pp. 1471–1476.

17. Hartley, T.T.; Lorenzo, C.F. The initialization response of linear fractional order system with constant History function. In Proceedings of the 2009 ASME/IDETC, San Diego, CA, USA, 30 August–2 September 2009.

18. Hartley, T.T.; Lorenzo, C.F. The initialization response of linear fractional order system with ramp History function. In Proceedings of the 2009 ASME/IDETC, San Diego, CA, USA, 30 August–2 September 2009.

19. Lorenzo, C.F.; Hartley, T.T. Initialization of fractional differential equation. *ASME J. Comput. Nonlinear Dyn.* **2007**, *4806*, 1341–1347.

20. Lorenzo, C.; Hartley, T. Time-varying initialization and Laplace transform of the Caputo derivative: With order between zero and one. In Proceedings of the IDETC/CIE FDTA’2011 Conference, Washington DC, USA, 14–17 August 2011.

21. Hartley, T.T.; Lorenzo, C.F.; Trigeassou, J.C.; Maamri, N. Equivalence of history function based and infinite dimensional state initializations for fractional order operators. *ASME J. Comput. Nonlinear Dyn.* **2013**, *8*, 041014. [CrossRef]

22. Trigeassou, J.C.; Maamri, N.; Oustaloup, A. The infinite state approach: Origin and necessity. *Comput. Math. Appl.* **2013**, *66*, 892–907. [CrossRef]

23. Tari, M.; Maamri, N.; Trigeassou, J.C. Initial conditions and initialization of fractional systems. *ASME J. Comput. Nonlinear Dyn.* **2016**, *11*, 041014. [CrossRef]

24. Maamri, N.; Tari, M.; Trigeassou, J.C. Improved initialization of fractional order systems. In Proceedings of the 20th World IFAC Congress, Toulouse, France, 14 July 2017; pp. 8567–8573.

25. Khalil, R.; Yousef, A.; Sababdeh, M. A new definition of fractional derivative. *J. Comput. Appl. Math.* **2014**, *264*, 65–70. [CrossRef]

26. Caputo, M.; Fabrizio, M. A new definition of fractional derivative without singular kernel. *Prog. Fract. Differ. Appl.* **2015**, *1*, 73–85.

27. Atanackovic, T.M.; Pilipovic, S.; Zorica, D. Properties of the Caputo-Fabrizio fractional derivative and its distributional setting. *Fract. Calc. Appl. Anal.* **2018**, *21*, 29–44. [CrossRef]

28. Martinez, F.; Othman Mohammed, P.; Napoles Valdes, J.E. Non-conformable fractional Laplace transform. *Kragujev. J. Math.* **2022**, *46*, 341–354. [CrossRef]

29. Ortigueira, M.D.; Machado, J.T. A critical analysis of the Caputo-Fabrizio operator. *Commun. Nonlinear Sci Numer. Simulat.* **2018**, *59*, 608–611. [CrossRef]

30. Tarasov, V.E. No nonlocality, no fractional derivative. *Commun. Nonlinear Sci Numer. Simulat.* **2018**, *62*, 157–163. [CrossRef]

31. Montseny, G. Diffusive Representation of Pseudo Differential Time Operators. *ESAIM Proc.* **1998**, *5*, 159–175. [CrossRef]

32. Heleschewitz, D.; Matignon, D. Diffusive realizations of fractional integro-differential operators: Structural analysis under approximation. In Proceedings of the Conference IFAC, System, Structure and Control, Nantes, France, 8–10 July 1998; Volume 2, pp. 243–248.

33. Trigeassou, J.C.; Maamri, N.; Sabatier, J.; Oustaloup, A. Transients of fractional order integrator and derivatives. In *Special Issue: “Fractional Systems and Signals” of Signal, Image and Video Processing*; Springer: Berlin/Heidelberg, Germany, 2012; Volume 6, pp. 359–372.
34. Trigeassou, J.C.; Maamri, N. *Analysis, Modeling and Stability of Fractional Order Differential Systems: The Infinite State Approach*; John Wiley and Sons: Hoboken, NJ, USA, 2019; Volumes 1 and 2.

35. Trigeassou, J.C.; Poinot, T.; Lin, J.; Oustaloup, A.; Levron, F. Modeling and identification of a non integer order system. In Proceedings of the ECC'99 European Control Conference, Karlsruhe, Germany, 31 August–3 September 1999.

36. Lin, J.; Poinot, T.; Trigeassou, J.C. Parameter estimation of fractional systems: Application to the modeling of a lead-acid battery. In Proceedings of the 12th IFAC Symposium on System Identification, Santa Barbara, CA, USA, 21–23 June 2000.

37. Benchellal, A.; Poinot, T.; Trigeassou, J.C. Approximation and identification of diffusive interface by fractional models. *Signal Process.* **2006**, *86*, 2712–2727. [CrossRef]

38. Trigeassou, J.C.; Maamri, N.; Sabatier, J.; Oustaloup, A. State variables and transients of fractional order differential systems. *Comput. Math. Appl.* **2012**, *64*, 3117–3140.

39. Trigeassou, J.C.; Maamri, N.; Sabatier, J.; Oustaloup, A. A Lyapunov approach to the stability of fractional differential equations. *Signal Process.* **2011**, *91*, 437–445. [CrossRef]

40. Du, B.; Wei, Y.; Liang, S.; Wang, Y. Estimation of exact initial states of fractional order systems. *Nonlinear Dyn.* **2016**, *86*, 2061–2070. [CrossRef]

41. Yuan, J.; Zhang, Y.; Liu, J.; Shi, B. Equivalence of initialized fractional integrals and the diffusive model. *ASME J. Comput. Nonlinear Dyn.* **2018**, *13*, 034501. [CrossRef]

42. Zhao, Y.; Wei, Y.; Chen, Y.; Wang, Y. A new look at the fractional initial value problem: The aberration phenomenon. *ASME J. Comput. Nonlinear Dyn.* **2018**, *13*, 121004. [CrossRef]

43. Yuan, J.; Shi, B.; Ji, W. Adaptive sliding mode control of a novel class of fractional chaotic systems. *Adv. Math. Phys.* **2013**, *13*. [CrossRef]

44. Wang, B.; Ding, J.; Wu, F.; Zhu, D. Robust finite time control of fractional order nonlinear systems via frequency distributed model. *Nonlinear Dyn.* **2016**, *85*, 2133–2142. [CrossRef]

45. Chen, Y.; Wei, Y.; Zhou, X.; Wang, Y. Stability for nonlinear fractional order systems: An indirect approach. *Nonlinear Dyn.* **2017**, *89*, 1011–1018. [CrossRef]

46. Wei, Y.; Sheng, D.; Chen, Y.; Wang, Y. Fractional order chattering free robust adaptive backstepping control technique. *Nonlinear Dyn.* **2019**, *95*, 2383–2394. [CrossRef]

47. Hinze, M.; Schmidt, A.; Leine, R.L. Numerical solution of fractional order ordinary differential equations using the reformulated infinite state representation. *Fract. Calc. Appl. Anal.* **2019**, *22*, 1321–1350. [CrossRef]

48. Hinze, M.; Schmidt, A.; Leine, R.L. The direct method of Lyapunov for nonlinear dynamical systems with fractional damping. *Nonlinear Dyn.* **2020**, *102*, 2017–2037. [CrossRef]

49. Maamri, N.; Trigeassou, J.C. Integration of fractional differential equations without fractional derivatives. In Proceedings of the ICSC20 Conference, Caen France, 24–26 November 2021.

50. Lindelöf, E. Sur l'application de la méthode des approximations successives aux équations différentielles ordinaires du premier ordre. *C.R.A.C.* **1894**, *116*, 454–457.

51. Picard, E. Mémoire sur la théorie des équations aux dérivées partielles et la méthodes des approximations successives. *J. Math. Pures Appl.* **1890**, *4*, 145–210.

52. Kailath, T. *Linear Systems*; Prentice Hall Inc.: Englewood Cliffs, NJ, USA, 1980.

53. Monje, C.A.; Chen, Y.Q.; Vinagre, B.M.; Xue, D.; Feliu, V. *Fractional Order Systems and Control*; Springer: London, UK, 2010.

54. Curtain, R.F.; Zwart, H.J. *An Introduction to Infinite Dimensional Linear Systems Theory*; Springer: New York, NY, USA, 1995.