Novel Techniques for a Verified Simulation of Fractional-Order Differential Equations

Abstract

Verified simulation techniques have been investigated intensively by researchers who are dealing with ordinary and partial differential equations. Tasks that have been considered in this context are the solution to initial value problems and boundary value problems, parameter identification, as well as the solution of optimal control problems in cases in which bounded uncertainty in parameters and initial conditions are present. In contrast to system models with integer-order derivatives, fractional-order models have not yet gained the same attention if verified solution techniques are desired. In general, verified simulation techniques rely on interval methods, zonotopes, or Taylor model arithmetic and allow for computing guaranteed outer enclosures of the sets of solutions. As such, not only the influence of uncertain but bounded parameters can be accounted for in a guaranteed way. In addition, also round-off and (temporal) truncation errors that inevitably occur in numerical software implementations can be considered in a rigorous manner. This paper presents novel iterative and series-based solution approaches for the case of initial value problems to fractional-order system models, which will form the basic building block for implementing state estimation schemes in continuous-discrete settings, where the system dynamics is assumed as being continuous but measurements are only available at specific discrete sampling instants.

1. Introduction

As mentioned in the abstract of this paper, verified (often interval-based) simulation techniques are widely known for sets of ordinary differential equation when computing the solution to initial value problems [1]. Throughout this paper, we employ the widely used notion of verification from the community of interval analysis, where a verified method deals with solution procedures that by construction contain all possible results in a guaranteed manner. Validation instead means that a mathematical model for a real-life application is consistent with experimental observations [2].

The corresponding research activities resulted in many (often freely available) software implementations such as AWA [3,4], CAPD [5], Cosy-VI [6,7], DynIBEX [8,9], ValEncIA-IVP [10,11], VNODE-LP [12,13], and VSPODE [14]. Besides using intervals and Taylor models as the underlying solution representations, also zonotopic descriptions were investigated in these works [15]. Concerning the approaches implemented in AWA and Cosy-VI, also an up-to-date Matlab implementation should be pointed out [16]. The possible use of such verified simulation techniques ranges from parameter identification for complex dynamic systems [17,18], over the design and verification of optimal control procedures [19], reachability analysis [20,21], actuator fault analysis [22], observability analysis [23] to generalizations for systems with piecewise defined right-hand sides and other hybrid system models [24,25,26,27].

Although the dynamics of many system models in engineering applications can be analyzed with the methods mentioned above, the verified analysis of fractional-order systems is an emerging direction of research that is not yet fully solved. Fractional-order differential equations are characterized by the fact that not only integer orders of time derivatives appear in the state equations but also corresponding non-integer values are take into consideration [28,29,30,31]. If practical applications in (control) engineering are concerned, the applications for such dynamic system representations cover topics related to

• Battery modeling and state-of-charge estimation [32,33,34];
• Fractional-order modeling and control of robotic manipulators aiming at enhanced loop-shaping [35,36];
• Modeling the formation of multi-robot systems, where each agent has integer-order dynamics but the overall behavior may result in a fractional-order behavior [37];
• Advanced modeling of visco-elastic damping [38];
• The generalization of classical PID and integral sliding-mode controllers [39,40,41,42];
• Modeling phenomena that cannot be explained fully by classical stress–strain relations (Newton’s vs. Hooke’s law) or diffusion vs. wave propagation phenomena [43,44].

If set-valued simulation techniques for fractional-order system models are concerned, especially the Mittag–Leffler function approach for computing enclosures for systems with stable dynamics should be pointed out [45,46]. It is based on generalizing the exponential enclosure technique [47,48]. This technique was originally developed by A. Rauh and E. Auer for integer-order dynamics. It is now generalized toward the use in fractional system models, where the Mittag–Leffler function is an extension of the exponential function. Both the exponential and Mittag–Leffler functions naturally represent closed-form solutions of linear either integer-order or fractional-order system models that are specified by linear state equations [29,49,50,51,52].

In this paper, Section 2 gives an overview of preliminaries and selected state-of-the-art techniques for fractional-order differential equations. These techniques are extended in Section 3 toward novel interval-based routines. The specific properties of these novel approaches are discussed in Section 4 for selected linear and nonlinear case studies. Finally, conclusions and an outlook on future work are given in Section 5.

2. Preliminaries and State-of-the-Art

2.1. System Models under Consideration

In this paper, we consider the simulation of in general nonlinear, commensurate fractional-order autonomous differential equations

$${\mathbf{x}}^{\left(\nu \right)}\left(t\right)=\mathbf{f}\left(\mathbf{x}\left(t\right)\right),\mathbf{f}:{\mathbb{R}}^n\mapsto {\mathbb{R}}^n,$$

(1)

where the right-hand side $\mathbf{f}\left(\mathbf{x}\left(t\right)\right)$ is assumed to be given by a continuous function.

Furthermore, assume that $0<\nu \le 1$ holds. Then, in the sense of a fractional-order system model of Caputo type [28,29], initial conditions are specified at the time instant $t=0$. Uncertainty in the initial states $\mathbf{x}\left(0\right)$ is described in terms of the interval vector

$$\mathbf{x}\left(0\right)\in \left[\mathbf{x}\right]\left(0\right):=\left[\underset{¯}{\mathbf{x}}\left(0\right);\overline{\mathbf{x}}\left(0\right)\right],$$

(2)

for which ${\underset{¯}{x}}_i\left(0\right)\le {\overline{x}}_i\left(0\right)$ holds for each vector component $i\in \{1,\dots ,n\}$.

2.2. Grünwald–Letnikov Approximation of Fractional Derivatives

One of the most classical approximations for the solution of fractional-order differential equations is the use of the Grünwald–Letnikov definition of fractional derivatives and, on its basis, the representation of the solution in terms of the infinite series

$$\mathbf{x}\left(t_{k+1}\right)=\nu \mathbf{I}·\mathbf{x}\left(t_k\right)+\Delta T_k·\mathbf{f}\left(\mathbf{x}\left(t_k\right)\right)-\sum _{i=2}^{\infty }{\left(-1\right)}^i·\left(\genfrac{}{}{0pt}{}{\nu }{i}\right)·\mathbf{x}\left(t_{k+1-i}\right)$$

(3)

with the step size $\Delta T_k=t_{k+1}-t_k$.

Here, the term $\left(\genfrac{}{}{0pt}{}{\nu }{i}\right)$ represents the Newton binomial coefficient

$$\left(\genfrac{}{}{0pt}{}{\nu }{i}\right)=\frac{\Gamma \left(\nu +1\right)}{\Gamma \left(i+1\right)·\Gamma \left(\nu -i+1\right)},$$

(4)

where the gamma function $\Gamma \left(\nu \right)$ serves as a generalization of the factorial for the case of real-valued arguments according to

$$\Gamma \left(\nu \right)=\underset{0}{\overset{\infty }{\int }}{\xi }^{\nu -1}{e}^{-\xi }\mathrm{d}\xi .$$

(5)

The severe disadvantage of this series definition in practical applications is its low rate of convergence. Hence, a large number of summands need to be accounted for in the approximation of the sum in Equation (3) to get sufficiently close to the true solution of the system model (1).

Note, in a verified setting, the truncated part of the infinite series cannot be ignored. In principle, however, it could be overapproximated with the help of the binomial series, cf. [53,54], dating back to Isaac Newton (around 1665) according to

$${\left(1+x\right)}^{\nu }=\sum _{i=0}^{\infty }\left(\genfrac{}{}{0pt}{}{\nu }{i}\right)x^i,$$

(6)

or, more generally,

$${\left(x+y\right)}^{\nu }=\sum _{i=0}^{\infty }\left(\genfrac{}{}{0pt}{}{\nu }{i}\right)x^{\nu -i}{y}^i\mathrm{with}\left|x\right|>\left|y\right|,$$

(7)

which would then have to be split up into a finite sum and into a remainder, where the values of x can be set to fixed values or interval bounds (often the initial condition).

Due to the fact that (3) converges too slowly in practice and that bounding the remainder term as sketched above is not trivial and especially prone to overestimation in an interval setting, this option will not be considered further in this paper.

2.3. Single Time Step Picard Iteration

Further, quite recently investigated solution techniques for fractional-order differential equations are given by a Picard iteration procedure. As summarized in the two following theorems, taken from the cited literature without explicitly stating the proof in this paper, either point-valued approximations of the solution (Theorem 1) or interval-valued bounding boxes enclosing all solutions over a finitely long time span (Theorem 2) can be distinguished.

The following section of this paper will make use of a combination of the ideas of both theorems in the sense that a novel formulation which is close to Theorem 2 will be turned into an interval-valued iteration that produces bounding boxes of states for multiple subsequent time slices.

Theorem 1

([55,56] Integral formulation of Picard iterations for fractional-order differential equations). Let $\mathbf{f}\left(\mathbf{x}\left(t\right)\right)$ be a continuous Lipschitzian function on a bounded state and time domain. The solution to time-invariant fractional-order differential equations at the point of time $T>0$ can be computed iteratively according to the fixed-point iteration

$${\mathbf{x}}^{\langle \kappa +1\rangle }\left(T\right):=\mathbf{x}\left(0\right)+\frac{1}{\Gamma \left(\nu \right)}·\underset{0}{\overset{T}{\int }}{\left(T-s\right)}^{\nu -1}·\mathbf{f}\left({\mathbf{x}}^{\langle \kappa \rangle }\left(s\right)\right)\mathrm{d}s,\kappa \in {\mathbb{N}}_0,$$

(8)

with the initialization ${\mathbf{x}}^{\langle 0\rangle }:=\mathbf{x}\left(0\right)$ at the iteration step $\kappa =0$.

The convergence properties of the iteration in Theorem 1 are analyzed in detail in [56].

Theorem 2

([55] Interval-based, single time step Picard iterations for fractional-order differential equations). The iteration in Theorem 1 generalizes to interval bounded initial conditions $\mathbf{x}\left(0\right)\in \left[\mathbf{x}\right]\left(0\right)=\left[{\mathbf{x}}_0\right]$ according to

$${\left[\mathbf{x}\right]}^{\langle \kappa +1\rangle }:=\left[{\mathbf{x}}_0\right]+\frac{1}{\Gamma (\nu +1)}·\left[0;T^{\nu }\right]·\mathbf{f}\left({\left[\mathbf{x}\right]}^{\langle \kappa \rangle }\right),\kappa \in {\mathbb{N}}_0,$$

(9)

where convergence requires ${\left[\mathbf{x}\right]}^{\langle \kappa +1\rangle }\subseteq {\left[\mathbf{x}\right]}^{\langle \kappa \rangle }$, leading to $\mathbf{x}\left(t\right)\in {\left[\mathbf{x}\right]}^{\langle \kappa +1\rangle }$ for all $t\in \left[0;T\right]$.

As mentioned above, the numerical evaluation of Theorem 2 provides time-invariant bounds ${\left[\mathbf{x}\right]}^{\langle \kappa +1\rangle }$ containing all possible solutions $\mathbf{x}\left(t\right)$ over the time interval $t\in \left[0;T\right]$. This fact limits the practical applicability to quite small values of T, which becomes obvious if, for example, $\nu =1$ was considered. Then, (9) would be nothing else than an explicit Euler evaluation of the solution with its corresponding well-known restrictions to small admissible values for T.

Remark 1.

Throughout this paper, the evaluation of functions for interval arguments such as in (9) denotes an arbitrary interval extension of the corresponding expressions. To reduce pessimism, these interval extensions should result in tight interval boxes without significant overestimation. This can be obtained by means of the exploitation of higher-order enclosure techniques or by making use of certain monotonicity properties of the functions to be evaluated (cf. [57,58,59]).

2.4. Special Case: “Closed-Form” Solutions

Consider—as a special linear case—the scalar fractional-order differential equation

$$x^{\left(\nu \right)}\left(t\right)=\lambda ·x\left(t\right)\mathrm{with}x\left(0\right)=x_0$$

(10)

and the parameter $\lambda \in \mathbb{R}$. According to [50,60], the exact solution of this linear system model of Caputo type is given by the expression

$$x\left(t\right)=E_{\nu ,1}\left(\lambda t^{\nu }\right)·x\left(0\right),$$

(11)

where $E_{\nu ,\beta }\left(\zeta \right)$ is the two-parameter Mittag–Leffler function

$$E_{\nu ,\beta }\left(\zeta \right)=\sum _{i=0}^{\infty }\frac{{\zeta }^i}{\Gamma \left(\nu i+\beta \right)}$$

(12)

with the general argument $\zeta \in \mathbb{C}$, the gamma function $\Gamma \left(\nu i+\beta \right)$ as defined in (5), the derivative order $\nu$, and $\beta =1$ according to (11).

An accurate floating-point evaluation of the function (12) becomes possible by the Matlab implementation by R. Garrappa [51,61]. For the sake of completeness, although the procedure cannot be turned directly into a verified implementation using interval analysis, it should be pointed out at this position that numerical approximations of the solution of fractional-order differential equations can be obtained effectively with a predictor–corrector method implemented in the Matlab routine fde12 [62] published also by R. Garrappa.

3. Verified Simulation Techniques

In this section, three novel procedures are introduced for a guaranteed, interval-based state enclosure technique that is applicable to fractional-order differential equations. Note, in the limit case $\nu =1$, they turn into results that are well-known from the field of verified simulation of classical sets of ordinary differential equations.

3.1. Mittag–Leffler Type State Enclosures for Fractional-Order Differential Equations

As stated in Section 2.4, closed-form solutions of linear fractional-order differential equations can be stated in terms of Mittag–Leffler functions. Therefore, the first investigated solution technique generalizes this solution representation towards the nonlinear case by computing parameters $\lambda$ in (11) that are no longer point-valued but defined in terms of intervals with a finitely large width. In such a way, the solution representation according to Definition 1 can be employed for a verified simulation of fractional-order models that may include nonlinearities and bounded uncertainty in parameters and initial conditions. This solution approach makes use of an interval extension of Mittag–Leffler functions for real-valued arguments that was first published in [63,64].

Definition 1

(Mittag–Leffler type state enclosure). The time-dependent Mittag–Leffler type enclosure function

$${\mathbf{x}}^{*}\left(t\right)\in \left[{\mathbf{x}}_e\right]\left(t\right)={\mathbf{E}}_{\nu ,1}\left(\left[\Lambda \right]·t^{\nu }\right)·\left[{\mathbf{x}}_e\right]\left(0\right),\left[{\mathbf{x}}_e\right]\left(0\right)=\left[{\mathbf{x}}_0\right]$$

(13)

with the diagonal parameter matrix $\left[\Lambda \right]:=\mathrm{diag}\left\{\left[{\lambda }_i\right]\right\}$, $i\in \{1,\dots ,n\}$, is denoted as a verified Mittag–Leffler type state enclosure for the system model (1) with (2) if it is determined according to Theorem 3.

Theorem 3

([45,46,64] Iteration for Mittag–Leffler type enclosures). The Mittag–Leffler type state enclosure (13) contains all reachable states ${\mathbf{x}}^{*}\left(T\right)$ at the point of time $t=T>0$ according to

$${\mathbf{x}}^{*}\left(T\right)\in \left[{\mathbf{x}}_e\right]\left(T\right)={\mathbf{E}}_{\nu ,1}\left(\left[\Lambda \right]·T^{\nu }\right)·\left[{\mathbf{x}}_e\right]\left(0\right)$$

(14)

with certainty if the diagonal elements of $\left[\Lambda \right]$ are set to the outcome of the converging iteration

$${\left[{\lambda }_i\right]}^{\langle \kappa +1\rangle }:=\frac{f_i\left({\mathbf{E}}_{\nu ,1}\left({\left[\Lambda \right]}^{\langle \kappa \rangle }·{\left[t\right]}^{\nu }\right)·\left[{\mathbf{x}}_e\right]\left(0\right)\right)}{E_{\nu ,1}\left({\left[{\lambda }_i\right]}^{\langle \kappa \rangle }·{\left[t\right]}^{\nu }\right)·\left[x_{e,i}\right]\left(0\right)},$$

(15)

$i\in \{1,\dots ,n\}$, with the prediction horizon $\left[t\right]=\left[0;T\right]$. The value $x_i^{*}=0$ must not belong to the solution for any vector component $i\in \{1,\dots ,n\}$.

Proof. 

According to [45,46], Picard iteration procedures for fractional-order systems cannot only be specified in integral form as in Theorems 1 and 2 but also in the equivalent differential form

$${\left[{\mathbf{x}}^{\left(\nu \right)}\right]}^{\langle \kappa \rangle }\left(\left[0;T\right]\right)\supset {\left[{\mathbf{x}}^{\left(\nu \right)}\right]}^{\langle \kappa +1\rangle }\left(\left[0;T\right]\right)=\mathbf{f}\left({\left[{\mathbf{x}}^{\left(\nu \right)}\right]}^{\langle \kappa \rangle }\left(\left[0;T\right]\right)\right),$$

(16)

where ${\left[{\mathbf{x}}^{\left(\nu \right)}\right]}^{\langle \kappa \rangle }\left(\left[0;T\right]\right)$ is an interval extension of the time derivative of order $\nu$ of the time-dependent function ${\left[\mathbf{x}\right]}^{\langle \kappa \rangle }\left(t\right)$ over the time interval $t\in \left[t\right]=\left[0;T\right]$ in the iteration step $\kappa$. Substituting the Mittag–Leffler type enclosure of Definition 1 into (16) yields the expression

$$\begin{array}{cc}\hfill {\mathbf{x}}^{\left(\nu \right)}\left(t\right)& \in \left({\left[\Lambda \right]}^{\langle \kappa +1\rangle }\right)·{\mathbf{E}}_{\nu ,1}\left({\left[\Lambda \right]}^{\langle \kappa +1\rangle }·t^{\nu }\right)·\left[{\mathbf{x}}_e\right]\left(0\right)\hfill \\ \hfill & =\mathbf{f}\left({\mathbf{E}}_{\nu ,1}\left({\left[\Lambda \right]}^{\langle \kappa \rangle }·t^{\nu }\right)·\left[{\mathbf{x}}_e\right]\left(0\right)\right)\mathrm{for}\mathrm{all}t\in \left[t\right].\hfill \end{array}$$

(17)

Overapproximating the Mittag–Leffler type state enclosure ${\mathbf{E}}_{\nu ,1}\left({\left[\Lambda \right]}^{\langle \kappa +1\rangle }·t^{\nu }\right)·\left[{\mathbf{x}}_e\right]\left(0\right)$ in the iteration step $\kappa +1$ on the first line of (17) by the enclosure ${\left[{\mathbf{x}}_e\right]}^{\langle \kappa \rangle }\left(\left[t\right]\right)$ in the case of convergence, i.e., using the relation

$${\left[{\mathbf{x}}_e\right]}^{\langle \kappa +1\rangle }\left(\left[t\right]\right)\subset {\left[{\mathbf{x}}_e\right]}^{\langle \kappa \rangle }\left(\left[t\right]\right),$$

(18)

which is based on the assumption

$${\left[{\lambda }_i\right]}^{\langle \kappa +1\rangle }\subset {\left[{\lambda }_i\right]}^{\langle \kappa \rangle },$$

(19)

leads to

$$\begin{array}{cc}\hfill & \mathrm{diag}\left\{{\left[{\tilde{\lambda }}_i\right]}^{\langle \kappa +1\rangle }\right\}·{\left[{\mathbf{x}}_e\right]}^{\langle \kappa \rangle }\left(\left[t\right]\right)=\mathbf{f}\left({\left[{\mathbf{x}}_e\right]}^{\langle \kappa \rangle }\left(\left[t\right]\right)\right),\hfill \end{array}$$

(20)

where

$$\mathrm{diag}\left\{{\left[{\tilde{\lambda }}_i\right]}^{\langle \kappa +1\rangle }\right\}\supseteq \mathrm{diag}\left\{{\left[{\lambda }_i\right]}^{\langle \kappa +1\rangle }\right\}.$$

(21)

The relation (21) then holds naturally due to inclusion monotonicity of interval analysis [57,58,59].

Solving the expression (20) for ${\left[{\tilde{\lambda }}_i\right]}^{\langle \kappa +1\rangle }$ with subsequently renaming this parameter into ${\left[{\lambda }_i\right]}^{\langle \kappa +1\rangle }$ completes the proof of Theorem 3. □

Remark 2.

A typical initialization of ${\left[{\lambda }_i\right]}^{\langle 0\rangle }$ is given by the real parts of the eigenvalues of the Jacobian of $\mathbf{f}$ at the midpoint of $\left[{\mathbf{x}}_e\right]\left(0\right)$, inflated to small intervals with non-zero diameter. The inflation — currently implemented as the interval hull over ${\left[{\lambda }_i\right]}^{\langle \kappa \rangle }$ and ${\left[{\lambda }_i\right]}^{\langle \kappa +1\rangle }$ — of these bounds is continued up to the point, where the condition (19) is satisfied.

Remark 3.

To achieve tight interval bounds, it is recommended to transform the state equations according to [45,46] into a coordinate frame in which they are decoupled as far as possible. For a time-invariant change of coordinates, the matrix of the eigenvectors of the Jacobian of $\mathbf{f}$ at the midpoint of $\left[{\mathbf{x}}_e\right]\left(0\right)$ can be used if all of them are real-valued, see also [47,48].

Corollary 1.

For quasi-linear state-space representations of fractional-order differential equations in the form

$${\mathbf{x}}^{\left(\nu \right)}\left(t\right)=\mathcal{𝓐}\left(\mathbf{x}\left(t\right)\right)·\mathbf{x}\left(t\right)\mathrm{with}0<\nu \le 1,$$

(22)

which replace the general dynamics (1) for systems with an equilibrium at $\mathbf{x}=\mathbf{0}$ and well-defined, finite entries in the state-dependent matrix $\mathcal{𝓐}\left(\mathbf{x}\left(t\right)\right)$ and considering the initial conditions (2), Theorem 3 simplifies to the iteration scheme

$$\begin{array}{cc}\hfill {\left[{\lambda }_i\right]}^{\langle \kappa +1\rangle }:=& a_{ii}\left({\left[{\mathbf{x}}_e\right]}^{\langle \kappa \rangle }\left(\left[t\right]\right)\right)+\sum _{\genfrac{}{}{0pt}{}{j=1}{j\ne i}}^n\left\{a_{ij}\left({\left[{\mathbf{x}}_e\right]}^{\langle \kappa \rangle }\left(\left[t\right]\right)\right)·\frac{E_{\nu ,1}\left({\left[{\lambda }_j\right]}^{\langle \kappa \rangle }·{\left[t\right]}^{\nu }\right)}{E_{\nu ,1}\left({\left[{\lambda }_i\right]}^{\langle \kappa \rangle }·{\left[t\right]}^{\nu }\right)}·\frac{\left[x_{e,j}\right]\left(0\right)}{\left[x_{e,i}\right]\left(0\right)}\right\}.\hfill \end{array}$$

(23)

Remark 4.

Note that for $\nu \ne 1$ the quotient of two Mittag–Leffler functions in (23) can usually not be simplified further. Overestimation in the interval evaluation due to multiple dependencies on common interval parameters can be reduced by exploiting the monotonicity properties for Mittag–Leffler functions that were analyzed in detail in [64].

Remark 5.

For the order $\nu =1$, the iteration formulas in Theorem 3 and Corollary 1 become identical to the exponential enclosure technique published in [47,48] due to $E_{1,1}\left(x\right)\equiv e^x$.

3.2. Multi Time Step Picard Iteration

When using Theorem 1 to compute approximations of the solution to the system model (1) with (2), the integration that is desired for the horizon $t\in \left[0;T\right]$ can be performed in intermediate steps for shorter ranges $t\le T$ according to

$${\mathbf{x}}^{\langle \kappa +1\rangle }\left(t\right):=\mathbf{x}\left(0\right)+\frac{1}{\Gamma \left(\nu \right)}·\underset{0}{\overset{t}{\int }}{\left(t-\tau \right)}^{\nu -1}·\mathbf{f}\left({\mathbf{x}}^{\langle \kappa \rangle }\left(\tau \right)\right)\mathrm{d}\tau .$$

(24)

This reformulation leads to the interval formulation of the Picard iteration on the temporal discretization mesh $0=t_0<t_1<\dots <t_n=T$ according to the following theorem.

Theorem 4

(Interval-valued multi time step Picard iteration). A guaranteed interval enclosure of the state vector $\mathbf{x}\left(t_n\right)$ is obtained by the iteration

$${\mathbf{x}}^{\langle \kappa +1\rangle }\left(t_n\right)\in \mathbf{x}\left(0\right)+\frac{1}{\Gamma \left(\nu \right)}·\sum _{i=0}^{n-1}{\left[{\mathbf{f}}^{\langle \kappa \rangle }\right]}_i^{i+1}·\frac{{\left(t_n-t_i\right)}^{\nu }-{(t_n-t_{i+1})}^{\nu }}{\nu }$$

(25)

with

$${\left[{\mathbf{f}}^{\langle \kappa \rangle }\right]}_i^{i+1}=\mathbf{f}\left({\left[\mathbf{x}\right]}^{\langle \kappa \rangle }\left(\left[t_i;t_{i+1}\right]\right)\right),$$

(26)

if the temporally piecewise constant defined state bounds satisfy the relation

$${\left[\mathbf{x}\right]}^{\langle \kappa +1\rangle }\left(\left[t_i;t_{i+1}\right]\right)\subset {\left[\mathbf{x}\right]}^{\langle \kappa \rangle }\left(\left[t_i;t_{i+1}\right]\right).$$

(27)

Proof. 

Using (24), Formula (8) is cast into the interval-valued Picard iteration

$${\mathbf{x}}^{\langle \kappa +1\rangle }\left(t_n\right)\in \mathbf{x}\left(0\right)+\frac{1}{\Gamma \left(\nu \right)}·\sum _{i=0}^{n-1}\left(\underset{t_i}{\overset{t_{i+1}}{\int }}{\left(t_n-\tau \right)}^{\nu -1}·\underset{\subseteq {\left[{\mathbf{f}}^{\langle \kappa \rangle }\right]}_i^{i+1}}{\underbrace{\mathbf{f}\left(\left[{\mathbf{x}}^{\langle \kappa \rangle }\right]\left(\tau \right)\right)}}\mathrm{d}\tau \right).$$

(28)

Assuming temporally piecewise constant defined interval bounds ${\left[{\mathbf{f}}^{\langle \kappa \rangle }\right]}_i^{i+1}$ as it was done, for example, in [65], the integral in (28) can be bounded by

$${\mathbf{x}}^{\langle \kappa +1\rangle }\left(t_n\right)\in \mathbf{x}\left(0\right)+\frac{1}{\Gamma \left(\nu \right)}·\sum _{i=0}^{n-1}\left({\left[{\mathbf{f}}^{\langle \kappa \rangle }\right]}_i^{i+1}·\underset{t_i}{\overset{t_{i+1}}{\int }}{\left(t_n-\tau \right)}^{\nu -1}\mathrm{d}\tau \right).$$

(29)

Solving the remaining integral expression in (29) in closed form completes the proof of Theorem 4. □

Remark 6.

For the integer-order case with $\nu =1$, this formulation becomes identical to the fundamental iteration scheme ofValEncIA-IVP, when setting the approximate solution introduced in [11] to zero.

3.3. Truncated Temporal Series Expansion

Theorem 5

(Truncated series enclosure for fractional-order differential equations). A truncated temporal series of length L with interval coefficients according to

$$\mathbf{x}\left(t\right)\in \sum _{i=0}^L\frac{\left[{\mathbf{x}}_{i·\nu }\right]·t^{i·\nu }}{\Gamma (1+i·\nu )}+\left[\mathbf{R}\right]$$

(30)

with the interval enclosure of the corresponding time derivative

$${\mathbf{x}}^{\left(\nu \right)}\left(t\right)\in \sum _{i=0}^{L-1}\frac{\left[{\mathbf{x}}_{(i+1)·\nu }\right]·t^{i·\nu }}{\Gamma (1+i·\nu )}+\left[{\mathbf{R}}^{\prime }\right],$$

(31)

where

$$\left[\mathbf{R}\right]=\left[-1;1\right]·\left[{\mathbf{R}}^{\prime }\right]·\frac{\left[0;t^{\nu }\right]}{\Gamma (1+\nu )},$$

(32)

is a guaranteed enclosure of all reachable states at $t>0$ if it satisfies the relation

$$\sum _{i=0}^{L-1}\frac{\left[{\mathbf{x}}_{(i+1)·\nu }\right]·t^{i·\nu }}{\Gamma (1+i·\nu )}+\left[{\mathbf{R}}^{\prime }\right]\subset \mathbf{f}\left(\sum _{i=0}^L\frac{\left[{\mathbf{x}}_{i·\nu }\right]·t^{i·\nu }}{\Gamma (1+i·\nu )}+\left[\mathbf{R}\right]\right)\mathrm{for}\mathrm{all}t\in \left[0;T\right].$$

(33)

Proof. 

The proof of this theorem consists of two parts. First, it is verified that differentiating (30) yields the expression (31). According to [28], the fractional-order derivative of the time-dependent polynomial ${\left(t-t_0\right)}^i$, $i\in \mathbb{N}$, is given by

$${}_{t_0}^{}{\mathcal{D}}_t^{\nu }{\left(t-t_0\right)}^i=\frac{\Gamma \left(1+i\right)}{\Gamma \left(1+i-\nu \right)}·{\left(t-t_0\right)}^{i-\nu }.$$

(34)

Here, $t_0$ is the point of time with respect to which the derivative is defined. Setting $t_0=0$ equal to the initial point of time at which the initial conditions are defined for the fractional-order system models in this paper, correctness of (31) is proven.

The remainder of the proof is a direct consequence of applying the differential formulation (16) of the Picard iteration into which both (30) and (31) are substituted. This is equivalent to the relation (33) and thus completes the proof of Theorem 5. □

Remark 7.

For general nonlinear (non-polynomial) state Equations (1), interval Newton methods [57,66,67] can be employed to determine the coefficients $\left[{\mathbf{x}}_{(i+1)·\nu }\right]$ and $\left[{\mathbf{R}}^{\prime }\right]$ under consideration of (2) as additional constraint. For that purpose, either a temporal series expansion on both sides of (33) in the powers $t^{i·\nu }$ is determined with subsequently extracting equalities for terms of identical orders on the left- and right-hand sides, or by determining a set of nonlinear algebraic equations after differentiating both sides with respect to the unknown coefficients.

Remark 8.

For a scalar polynomial differential equation, all interval coefficients are obtained in closed form as functions of the initial condition $\left[{\mathbf{x}}_0\right]=\left[\mathbf{x}\right]\left(0\right)$. Note, all interval coefficients $\left[{\mathbf{x}}_{(i+1)·\nu }\right]$ are then independent of the time horizon T for which the solution is computed. Only the remainder term $\left[{\mathbf{R}}^{\prime }\right]$ depends on T in this special case.

In future work, it will be investigated whether this truncated temporal series expansion can be used efficiently to determine interval enclosures for Mittag–Leffler functions with complex arguments. Such evaluations would be necessary if Theorem 3 was to be extended toward oscillatory systems in analogy to the complex-valued solution procedure that was introduced for integer-order dynamics in [48].

Remark 9.

For the case $\nu =1$, Theorem 5 represents an L-th order temporal Taylor series approximation of the solution of the considered initial value problem.

4. Numerical Case Studies

In this section, numerical examples are discussed for selected linear and nonlinear benchmark problems. Although only scalar differential equations are considered in Section 4.1, Section 4.2, Section 4.3, Section 4.4, straightforward extensions to higher-order models are possible as shown in [46] for the Mittag–Leffler type enclosure technique. Moreover, a nonlinear higher-dimensional system model is considered in Section 4.5.

To make the Mittag–Leffler enclosure technique (Theorem 3) and the series expansion approach (Theorem 5) applicable without any further restrictions, it is assumed that the fractional derivative order $\nu$ is temporally constant in all cases in which it is specified as an interval parameter. Theorem 4, however, would also be valid directly for temporally varying, but interval-bounded derivative orders due to the conservative overapproximation of the temporal integral according to Equations (28) and (29). Note that the implementation of Theorems 3 and 5 uses an equidistant interval splitting in all function evaluations to reduce the effect of the interval-related dependency effect (cf. [57]), where 5–10 subdivisions were (heuristically) chosen for each interval quantity. This is not necessary for the considered examples when using Theorem 4, where overestimation can be countered at least partially by choosing small integration step sizes.

All following simulation results are obtained by using the Matlab toolbox IntLab [68] which provides the required functionalities for interval analysis including verified implementations of standard functions and directed rounding.

4.1. Linear System with Point-Valued Parameters

As the first benchmark scenario, consider the linear fractional system model

$$x^{\left(0.5\right)}\left(t\right)=-2x\left(t\right)\mathrm{with}x\left(0\right)=x_0=1$$

(35)

with precisely known initial conditions and parameters.

The application of Theorem 3 naturally leads to the point interval $\left[\lambda \right]={\left[\lambda \right]}^{\langle \kappa \rangle }={\left[\lambda \right]}^{\langle \kappa +1\rangle }\equiv -\left[2;2\right]$ as the outcome of the iteration procedure if its simplification according to Corollary 1 is applied. Hence, the solution according to

Table 1. Comparison of all verified solution techniques for the example (35).

	Theorem 3		Theorem 4		Theorem 5
$t$	$\underset{¯}{x}\left(t\right)$	$\overline{x}\left(t\right)$	$\underset{¯}{x}\left(t\right)$	$\overline{x}\left(t\right)$	$\underset{¯}{x}\left(t\right)$	$\overline{x}\left(t\right)$
$0.1$	$0.55360625378487$	$0.55360625378489$	$0.55060278300222$	$0.55656833284754$	$0.55360625376894$	$0.55360625381019$
$0.2$	$0.45824602279222$	$0.45824602279224$	$0.45379129319789$	$0.46269008555747$	$0.45824599929468$	$0.45824605902117$
$0.3$	$0.40284721803966$	$0.40284721803968$	$0.39621103307958$	$0.40948577622286$	$0.40284553806290$	$0.40284975624289$
$0.4$	$0.36473273958222$	$0.36473273958223$	$0.35483783941352$	$0.37464047354682$	$0.36469796639549$	$0.36478445957060$
$0.5$	$0.33620400244633$	$0.33620400244635$	$0.32144602599869$	$0.35098668635366$	$0.33583913946883$	$0.33673974990786$
$0.6$	$0.31371669977514$	$0.31371669977516$	$0.29170325652201$	$0.33577066915368$	$0.31122583934236$	$0.31733441644716$
$0.7$	$0.29535370875659$	$0.29535370875660$	$0.26251634384318$	$0.32825411048854$	$0.28271413071347$	$0.31353758812897$
$0.8$	$0.27996612725009$	$0.27996612725010$	$0.23096027944650$	$0.32905324213473$	$0.22834784324700$	$0.35360334022652$
$0.9$	$0.26681423461238$	$0.26681423461239$	$0.19361581561880$	$0.34009354592870$	$0.08823371192553$	$0.51965173080075$
$1.0$	$0.25539567631050$	$0.25539567631051$	$0.14603825143031$	$0.36483039437027$	$-0.28664622848729$	$1.01757700069746$

is given by the evaluation of the interval-valued Mittag–Leffler function [64] with

$$x\left(t\right)\in \left[E_{0.5,1}\right]\left(\left[\lambda \right]·t^{0.5}\right)·x_0.$$

(36)

Table 1 further contains the guaranteed lower and upper bounds computed by Theorem 4 at the same sampling instants after the maximum number of $\kappa =100$ iterations as well as for the evaluation of Theorem 5 with the fixed order $L=20$.

In addition, the decrease of the interval width, $\mathrm{diam}\left\{\left[x\right]\left(t\right)\right\}=\overline{x}\left(t\right)-\underset{¯}{x}\left(t\right)$, at $t=0.1$ is depicted in Figure 1a in dependence of the iteration parameter $\kappa$, while Figure 1b represents the dependence of the solution width on the approximation order L at the same point of time. Here, the interval-valued Picard iteration was evaluated on an equidistant discretization mesh with $t_{i+1}-t_i={10}^{-3}$. A plateau of approximately constant solution diameters is reached already after $\kappa =10$ iterations. Hence, further work can implement step size control strategies which reduce the width of the discretization mesh to further enhance the solution quality. It can be noticed that the series-based solution representation outperforms the Picard iteration procedure with relatively small orders $L<10$ for short integration horizons, while the Picard iteration is more accurate if longer time spans are accounted for. Especially for $t>0.5$, the most accurate results are obtained by a direct evaluation of the interval extension of the real-valued Mittag–Leffler function.

However, the reasonable accuracy of both alternatives from Theorems 4 and 5 motivates a possible direction for future research: In order to make the Mittag–Leffler type enclosures applicable for systems with complex eigenvalues, the Theorems 4 and 5 could be applied as already indicated in Section 3.3 to determine guaranteed enclosures of complex-valued Mittag–Leffler functions, for which the monotonicity properties derived in [64] no longer hold.

4.2. Linear System with Interval Parameters

The general observations from the previous subsection also hold for the uncertain linear system model

$$x^{\left(\nu \right)}\left(t\right)=px\left(t\right)$$

(37)

with

$$x\left(0\right)\in \left[0.9;1.1\right],p\in \left[-2;-1\right].$$

(38)

However,

Table 2. Comparison of all verified solution techniques for the example (37) with (38) and $\nu =0.5$.

	Theorem 3		Theorem 4		Theorem 5
$t$	$\underset{¯}{x}\left(t\right)$	$\overline{x}\left(t\right)$	$\underset{¯}{x}\left(t\right)$	$\overline{x}\left(t\right)$	$\underset{¯}{x}\left(t\right)$	$\overline{x}\left(t\right)$
$0.1$	$0.49824562840638$	$0.60896687916337$	$0.19954647605972$	$0.97724861388458$	$0.44689051530214$	$0.82173567606893$
$0.2$	$0.41242142051300$	$0.50407062507146$	$-0.12396317781484$	$1.06493841015708$	$0.26212612896212$	$0.76579771031657$
$0.3$	$0.36256249623570$	$0.44313193984365$	$-0.49545940985475$	$1.31513132817254$	$0.03262962704483$	$0.86164495107816$
$0.4$	$0.32825946562400$	$0.40120601354046$	$-1.00135342456102$	$1.74043523073777$	$-0.38579044832638$	$1.20299765864330$
$0.5$	$0.30258360220170$	$0.36982440269098$	$-1.72594112115522$	$2.40556040713722$	$-1.07052397618301$	$1.82458576007731$
$0.6$	$0.28234502979763$	$0.34508836975267$	$-2.78641216611905$	$3.41954132043425$	$-2.18713576624901$	$2.89195763210021$
$0.7$	$0.26581833788093$	$0.32488907963226$	$-4.35414527639700$	$4.94951257461380$	$-4.00860384422933$	$4.67585420464389$
$0.8$	$0.25196951452508$	$0.30796273997511$	$-6.68331750998884$	$7.24716113915103$	$-6.97229005895072$	$7.61747941911901$
$0.9$	$0.24013281115114$	$0.29349565807363$	$-10.15263321750647$	$10.68960971888492$	$-11.78474276629047$	$12.43931868631812$
$1.0$	$0.22985610867945$	$0.28093524394156$	$-15.32727203998743$	$15.84097397657622$	$-19.59289086573390$	$20.33009019218864$

and Figure 2 show the effect that the series expansion approach is more accurate for shorter time horizons as compared with the Picard iteration in a more pronounced way. Due to the fact that both are computationally efficient and simple to implement, it is advantageous to intersect the results with each other (admissible due to the fact that both represent guaranteed solution enclosures) to obtain more accurate simulation routines in future work. Such intersections can also be performed with the outcome of Theorem 3, which is recommended for either not fully decoupled linear models with $n>1$ or in the nonlinear case that is addressed in the following subsections.

To show the insensitivity of the Mittag–Leffler type enclosure technique on interval parameters $\left[\nu \right]$ for the fractional derivative order, see

Table 3. Solution enclosures according to Theorem 3 for the uncertain fractional derivative order $\nu \in \left[0.5;0.6\right]$ in the example (37) with (38).

	Theorem 3
$t$	$\underset{¯}{x}\left(t\right)$	$\overline{x}\left(t\right)$
$0.1$	$0.49041110243204$	$0.67592092033214$
$0.2$	$0.40081299111450$	$0.55302527456765$
$0.3$	$0.34877260547986$	$0.47914436188981$
$0.4$	$0.31308075871594$	$0.42774986528601$
$0.5$	$0.28647584215609$	$0.38918852078824$
$0.6$	$0.26559917303058$	$0.35884192008739$
$0.7$	$0.24862936318946$	$0.33415129010471$
$0.8$	$0.23447350669331$	$0.31355890592956$
$0.9$	$0.22242767195137$	$0.29605125757402$
$1.0$	$0.21201392800064$	$0.28093524394156$

. Note, due to multiple interval dependencies on $\left[\nu \right]$, the other two solution techniques would need to be extended by techniques for reduction of interval-related overestimation (such as interval bisectioning). This is out of the scope of this paper, because the Mittag–Leffler type procedure is still more accurate than the solutions of both other alternatives, even if they are purely evaluated for a point-valued parameter $\nu$.

4.3. Polynomial System with Point-Valued Parameters

For the point-valued nonlinear system model

$$x^{\left(0.5\right)}\left(t\right)=-2x^3\left(t\right)\mathrm{with}x\left(0\right)=1,$$

(39)

the multi time step Picard iteration is by far the most accurate approach according to the comparison in

Table 4. Comparison of Theorems 4 and 5 for the example (39).

	(a) Theorem 4
$\mathit{t}$	$\underset{¯}{\mathit{x}}\left(\mathit{t}\right)$	$\overline{\mathit{x}}\left(\mathit{t}\right)$
$0.1$	$0.70197801790413$	$0.70473417747045$
$0.2$	$0.65120112737833$	$0.65656646920538$
$0.3$	$0.62007265058730$	$0.62933643671742$
$0.4$	$0.59661079870802$	$0.61148450148314$
$0.5$	$0.57671384967976$	$0.59945057578296$
$0.6$	$0.55831232761783$	$0.59179684100157$
	(a) Theorem 5
$\mathit{t}$	$\underset{¯}{\mathit{x}}\left(\mathit{t}\right)$	$\overline{\mathit{x}}\left(\mathit{t}\right)$
$0.01$	$0.82899320426943$	$0.87799864529316$
$0.02$	$-0.10790209855635$	$2.10983099954968$
$0.03$	$-7.93809912730752$	$12.68771986974560$
$0.04$	$-42.31636078898693$	$58.04678160623410$
$0.05$	$-147.886460989286$	$194.549458404715$
$0.06$	$-407.974159547404$	$525.443781220712$

as well as Figure 3 and Figure 4. Moreover, it can be noted that the Mittag–Leffler type enclosures start to worsen for longer time horizons. This was already observed in [46], where a procedure for subdividing the integration horizon in combination with a guaranteed quantification of the arising truncation errors was introduced. Especially this routine could benefit from a combination of Theorem 3 with Theorem 4, where the latter would provide the reference solutions introduced in ([46], Equation (28)). The temporal series expansion approach alone, however, turns out to be inefficient in this case and is only applicable for very short integration horizons due to the strong nonlinearity. Therefore, the series definition in its current form should only be employed if the system under consideration has a dominant linear behavior.

4.4. Polynomial System with Interval Parameters

As a final benchmark, consider the nonlinear system model

$$x^{\left(\nu \right)}\left(t\right)=p{x}^3\left(t\right)$$

(40)

with the uncertain initial condition and parameter

$$x\left(0\right)\in \left[0.9;1.1\right],p\in \left[-2;-1\right].$$

(41)

The uncertainty in this model turns the Mittag–Leffler type enclosures into the most efficient approach for this setting. This holds also if the fractional derivative order becomes uncertain, see Figure 5 and Figure 6 as well as

Table 5. Comparison of Theorems 4 and 5 for the example (40) with (41).

	(a) Theorem 4
$\mathit{t}$	$\underset{¯}{\mathit{x}}\left(\mathit{t}\right)$	$\overline{\mathit{x}}\left(\mathit{t}\right)$
$0.1$	$0.00507273865349$	$1.08477651323397$
$0.2$	$-0.41385241125397$	$1.11290925222648$
$0.3$	−	−
$0.4$	−	−
$0.5$	−	−
$0.6$	−	−
	(a) Theorem 5
$\mathit{t}$	$\underset{¯}{\mathit{x}}\left(\mathit{t}\right)$	$\overline{\mathit{x}}\left(\mathit{t}\right)$
$0.01$	$0.49954165843437$	$1.37239730772828$
$0.02$	$-10.43838386099088$	$14.31034577538966$
$0.03$	$-96.233087027468$	$115.588474698140$
$0.04$	$-458.283614635841$	$538.426532130077$
$0.05$	$-1546.31209895609$	$1794.48003417495$
$0.06$	$-4196.99132289793$	$4816.00222505225$

.

However, this result also opens up a further idea for future investigations. In the preceding sections, only the option to enhance the Mittag–Leffler type enclosures by means of Theorem 4 was discussed. However, it could be possible that one obtaines solution enclosures after a first intersection of the Picard iteration scheme with the results of Theorem 3 that can then be used in a second step for two further purposes:

• Exploit advantages of Theorem 4 in early parts of the simulation (close to $t=0$) to further improve Theorem 3;
• Use this combination of both to further implement observer approaches, where measured data are only available at some discrete points of time.

4.5. Nonlinear Higher-Dimensional System Model

As a final nonlinear application scenario, consider the fractional-order epidemiological SEIR model with a compartmental structure according to [69] which can be used to forecast the development of infectious diseases within a population. The state equations of a corresponding fractional-order model can be stated by

$$\left[\begin{array}{c}{S}^{\left(\nu \right)}\left(t\right)\\ E^{\left(\nu \right)}\left(t\right)\\ I^{\left(\nu \right)}\left(t\right)\\ N^{\left(\nu \right)}\left(t\right)\end{array}\right]=\left[\begin{array}{c}b·N\left(t\right)-p·b·E\left(t\right)-\left(q·b+r·\frac{S\left(t\right)}{N\left(t\right)}\right)·I\left(t\right)-d·S\left(t\right)\\ (p·b-\beta -d)·E\left(t\right)+\left(q·b+r·\frac{S\left(t\right)}{N\left(t\right)}\right)·I\left(t\right)\\ \beta ·E\left(t\right)-\left(\theta +d\left(N\right(t\left)\right)+\gamma \right)·I\left(t\right)\\ N\left(t\right)·(b-d)-\theta ·I\left(t\right)\end{array}\right],\left[\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\\ x_3\left(t\right)\\ x_4\left(t\right)\end{array}\right]=\left[\begin{array}{c}S\left(t\right)\\ E\left(t\right)\\ I\left(t\right)\\ N\left(t\right)\end{array}\right].$$

(42)

Here, the state variables are specified as follows:

• $S\left(t\right)$ denotes the susceptible part of a population;
• $E\left(t\right)$ the exposed group of the population;
• $I\left(t\right)$ the infectious group;
• $N\left(t\right)$ the non-constant total population size.

An important feature of this system model is the consideration that many diseases have a significant incubation period during which individuals have already been infected but are not yet infectious themselves. Those individual belong to the compartment $E\left(t\right)$ before transitioning to $I\left(t\right)$.

For the following simulations, the fractional derivative order is set to $\nu =0.95$ with the natural birth rate $b=0.001555$, the vertical transmission parameters $p=0.8$ and $q=0.95$, the horizontal transmission rate $r=0.05$, the rate $\beta =0.05$ with which an exposed individual becomes infectious, the infection-related death rate $\theta =0.002$, the recovery rate $\gamma =0.003$ and the state-dependent natural death rate

$$d\left(N\right)=0.00001+0.000007·N.$$

(43)

For further details about this model and the choice of parameters, the reader is referred to [69] and the references therein.

Due to the non-negligible nonlinearity of this system model and the partially oscillatory dynamics observed in [69], the following simulations are restricted to the novel Picard iteration approach according to Theorem 4. The simulation case study takes into account three different sets of initial conditions:

Scenario 1: 

Consideration of point-valued initial conditions in Figure 7 with

$$\left[\begin{array}{c}S\left(0\right)\\ E\left(0\right)\\ I\left(0\right)\\ N\left(0\right)\end{array}\right]\in \left[\begin{array}{c}140·\left[1.0;1.0\right]\\ 0.01·\left[1.0;1.0\right]\\ 0.02·\left[1.0;1.0\right]\\ 141·\left[1.0;1.0\right]\end{array}\right];$$

(44)

Scenario 2: 

Consideration of a $\pm 5\%$ uncertainty in the initial conditions for the compartments of exposed and infectious individuals in Figure 8 with

$$\left[\begin{array}{c}S\left(0\right)\\ E\left(0\right)\\ I\left(0\right)\\ N\left(0\right)\end{array}\right]\in \left[\begin{array}{c}140·\left[1.0;1.0\right]\\ 0.01·\left[0.95;1.05\right]\\ 0.02·\left[0.95;1.05\right]\\ 141·\left[1.0;1.0\right]\end{array}\right];$$

(45)

Scenario 3: 

Consideration of a $\pm 1\%$ uncertainty in the initial conditions for the compartments of exposed and infectious individuals in Figure 9 together with $\pm 10\%$ uncertainty in the susceptible class and population size with

$$\left[\begin{array}{c}S\left(0\right)\\ E\left(0\right)\\ I\left(0\right)\\ N\left(0\right)\end{array}\right]\in \left[\begin{array}{c}140·\left[0.90;1.10\right]\\ 0.01·\left[0.99;1.01\right]\\ 0.02·\left[0.99;1.01\right]\\ 141·\left[0.90;1.10\right]\end{array}\right].$$

(46)

In all simulations, an equidistant discretization mesh with $t_{i+1}-t_i={10}^{-1}$ has been employed.

The first scenario in Figure 7 clearly indicates that uncertainty due to the temporal discretization mesh leads to a blowup of the widths of the state enclosure after approximately 1000 integration time steps ($t=100$). The state variables that are most affected by this blowup are the exposed and infectious compartments $E\left(t\right)$ and $I\left(t\right)$, respectively. It should be pointed out that the state equations according to (42) are implemented in such a way that overestimation due to multiple interval dependencies is reduced as far as possible by factoring out common state variables. However, the wrapping effect is not countered directly by the current implementation, which directly operates on the state equations in the form (42). Therefore, future work should try to reduce this type of pessimism by using solution representations that are less prone to the wrapping effect (for example, ellipsoidal enclosures [70,71]). Alternatively, procedures for a preconditioning of the state equations can be considered which aim at computing tighter solution bounds in a transformed set of coordinates [1,4]. Note, numerical simulations have shown that a further reduction of the integration step sizes does not lead to any noticeable decrease in the widths of the computed solutions for the scenario shown in Figure 7.

According to Figure 8, an increase in the uncertainty of the initial compartment sizes of exposed and infectious individuals firstly propagates to uncertainty in the group of susceptible persons. This kind of behavior can also be observed in Figure 9, where interval uncertainty was accounted for in all four state variables. To avoid the evaluation of the system models with states that are not reasonable, negative simulation results are cut off. This is justified by the proof given in [69] that all reasonable state variables lie with certainty in the interior of the positive orthant of the state space. Especially the two cases shown in Figure 8 and Figure 9 point out a further direction for future research. Using the approach of interval observes, cf. [72] and the following conclusions section tightening of the solution enclosures will become possible on the basis of measured data. This will allow for a more efficient forecast of the propagation of diseases together with the possibility to estimate system parameters and to identify characteristics such as the death rate (43) in an online manner.

5. Conclusions and Outlook on Future Work

In this paper, it has been shown that Picard iteration procedures and temporal series expansion techniques can be employed for the development of interval-based simulation routines for fractional-order differential equations with nonlinearities and uncertain parameters. As such, they complement the Mittag–Leffler type solution technique that was already developed in previous work.

Combinations of all three proposed options will open up novel possibilities for enhancing solution enclosures not only in pure open-loop settings but also in cases where an interval observer technique such as in [72] is extended towards scenarios in which measured data are not available in continuous-time form but only at specific sampling instants. Due to the fact that the combination of various simulation techniques, as proposed in this paper, inevitably increases the required computational effort, the following two ideas could provide interesting directions for future research. Automatic recommender systems such as VERICOMP [73] could be generalized to the fractional-order case in order to automatically propose the most suited verified solver with its most appropriate settings. In addition, the CADNA software [74] (Control of Accuracy and Debugging for Numerical Applications) could be employed to determine those maximum numbers of iterations after which no solution improvements become possible because all resulting error terms purely result from numerical round-off errors.

In addition, future work should also deal with the development of estimation scheme that rely on the simulation procedures developed in this paper in order to suitably reconstruct initialization functions for fractional-order systems. As indicated, for example in [72], this is yet an open problem which is of specific practical interest if state estimation and simulation routines are not initialized at a system’s equilibrium state.