Analytical solution of fractional variable order differential equations

Accepted Manuscript Analytical solution of fractional variable order differential equations W. Malesza, M. Macias, D. Sierociuk

PII: DOI: Reference:

S0377-0427(18)30522-3 https://doi.org/10.1016/j.cam.2018.08.035 CAM 11872

To appear in:

Journal of Computational and Applied Mathematics

Received date : 6 April 2018 Revised date : 21 August 2018 Please cite this article as: W. Malesza, M. Macias, D. Sierociuk, Analytical solution of fractional variable order differential equations, Journal of Computational and Applied Mathematics (2018), https://doi.org/10.1016/j.cam.2018.08.035 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Analytical solution of fractional variable order differential equations✩ W. Maleszaa,∗, M. Maciasa , D. Sierociuka a Institute

of Control and Industrial Electronics, Warsaw University of Technology, Koszykowa 75, 00-662 Warsaw, Poland

Abstract The aim of this paper is to introduce an approach for solving fractional variable order linear differential equations. The approach is based on switching schemes that realize different types of variable order derivatives. Obtained analytical solutions are compared with numerical results. Achieved methods can be helpful, e.g., to validate existing or developed numerical algorithms for solving these types of equations. Keywords: fractional derivatives and integrals, fractional differential equations, variable order 2010 MSC: 26A33, 34A08, -

1. Introduction Fractional calculus is a generalization of traditional integer order integration and differentiation on non-integer orders. This generalization assumes that the traditional integer order calculus is only a special case of fractional order dif5

ferential calculus. The history of fractional order calculus is nearly as old as traditional differential calculus, because the first time when possibility of defining the fractional order derivative has been mentioned was in 1695 by Leibniz and L’Hˆ ospital. ✩ This work was supported by the Polish National Science Center with the decision number UMO-2014/15/B/ST7/00480. ∗ Corresponding author Email addresses: [email protected] (W. Malesza), [email protected] (M. Macias), [email protected] (D. Sierociuk)

Preprint submitted to Journal of Computational and Applied Mathematics August 28, 2018

There are many applications in scientific and engineering areas which present 10

efficiency of using fractional calculus in modelling and control. One of the examples of application is modelling of diffusive systems. For example, in heat transfer process in semi-infinite beam, the analytical solution in Laplace do√ main includes a term 1/ s which can be interpreted in time domain as a half order integrator. In [1], results of successful modeling of heat transfer process in

15

solid material were presented. Moreover, in [2], similar results for heat transfer in heterogeneous materials, described by anomalous diffusion using fractional order partial differential equation, were shown. Memory effect within a fractional calculus approach has been investigated in various areas, for example in nonlinear heat conduction [3] and viscoelastic [4, 5]. Another example of

20

successful exploitation of fractional calculus are ultracapacitors. Modeling of these devices based on anomalous diffusion (fractional order model) was presented in [6, 7, 8, 9, 10]. In [11], analysis of resonance phenomenon in circuit with ultracapacitor was investigated. This analysis clearly shows that modelling of ultracapacitors based only on integer order calculus causes unacceptable er-

25

rors in obtaining value of resonance frequency. It was also presented that the fractional order model can describe such a phenomenon, where resonance frequency is different from maximum current frequency. Not only in modeling, fractional calculus was also found interesting in signal processing. For example, in control theory, fractional calculus is used to obtain more efficient algorithms

30

like fractional order PID controllers [12] and variable order PID controllers [13]. More theoretical background and applications of fractional calculus can be found in [14, 15, 16, 17, 18, 19, 20, 21]. When the fundamental properties of a system or its structure are changing in time, a variation of the system’s order is to be observed. This behavior

35

can be described as state vector-dependent or time-dependent (or both timeand state vector-dependent). In this paper, time-dependent variable order case is taken into consideration. Such a behavior can be met, e.g. in chemistry (when the properties of the system are changing due to chemical reactions), medicine, electrochemistry, material science and other areas. In [22], variable 2

40

order differential equations were used to model a memory behavior of shapememory polymer (SMP). Presented results clearly showed that the variable order model is more suitable than constant order. In [23] usage of variable order fractional calculus to describe the evolution of protein lateral diffusion ability in cell membranes is presented. In [24], process of an electrochemical

45

device with parameters changing in time was presented. Authors have shown that the order of used model depends on time, which implies that for describing long-time behavior, a variable order model has to be used. In [25], the variable order equations are used to describe the history of drag expression. Results of very accurate modelling of heat transfer process in a time-varying grid-hole

50

structure are presented in [26]. The description of the variable order systems is more complex than in the constant order case. In [27, 28], three general types of variable order derivative definitions can be found, however, these definitions were given without interpretation and derivation. Some other modifications of the mentioned definitions

55

are presented in [29]. Interpretations of some of those definitions in the form of switching schemes were presented in [30, 31, 32]. Based on those switching schemes it is possible to categorize fractional order derivatives not only according to equation’s aspects but also according to its behavior and practical properties. Moreover, based on proposed switching schemes, it is possible to

60

build analog models of variable order systems. Recursive approximations for variable order fractional operators with applications are presented in [33]. Experimental validation of variable order analog models was presented in [34, 35]. The duality between variable order definitions, that is a very important tool for analysis, was introduced in [36]. The existence of solution form of variable order

65

differential equations is discussed in [37, 38]. In [39], a fractional variable order derivative model for chloride ions sub-diffusion, was used. To the best of our knowledge, there are no general methods for finding an analytical solution of linear fractional variable order equations (in contradiction to constant fractional order case [14, 17, 40]). In [41], preliminary results for

70

finding an analytical solution of a particular case of variable order derivative 3

definition, are presented. It is also an explanation, why numerical and analog modelling methods, presented e.g. in [41, 42, 43, 44, 45], are so important. In this paper, a method for solving fractional variable order linear differential equations for different types of variable order derivative definitions (e.g. A75

type, B-type, D-type, E-type) are presented. Obtained methods are validated by comparison with numerical results. The paper is organized as follows. In Section 2 variable order derivative definitions are recalled and explained. In Section 3, main results of this paper— the methods for finding the analytical solutions of linear fractional variable order

80

equations (for D-type, A-type, B-type and E-type definitions, respectively) and their numerical validations are presented. Finally, in Section 4 we provide some concluding remarks and an outlook on future research. 2. Variable order derivatives As a base of generalization of the constant fractional order derivative on

85

variable order case, the following Gr¨ unwald-Letnikov (G-L) definition is taken into consideration: Definition 1. Fractional constant order derivative is defined as follows:   n 1 X α r α f (t − rh), (−1) 0 Dt f (t) = lim α h→0 h r r=0

(1)

where α ∈ R, h > 0 is a step time, and n = ⌊t/h⌋.

For the case of order changing in time (variable order case), the following types of definitions can be found in the literature [27, 28, 30, 31, 32]. 90

The first one is obtained by replacing in (1) a constant order α by variable order α(t). In that approach, all of the coefficients for past samples are obtained for present values of the order. Definition 2 ([31]). The A-type of fractional variable order derivative is defined as follows: A α(t) f (t) 0 Dt

= lim

h→0

1 hα(t)

n X r=0

4

(−1)r



 α(t) f (t − rh). r

The second type of definition assumes that coefficients for past samples are obtained for order that was present for these samples. In this case, the definition 95

has the following form: Definition 3 ([30]). The B-type of fractional variable order derivative is defined as follows: B α(t) f (t) 0 Dt

  n X (−1)r α(t − rh) = lim f (t − rh). h→0 r hα(t−rh) r=0

The third and fourth type of variable order derivatives are defined recursively and have the following forms: Definition 4 ([31]). The D-type fractional variable order derivative is defined as follows: D α(t) f (t) 0 Dt



n X



−α(t) f (t) (−1)j = lim  α(t) − h→0 j h j=1





D α(t)  0 Dt−jh f (t) .

Definition 5 ([32]). The E-type fractional variable order derivative is defined as follows: E α(t) f (t) 0 Dt



n X





α(t−jh)

f (t) −α(t − jh) h = lim  α(t) − (−1)j h→0 j h hα(t) j=1



E α(t)  0 Dt−jh f (t) .

Remark 1 ([31]). For a constant order α(t) = const, the same results as for constant order derivative definitions are obtained, i.e., A α(t) f (t) 0 Dt

α D α E α α =B 0 Dt f (t) = 0 Dt f (t) = 0 Dt f (t) = 0 Dt f (t).

In further considerations, we will assume that the fractional variable order α(t) is a piecewise constant function of time, i.e., α(t) = αi for ti−1 ≤ t < ti , 100

i = 1, . . . , N , (usually we take t0 = 0), where αi is a real constant, and N is a natural finite number. In contrast to the variable order derivatives, the Laplace transform of constant order G-L derivative of function f (t) (with zero initial conditions) can be obtained, which is given as follows [46] α L{0 Dα t f (t)} = s F (s),

5

(2)

where L is a symbol of Laplace transform (inverse Laplace transform will be denoted by L−1 ).

2.1. Switching schemes of variable order derivatives 105

2.1.1. Switching scheme of A-type derivative Let us consider the following so-called output-reductive switching order scheme in the serial form (serial o-r scheme), denoted o-r ser Ξ{α(t)} and presented in Fig. 1, where αN -block and α ˆ i -blocks, i = 1, . . . , N − 1, denote fractional transfer func-

ˆ i (s) = sαˆ i , respectively. Moreover, tions GN (s) = sαN and G i = 1, . . . , N − 1.

αi = αi+1 + α ˆi,

(3)

The switches Si , i = 1, . . . , N , take the following positions   b for ti−1 ≤ t < ti , i = 1, . . . , N. Si =  a otherwise, o-r ser Ξ{α(t)}

f (t)

o-r ser y(t) b

αN

b

a

b

a

SN

a

α ˆ2

S3

b a

α ˆ1

S2

S1

Figure 1: Structure of output-reductive switching order scheme in serial form

o-r Ξ{α(t)} ser

(presented configuration: switching from α1 to α2 ).

Remark 2. As it was shown (in [31]) the A-type derivative can be realized by means of output-reductive switching order scheme. Next, let us consider the so-called output-reductive switching order scheme in the parallel form (parallel o-r scheme) o-r par Ξ{α(t)} depicted in Fig. 2, where, analogously as in the case of serial o-r scheme (see Fig. 1), αi -blocks, i = 1, . . . , N , denote fractional transfer functions Gi (s) = sαi . The switch takes the positions S=i

for

ti−1 ≤ t < ti , 6

i = 1, . . . , N.

o-r par y(t)

1

y2

α2

f (t)

o-r par Ξ{α(t)}

y1

α1

2

S

N

yN

αN

Figure 2: Structure of output-reductive switching order scheme in parallel form

o-r Ξ{α(t)} par

(presented configuration: switching from α1 to α2 ).

Lemma 1. From the input-output point of view the serial and parallel o-r schemes are equivalent. Proof. First, we calculate the output signal (solution) of the serial o-r scheme. Thus, during the time interval ti−1 ≤ t < ti , i = 1, . . . , N , we have (in s-domain)   N −i Y o-r ˆ N −j (s) F (s) GN (s) G ser Yi (s) = j=1

= Gi (s)F (s),

(4)

since (and by (3)) s αN

N −i Y

sαˆ N −j = sαN +

j=1

PN −i j=1

ˆ N −j sα

= s αi .

Therefore, in time domain, the solution of the serial o-r scheme is o-r ser y(t)

= L−1 {o-r ser Yi (s)}

for

ti−1 ≤ t < ti , i = 1, . . . , N.

(5)

Next, the solution of parallel o-r scheme is the following, that is, in s-domain, during the time interval ti−1 ≤ t < ti , i = 1, . . . , N , we have o-r par Yi (s)

= Gi (s)F (s),

(6)

and thus, in time domain, o-r par y(t)

= L−1

o-r

par Yi (s)



for 7

ti−1 ≤ t < ti , i = 1, . . . , N.

(7)

Thus, the serial and parallel o-r scheme are equivalent, because the solutions given by (5) and (7) are the same (for the same f (t) and αi , i = 1, . . . , N ,), i.e., o-r ser y(t)

110

≡ o-r par y(t).

(8)

Remark 3. From Lemma 1 and Remark 2 it follows that A-type derivative can be realized by means of both serial- and parallel o-r scheme. 2.1.2. Switching scheme of D-type derivative Let us consider the following so-called input-reductive switching order scheme in the serial form (serial i-r scheme), denoted i-r ser Ξ{α(t)} and presented in Fig. 3, where αN -block and α ˆ i -blocks, i = 1, . . . , N − 1, denote fractional transfer func-

ˆ i (s) = sαˆ i , respectively. Moreover, tions GN (s) = sαN and G i = 1, . . . , N − 1.

αi = αi+1 + α ˆi,

(9)

The switches Si , i = 1, . . . , N , take the following positions   b for ti−1 ≤ t < ti , Si = i = 1, . . . , N.  a otherwise, i-r ser Ξ{α(t)}

f (t) b a

b

α ˆ1 S1

a

b

b a

S2

SN −1

α ˆ N −1

a

αN

i-r ser y(t)

SN

Figure 3: Structure of input-reductive switching order scheme in serial form

i-r Ξ{α(t)} ser

(pre-

sented configuration: switching from α1 to α2 ).

Remark 4. As it was shown (in [31]) the D-type derivative can be realized by means of input-reductive switching order scheme. 115

2.1.3. Switching scheme of B-type derivative Let us consider the following so-called input-additive switching order scheme in the serial form (serial i-a scheme), denoted i-a ser Ξ{α(t)} and presented in Fig. 4, 8

where α1 -block and α ¯ i -blocks, i = 1, . . . , N , denote fractional transfer functions ¯ i (s) = sα¯ i , respectively. Moreover, G1 (s) = sα1 and G αi = αi−1 + α ¯i,

i = 2, . . . , N.

(10)

The switches Si , i = 1, . . . , N , take the following positions   b for ti−1 ≤ t < ti , i = 1, . . . , N. Si =  a otherwise, i-a ser Ξ{α(t)}

f (t) b a

b

b

α ¯N SN

a

a

SN −1

b

α ¯2 S2

a

α1

i-a ser y(t)

S1

Figure 4: Structure of input-additive switching order scheme in serial form

i-a Ξ{α(t)} ser

(pre-

sented configuration: switching from α1 to α2 ).

Remark 5. As it was shown (in [30]) the B-type derivative can be realized by means of input-additive switching order scheme. Next, consider the so-called input-additive switching order scheme in the parallel form (parallel i-a scheme)

i-a par Ξ{α(t)}

depicted in Fig. 5, where, analogously as

in the case of serial i-a scheme (see Fig. 4), αi -blocks, i = 1, . . . , N , denote fractional transfer functions Gi (s) = sαi . The switch S takes the positions S=i

for

ti−1 ≤ t < ti ,

9

i = 1, . . . , N.

1

f (t)

S

y2

α2

2

i-a par Ξ{α(t)}

y1

α1

+

i-a par y(t)

N

yN

αN

Figure 5: Structure of input-additive switching order scheme in parallel form

i-a Ξ{α(t)} par

(presented configuration: switching from α1 to α2 ).

Lemma 2. From the input-output point of view the serial and parallel i-a schemes are equivalent. Proof. First, we calculate the output signal (solution) of the serial i-a scheme (see [47]). It is helpful to define the rectangular unit function wi (t) = H(t − ti−1 ) − H(t − ti ), 120

where H(t) is a Heaviside step function. Thus, during ith time interval ti−1 ≤ t < ti , i = 1, . . . , N , the solution (in s-domain) is i-a ser Yi (s)



= G1 (s) L{u(t)w1 (t)} + = G1 (s)L{u(t)w1 (t)} +

i X j=2

i X

j Y

k=2

!



ˆ k (s) L{u(t)wj (t)} G

Gj (s)L{u(t)wj (t)},

(11)

j=2

since (and by (10)) s

α1

i Y

sα¯ j = sα1 +

j=2

Pi

j=2

¯j sα

= s αi .

Therefore, in time domain, the solution of the serial i-a scheme is i-a ser y(t)

= L−1

i-a

ser Yi (s)



for

10

ti−1 ≤ t < ti , i = 1, . . . , N.

(12)

Next, the solution of parallel i-a scheme is the following, that is, during the time interval ti−1 ≤ t < ti , i = 1, . . . , N , we have (in s-domain) i-a par Yi (s)

= Yi−1 (s) + Gi (s)L {u(t)wi (t)} ,

(13)

where Y1 (s) = G1 L {u(t)w1 (t)}, and thus, in time domain, i-a par y(t)

= L−1

i-a

par Yi (s)



for

ti−1 ≤ t < ti , i = 1, . . . , N.

(14)

Thus, the serial and parallel i-a scheme are equivalent, because solutions (12) and (14) are the same (for fixed f (t) and αi , i = 1, . . . , N ,), i.e., ≡ i-a par y(t).

i-a ser y(t)

(15)

Remark 6. From Lemma 2 and Remark 5 it follows that B-type derivative can be realized by means of both serial- and parallel i-a scheme. 2.1.4. Switching scheme of E-type derivative Let us consider the following so-called output-additive switching order scheme in the serial form (serial o-a scheme), denoted o-a ser Ξ{α(t)} and presented in Fig. 6, where α1 -block and α ¯ i -blocks, i = 2, . . . , N , denote fractional transfer functions ¯ i (s) = sα¯ i , respectively. Moreover, G1 (s) = sα1 and G i = 1, . . . , N − 1.

αi+1 = αi + α ¯ i+1 ,

(16)

The switches Si , i = 1, . . . , N , take the following positions   b for ti−1 ≤ t < ti , i = 1, . . . , N. Si =  a otherwise, o-a ser Ξ{α(t)}

f (t)

o-a ser y(t) b

α1

a S1

α ¯2

b

b

b a

a

α ¯N

SN −1

S2

Figure 6: Structure of output-additive switching order scheme in serial form sented configuration: switching from α1 to α2 ).

11

a SN

o-a Ξ{α(t)} ser

(pre-

Remark 7. As it was shown (in [32]) the E-type derivative can be realized by 125

means of output-additive switching order scheme.

3. Solution of linear switching order equations—switching schemes approach 3.1. Solving variable (switching) order differential equation Let us define two 4-tuples T = {D, E, A, B} and T˜ = {A, B, D, E}, where

130

Tℓ and T˜ℓ denote the ℓ-th elements of T and T˜ , respectively. Analogously,

we define the 4-tuple X = {o-r, i-a, i-r, o-a} and 2-tuple F = {ser, par}. It

is already known [31, 36] that the T˜ℓ -type derivative is equivalent to (can be realized by means of) a corresponding switching scheme, which will be denoted Xℓ Fj Ξ{α(t)}, 135

for ℓ ∈ {1, . . . , 4}, j ∈ {1, 2}. The elements of X denote types of

switching scheme: output-reductive (o-r), input-additive (i-a), input-reductive (i-r), output-additive (o-a); and the elements of F denote form of switching scheme which can be serial or parallel (since for D-type and E-type derivatives parallel forms of corresponding switching systems are not known, the index j equals 1 for ℓ ∈ {3, 4}).

140

In this paper, the switching schemes realizing variable (switching) order derivatives will be used to solve analytically fractional variable order differential equations. Consider the following fractional Tℓ -type, ℓ ∈ {1, . . . , 4}, variable order differential equation (shortly Tℓ -differential equation) Tℓ α(t) x 0 Dt

= λ(t) (a(t)γ(t)x + u) ,

x(0) = 0,

ℓ ∈ {1, . . . , 4},

y(t) = γ(t)x(t),

(17a) (17b)

where x = x(t) is the real valued unknown function, u = u(t) is a real valued known function, and the following piecewise constant functions, all defined on ti ≤ t < ti+1 , for i = 0, . . . , N , where t0 = 0: α(t) = αi ∈ R+ ,

λ(t) = λi ∈ R,

a(t) = ai ∈ R, 12

γ(t) = γi ∈ R.

(18)

Using the duality property [31, 36], i.e.,

T˜ℓ −α(t) 0 Dt

to (17), we get ˜

−α(t)

y(t) = γ(t)T0 ℓ Dt

(λ(t)(a(t)γ(t)x + u)) ,



Tℓ α(t) x 0 Dt



= x, applied

ℓ ∈ {1, . . . , 4}.

(19)

Realization of (17), and thereby of (19), is depicted in Fig. 7, where the integration block

T˜ℓ −α(t) 0 Dt

u(t)

is replaced by switching scheme

Xℓ Fj Ξ{−α(t)}.

Tℓ α(t) x Xℓ 0 Dt Ξ{−α(t)}

y(t)

+

λ(t)

Fj

x

γ(t)

a(t)

Figure 7: Realization of Tℓ -type differential equation (17).

System (17) can also be rewritten as Tℓ α(t) x 0 Dt

= a(t)λ(t)γ(t)x + u ˜(t),

x(0) = 0,

ℓ ∈ {1, . . . , 4},

y(t) = γ(t)x(t),

(20a) (20b)

where u ˜(t) = λ(t)u(t),

(20c)

and then (19) modified to ˜

−α(t)

y(t) = γ(t)T0 ℓ Dt 145

(λ(t)a(t)γ(t)x + u ˜) ,

ℓ ∈ {1, . . . , 4}.

(20d)

Thus, the scheme in Fig. 7 can be obviously presented, what is evident, in an equivalent alternative form depicted in Fig. 8. u(t)

λ(t)

u ˜(t)

+

Tℓ α(t) x Xℓ 0 Dt Ξ{−α(t)} Fj

λ(t)

a(t)

x(t)

γ(t)

y(t)

γ(t)

Figure 8: Realization of Tℓ -type differential equation (17) in an alternative form.

13

Finding a solution of system depicted in Fig. 7 (or Fig. 8) will depend, and be based, on particular type of switching scheme

Xℓ Fj Ξ{−α(t)},

what will be

meticulously described in subsequent sections. 150

3.1.1. Solution of D-type differential equation Let us consider the Tℓ -type differential equation (17) for ℓ = 1, i.e., D α(t) x 0 Dt

= λ(t) (a(t)γ(t)x + u) ,

x(0) = 0,

y(t) = γ(t)x(t).

(21a) (21b)

By (19), it can be rewritten as −α(t)

y(t) = γ(t)A 0 Dt

(λ(t)(a(t)γ(t)x + u)) .

(22)

We will use the following objects (for i = 1, . . . , N ): Gi (s) =

γi , s αi − a i λ i γ i

˘ i (s) = γi G s αi

(23)

and wi (t) = H(t − ti−1 ) − H(t − ti ),

w ¯i (t) = H(t) − H(t − ti ),

where H(t) denotes Heaviside step function, and for i = 1 we have w1 = w ¯1 . Therefore, we have the following result. Proposition 1. Solution of D-type differential equation (21) is the following y(t) =

N X

yi (t)wi (t),

(24a)

i=1

where yi (t) = L−1 {Gi (s)L{ui (t)}},

i = 1, . . . , N,

ui (t) = ui−1 (t)w ¯i−1 (t) + u ˜(t)wi (t),

u1 = u ˜(t)wi (t),

(24b) (24c)

ui−1 (t) = ei−1 (t) − ai λi y˘i (t)

(24d)

ei−1 (t) = ui−1 (t) + ai−1 λi−1 yi−1 (t)

(24e)

˘ i (s)L{ei−1 (t)}}, y˘i (t) = L−1 {G

(24f)

where u ˜(t) is given by (20c). 14

Proof. In (22), the 155

schemes

A −α(t) 0 Dt

A Fj Ξ{−α(t)},

integral can be realized by means of switching

j = 1, 2. Thus, finding the solution of (21) gives rise

to the problem of finding the solution of the system presented in Fig. 8 with integration action

A −α(t) 0 Dt

by parallel o-r scheme

realized, without loss of generality (by Lemma 1),

o-r par Ξ{−α(t)}

(for ℓ = 1 and j = 2, i.e., X1 = o-r and

F2 = par) depicted in Fig. 2. Therefore, in time interval 0 ≤ t < t1 , we have the following: Y1 (s) = G1 (s)U1 (s),

U1 (s) = L{˜ u(t)w1 (t)}

˘ 2 (s)E1 (s), Y˘2 (s) = G

E1 (s) = U1 (s) + a1 λ1 Y1 (s)

Y (s) = Y1 (s). During time interval 0 ≤ t < t2 , regarding α2 -block, we artificially assume that we are concerned with G2 (s) (closed system of order α2 ). Knowing that during time 0 ≤ t < t1 the α2 -block was fed with signal E1 (s), we can determine dummy input providing this, i.e., U 1 (s) = E1 (s) − a2 λ2 Y˘2 (s). During the whole time interval 0 ≤ t < t2 , the input signal is u2 (t) = L−1 {U 1 (s)}w ¯1 (t) + u ˜(t)w2 (t), and then Y2 (s) = G2 (s)U2 (s) ˘ 3 (s)E2 (s), Y˘3 (s) = G

E2 (s) = U2 (s) + a2 λ2 Y2 (s)

y(t) = y1 (t)w1 (t) + y2 (t)w2 (t). During time interval 0 ≤ t < t3 , regarding α3 -block, we artificially assume that we are concerned with G3 (s) (closed system of order α3 ). Knowing that during time 0 ≤ t < t2 the α3 -block was fed with signal E2 (s), we can determine dummy input providing this, i.e., U 2 (s) = E2 (s) − a3 λ3 Y˘3 (s). 15

During the whole time interval 0 ≤ t < t3 , the input signal is u3 (t) = L−1 {U 2 (s)}w ¯2 (t) + u ˜(t)w3 (t), and then Y3 (s) = G3 (s)U3 (s) y(t) = y1 (t)w1 (t) + y2 (t)w2 (t) + y3 (t)w3 (t). 160

For the consecutive switches in time, we proceed analogously, obtaining, in general, iterated algorithm (24). Example 1. Using (24), we calculate the solution of D-differential equation (21), where u(t) = H(t), λ(t) = 1, and          2 for t ∈ [0, 1), 1 for t ∈ [0, 1), −1 for t ∈ [0, 1),          α(t) = 1 for t ∈ [1, 2), γ(t) = 2 for t ∈ [1, 2), a(t) = −2 for t ∈ [1, 2),               −3 for t ∈ [2, 10), 2 for t ∈ [2, 10), 3 for t ∈ [2, 10), and then G1 (s) =

1 , s2 + 1

G2 (s) =

2 , s+4

˘ 2 (s) = 2 , G s

G3 (s) =

The solution of (21) is     1 − cos(t),    
 1   + e4−4t 2 sin (1) − 12 ,  2  
y(t) = 1 + 1 cos (3t − 7) e−4 − 3 + cos (3t − 6) 3 2   
   + 38 1 − e−4 sin (3t − 5) − sin (3t − 7) +     
 + 1 cos (3t − 5) 3 + e−4 2

3 , s2 + 9

˘ 3 (s) = 3 . G s2

for t ∈ [0, 1), 155 48 2 3

+

3 −4 16 e







for t ∈ [1, 2), +

sin (3t − 6) − for t ∈ [2, 10).

Plots of analytical and numerical solutions of differential equation (21) are

depicted in Fig. 9, and the difference plot between these solutions is presented in Fig. 10, where the numerical step size is 0.002.

16

Figure 9: Solution plots (solid line—analytical, dotted line—numerical (for step size 0.002)) of differential equation (21) from Ex. 1.

Figure 10: Plot of difference between analytical and numerical solutions from Ex. 1.

Taking a(t) = 0 in (21), we obtain D α(t) x 0 Dt

= λ(t)u,

x(0) = 0,

y(t) = γ(t)x(t),

(25a) (25b)

which yields an integral equation y(t) = γ(t)A D−α(t) (λ(t)u) .

(26)

165

Corollary 1. The solution of (26), and thereby of (25), is given by (24a),

17

where yi = L−1 {Gi L{ui }},

i = 1, . . . , N,

(27a)

ui = ui−1 w ¯i−1 + u ˜(t)wi

(27b)

=u ˜(t), and u ˜(t) is given by (20c). Thus, (27) yields the same as (7), and thereby, thanks to (8), as (5). Example 2. Consider integral system (26), where λ(t) = 1, γ(t) = 1, and

α(t) =

   0.5       1   6       1.5       1

   0.25        2        0.5      1 6

for t ∈ [0, 1), for t ∈ [1, 2), for t ∈ [2, 3), for t ∈ [3, 4), for t ∈ [4, 5), for t ∈ [5, 6), for t ∈ [6, 7), for t ∈ [7, 8),

u(t) =

   1        1        0       0    1        0        1      1

for t ∈ [0, 1), for t ∈ [1, 2), for t ∈ [2, 3), for t ∈ [3, 4), for t ∈ [4, 5), for t ∈ [5, 6), for t ∈ [6, 7), for t ∈ [7, 8).

The solution is  √ 2 t    √π ,   √   3Γ(5/6) 6 t   ,  π    √   4t3/2 +(2−t) −2+t  √ ,  3 π   √ √ √ √  4 4 4    2 2Γ(3/4)( t− −2+t+ −4+t) , π y(t) =    − 13  2 + 3t,   √ √ √ √ √   2 t− −2+t+ −4+t− −5+t+ −6+t  √  ,  π   √ √ √ √  3Γ(5/6)( √ 6 6 6 6  t− −2+t+ −4+t− −5+t+ 6 −6+t)   ,  π   √ √ √ √ √  6 6 6 6 6 6  3Γ(5/6)( t− −2+t+ −4+t− −5+t+ −6+t− √ −8+t)  π

18

for t ∈ [0, 1), for t ∈ [1, 2), for t ∈ [2, 3), for t ∈ [3, 4), for t ∈ [4, 5), for t ∈ [5, 6), for t ∈ [6, 7), ,

for t ∈ [7, 8).

Plots of analytical and numerical solutions of differential equation (26) are depicted in Fig. 11, and the difference plot between these solutions is presented in 170

Fig. 10, where the numerical step size is 0.002.

Figure 11: Solution plots (solid line—analytical, dotted line—numerical (for step size 0.002)) of integral equation (26) from Ex. 2.

Figure 12: Plot of difference between analytical and numerical solutions from Ex. 2.

3.1.2. Solution of A-type differential equation Let us consider the Tℓ -type differential equation (17) for ℓ = 3, i.e., A α(t) x 0 Dt

= λ(t) (a(t)γ(t)x + u) ,

x(0) = 0,

y(t) = γ(t)x(t).

(28a) (28b)

By (19), it can be rewritten as −α(t)

y(t) = γ(t)D 0 Dt

(λ(t)(a(t)γ(t)x + u)) . 19

(29)

We will use the following objects (for i = 1, . . . , N , j = 1, . . . , N − 1): Gi (s) =

ˆ j (s) = 1 G sαˆ j

1 , s αi − a i λ i γ i

(30)

and wi = H(t − ti−1 ) − H(t − ti ),

w ¯i = H(t) − H(t − ti ),

where H(t) denotes Heaviside step function, and for i = 1 we have w1 = w ¯1 . Therefore, we have the following result. Proposition 2. Solution of A-type differential equation (28) is the following y(t) =

N X

γ i yi wi ,

(31a)

i=1

where yi (t) = L−1 {Gi L{ui (t)}},

i = 1, . . . , N,

ui (t) = ui−1 (t)w ¯i−1 (t) + u ˜(t)wi (t),

u1 = u ˜(t)wi (t),

(31b) (31c)

ui−1 (t) = yˆi−1 (t) − ai λi γi yi−1 (t)

(31d)

ˆ i−1 L{ei−1 (t)w yˆi−1 (t) = L−1 {G ¯i−1 (t)}}

(31e)

ei−1 (t) = ui−1 (t) + ai−1 λi−1 γi−1 yi−1 (t),

(31f)

where u ˜(t) is given by (20c). 175

Proof. Finding the solution of (28) gives rise to the problem of finding the solution of the system presented in Fig. 7 with integration action alized by serial i-r scheme

i-r par Ξ{−α(t)}

D −α(t) 0 Dt

re-

(for ℓ = 3 and j = 1, i.e., X3 = i-r and

F1 = ser) depicted in Fig. 3. In this proof, we will use an obvious alternative form of Fig. 7 with the γ(t)-block placed duplicated after the output-node 180

instead of single one placed ahead of it. Therefore, in time interval 0 ≤ t < t1 , we have the following: Y1 (s) = G1 (s)U1 (s),

U1 (s) = L{˜ u(t)w1 (t)}

Y (s) = γ1 Y1 (s) y(t) = γ1 y1 (t). 20

Moreover, ˆ 1 (s)Ew1 (s), Yˆ1 (s) = G

Ew1 (s) = L {(u(t) + a1 λ1 γ1 y1 (t))w ¯1 (t)} .

During time interval 0 ≤ t < t2 , we artificially assume (considering 0 ≤ t < t1 ) that we are concerned with G2 (s) (closed system of order α2 with already rejected α ˆ 1 -block). Knowing that during time 0 ≤ t < t1 the α ˆ 1 -block was

fed with signal Ew1 (s), and thereby the α ˆ 2 -block was fed with signal Yˆ1 (s) = ˆ 1 (s)Ew1 (s), we can determine dummy input providing this, i.e., G U 1 (s) = Yˆ1 (s) − a2 λ2 γ2 Y1 (s). Thus, during the whole time interval 0 ≤ t < t2 , the input signal is u2 (t) = L−1 {U 1 (s)}w ¯1 (t) + u ˜(t)w2 (t), and then Y2 (s) = G2 (s)U2 (s) y(t) = γ1 y1 (t)w1 (t) + γ2 y2 (t)w2 (t). During time interval 0 ≤ t < t3 , we artificially assume (considering 0 ≤ t < t2 ) that we are concerned with G3 (s) (closed system of order α3 , with already rejected α ˆ 1 -block and α ˆ 2 -block). Knowing, from the previous considerations, that during time 0 ≤ t < t2 the α ˆ 2 -block was fed with signal U2 (s), and ˆ 2 (s)Ew2 (s), where Ew2 (s) = thereby the α ˆ 3 -block was fed with signal Yˆ2 (s) = G

L {(u2 + a2 λ2 γ2 y2 (t))w ¯2 (t)}, we can determine dummy input providing this, i.e., U 2 (s) = Yˆ2 (s) − a3 λ3 γ3 Y2 (s). Thus, during the whole time interval 0 ≤ t < t3 , the input signal is u3 (t) = L−1 {U 2 (s)}w ¯2 (t) + u ˜(t)w3 (t), and then Y3 (s) = G3 (s)U3 (s) y(t) = γ1 y1 (t)w1 (t) + γ1 y2 (t)w2 (t) + γ1 y3 (t)w3 (t). 21

For the consecutive switchings in time, we proceed analogously, obtaining, in general, iterated algorithm (31). Example 3. Using (31), we calculate the solution of A-differential equation (28), where u(t) = H(t), λ(t) = 1, and          2 for t ∈ [0, 1), 1 for t ∈ [0, 1), −1 for t ∈ [0, 1),          α(t) = 1 for t ∈ [1, 2), γ(t) = 2 for t ∈ [1, 2), a(t) = −2 for t ∈ [1, 2),               −3 for t ∈ [2, 10), 2 for t ∈ [2, 10), 3 for t ∈ [2, 10), and then G1 (s) =

1 , s2 + 1

G2 (s) =

1 , s+4

G3 (s) =

1 , s2 + 9

ˆ 1 (s) = 1 , G s

The solution of (28) is     1 − cos(t),    
1   4−4t 3  e 2 − 2 cos (1) + 2 ,  
y(t) = 1 + 2e−4 sin (−5 + 3t) − 3 sin (−6 + 3t) + sin (−7 + 3t) 3 2     − 3 e−4 (cos (−5 + 3t) + cos (−7 + 3t))     4  
 + cos (−6 + 3t) 9 e−4 + 5 4 12

ˆ 2 (s) = s. G

for t ∈ [0, 1), for t ∈ [1, 2),

for t ∈ [2, 10).

Plots of analytical and numerical solutions of differential equation (28) are

depicted in Fig. 13, and the difference plot between these solutions is presented 185

in Fig. 14, where the numerical step size is 0.002.

22

Figure 13: Solution plots (solid line—analytical, dotted line—numerical (for step size 0.002)) of differential equation (28) from Ex. 3.

Figure 14: Plot of difference between analytical and numerical solutions from Ex. 3.

Taking a(t) = 0 in (28), we get A α(t) x 0 Dt

= λ(t)u,

x(0) = 0,

y(t) = γ(t)x(t),

(32a) (32b)

which yields an integral equation y(t) = γ(t)D D−α(t) (λ(t)u) .

23

(33)

Corollary 2. The solution of (33), and thereby of (32), is given by (31a), where yi (t) = L−1 {Gi (s)L{ui (t)}} , i = 1, . . . , N, n o ˆ i−1 (s)L{ui−1 (t)} w ui (t) = L−1 G ¯i−1 (t) + u ˜(t)wi (t),

(34a) u1 = u ˜(t)w1 (t), (34b)

where Gi ’s are taken with ai = 0, and u ˜(t) is given by (20c). Example 4. Consider integral system (33), where λ(t) = 1, γ(t) = 1, and      1 for t ∈ [0, 1), 1 for t ∈ [0, 1), 2 α(t) = u(t) =   1 for t ∈ [1, 2), 1 for t ∈ [1, 2). The solution is

y(t) =

 √   2√ t , π

  √2 + t − 1, π

for t ∈ [0, 1), for t ∈ [1, 2).

Plots of analytical and numerical solutions of differential equation (33) are depicted in Fig. 15, and the difference plot between these solutions is presented 190

in Fig. 16, where the numerical step size is 0.002.

Figure 15: Solution plots (solid line—analytical, dotted line—numerical (for step size 0.002)) of integral equation (33) from Ex. 4.

24

Figure 16: Plot of difference between analytical and numerical solutions from Ex. 4.

3.1.3. Solution of E-type differential equation Let us consider the Tℓ -type differential equation (17) for ℓ = 2, i.e., E α(t) x 0 Dt

= λ(t) (a(t)γ(t)x + u) ,

x(0) = 0,

y(t) = γ(t)x(t).

(35a) (35b)

By (19), it can be rewritten as −α(t)

y(t) = γ(t)B 0 Dt

(λ(t)(a(t)γ(t)x + u)) .

(36)

We will use the following objects (for i = 1, . . . , N ): Gi (s) =

1 , s αi − a i λ i γ i

˘ i (s) = 1 G s αi

and wi (t) = H(t − ti−1 ) − H(t − ti ), Therefore, we have the following result.

wγ (t) =

N X

(37)

γi wi (t).

i=1

Proposition 3. Solution of E-type differential equation (35), on time interval 0 ≤ t ≤ ti , i = 1, . . . , N , is the following   i−1 X y˘j (t) wγ (t), y(t) = yi (t) + j=1

25

(38a)

where yi (t) = L−1 {Gi (s)L{ui (t)}}   i−1 X ˜(t) + ai λi γi ui (t) = u y˘j (t) wi (t),

(38b) u1 (t) = u ˜(t)w1 (t)

(38c)

j=1

n o ˘ j (s)L{ej (t)} , y˘j (t) = L−1 G ej (t) =

u ˜(t) + aj λj γj

j = 1, . . . , i − 1 !! j−1 X yj (t) + y˘k (t) wj (t),

(38d) (38e)

k=1

where u ˜(t) is given by (20c). Proof. Finding the solution of (35) gives rise to the problem of finding the so195

−α(t)

lution of the system presented in Fig. 8 with integration action B 0 Dt without loss of generality (by Lemma 2), by parallel i-a scheme

realized,

i-a par Ξ{−α(t)}

(for

ℓ = 2 and j = 2, i.e., X2 = i-a and F2 = par) depicted in Fig. 4. Therefore, in time interval 0 ≤ t < t1 , we have the following: Y1 (s) = G1 (s)U1 (s) Y (s) = γ1 Y1 (s),

y(t) = γ1 y1 (t).

Moreover, since during time 0 ≤ t < t1 the α1 -block was fed with error signal coming from the sumator on the input, and after this time input signal equals 0, i.e., e1 (t) = (˜ u(t) + a1 λ1 γ1 y1 (t)) w1 (t), we have ˘ 1 (s)E1 (s). Y˘1 (s) = G During time interval 0 ≤ t < t2 , the input to and output from α2 -block are, respectively, u2 (t) = (˜ u(t) + a2 λ2 γ2 y˘1 (t))w2 (t),

y2 (t) = L−1 {G2 (s)L{u2 (t)}}.

Thus, during the time interval 0 ≤ t < t2 , the whole output is the sum of the outputs from α1 -block and α2 -block, multiplied: during 0 ≤ t < t1 by γ1 , and 26

during t1 ≤ t < t2 by γ2 , i.e., y(t) = (y2 (t) + y˘1 (t))wγ (t). For the consecutive switchings in time, we proceed analogously, obtaining, in general, iterated algorithm (38). Example 5. Using (38a), we calculate the solution of E-type differential equation (35), where u(t) = H(t), λ(t) = 1, and          2 for t ∈ [0, 1), 1 for t ∈ [0, 1), −1 for t ∈ [0, 1),          α(t) = 1 for t ∈ [1, 2), γ(t) = 2 for t ∈ [1, 2), a(t) = −2 for t ∈ [1, 2),               −3 for t ∈ [2, 10), 2 for t ∈ [2, 10), 3 for t ∈ [2, 10), and then

G1 (s) =

1 1 1 ˘ 1 (s) = 1 , G ˘ 2 (s) = 1 . , G2 (s) = , G3 (s) = 2 , G s2 + 1 s+4 s +9 s2 s

The solution of (35) is

200

   1 − cos(t),    
  1 2−2t  sinh (−2 + 2t) + e−4t+4 23 − 2 cos (1) ,  2 + sin (1) e    

  3 sin (1) 2e2−2t sinh (−2 + 2t) − e2 sinh (−4 + 2t) + 4 y(t) =

 3 −4t+4   1 − e4 + 21 cos (−7 + 3t) 1 − 3e−4 −  4 sin (1) e   

 1 1   cos (−5 + 3t) 3e−4 + 1 + 12 cos (−6 + 3t) 27e−4 + 5 +  2      3 1 − e−4
(sin (−5 + 3t) − sin (−7 + 3t)) + 1 8 3

for t ∈ [0, 1), for t ∈ [1, 2),

for t ∈ [2, 10).

Plots of analytical and numerical solutions of differential equation (35) are

depicted in Fig. 17, and the difference plot between these solutions is presented in Fig. 18, where the numerical step size is 0.005.

27

Figure 17: Solution plots (solid line—analytical, dotted line—numerical (for step size 0.005)) of differential equation (35) from Ex. 5.

Figure 18: Plot of difference between analytical and numerical solutions from Ex. 5.

Taking a(t) = 0 in (35), we get E α(t) x 0 Dt

= λ(t)u,

x(0) = 0,

y(t) = γ(t)x(t),

(39a) (39b)

which yields an integral equation y(t) = γ(t)B D−α(t) (λ(t)u) .

(40)

Corollary 3. The solution of (40), and thereby of (39), is given by (38a),

28

where ˘ i (s)L{˜ yi (t) = L−1 {G u(t)wi (t)}} n o ˘ j (s)L{˜ y˘j (t) = L−1 G u(t)wj (t)} ,

(41a) j = 1, . . . , i − 1,

(41b)

where u ˜(t) is given by (20c). 205

Remark 8. Solution (41) is equal to the solution of i-a scheme, both in serial and parallel form, given by (12) and (14), respectively. Example 6. Consider  1     2 α(t) = 3 4     1 6

The solution is

y(t) =

integral system (40), where λ(t) = 1, γ(t) = 1 and    t for t ∈ [0, 2), for t ∈ [0, 2),    u(t) = 0 for t ∈ [2, 6), for t ∈ [2, 6),      for t ∈ [6, 10), 1 for t ∈ [6, 10).

 3/2   43 t√π ,    3/2

for t ∈ [0, 2),

√

−(1+t) t−2 4t √ , 3 π   √ √ √  3/2   4 πt −(1+t)π t−2+ 49 πΓ( 56 ) 6 t−6 3

π 3/2

for t ∈ [2, 6), for t ∈ [6, 10).

Plots of analytical and numerical solutions of differential equation (40) are depicted in Fig. 19, and the difference plot between these solutions is presented in Fig. 20, where the numerical step size is 0.005.

Figure 19: Solution plots (solid line—analytical, dotted line—numerical (for step size 0.005)) of integral equation (40) from Ex. 6.

29

Figure 20: Plot of difference between analytical and numerical solutions from Ex. 6.

210

3.1.4. Solution of B-type differential equation Let us consider the Tℓ -type differential equation (17) for ℓ = 4, i.e., B α(t) x 0 Dt

= λ(t) (a(t)γ(t)x + u) ,

x(0) = 0,

y(t) = γ(t)x(t).

(42a) (42b)

By (19), it can be rewritten as −α(t)

y(t) = γ(t)E0 Dt

(λ(t)(a(t)γ(t)x + u)) .

(43)

We will use the following objects (for i = 1, . . . , N ): Gi (s) =

1 , s αi + a i λ i γ i

where αi = αi−1 + α ˆi, and ˆ i (s) = 1 G sαˆ i

w ˜i = H(t−ti−1 )−H(t−ti ),

wi = H(t−ti )

u ˜i (t) = λ(t)u(t)w˜i .

Therefore, we have the following result. Proposition 4. Solution of B-type differential equation (42), on time interval 0 ≤ t ≤ ti , i = 1, . . . , N , is the following y(t) =

N X i=1

30

yi (t),

(44a)

where yi (t) = γi y˘i (t)    j i−1   X Y ˜i (s) + ˆ j (s)−1  Gi (s) Yˆj (s) y˘i (t) = L−1 U G   j=1 k=1 n o yˆj (t) = L−1 Y˜j (s) wi (t) Yˆj (s) = L{ˆ yj (t)} Y˜j = L {ei−1 (t)}

j Y

ˆ k (s) G

(44b) (44c) (44d) (44e)

k=1

ei−1 (t) = (˜ ui−1 (t) + ai−1 λi−1 γi−1 y˘i−1 (t)) w ˜i−1 (t).

(44f)

Proof. Considering the system shown in Fig. 8 (for ℓ = 4 and j = 1, X4 = o-a and F1 = ser) containing the output-additive switching scheme (see Fig. 4) in forward path it is possible to find the solution of (42). At the beginning, in time interval 0 ≤ t < t1 , the system described above

ˆ 1 (s) in forward path, can be treated as a fractional constant order system with G and λ1 , a1 , γ1 in feedback loop. Then, the output is as following Y˘1 (s) = G1 (s)U1 (s), G1 (s) =

ˆ 1 (s) G , ˆ 1 (s)λ1 a1 γ1 1−G

Y1 (s) = γ1 Y˘1 (s),

y1 (t) = γ1 y˘1 (t).

ˆ 2 (s) integral block is added in series to G ˆ 1 (s) At time instant t1 , the second G block and starts to integrate with 0 initial value. In the same time, the first ˆ 1 (s) block keeps integrating from the value where it finished in the previous G ˆ 1 (s) time interval. To obtain the initial value of first fractional-order integrator G let us begin with e1 (t) = (˜ u1 + a1 λ1 γ1 y˘1 (t)) w ˜1 (t). ˆ 1 (s) block on time interval 0 ≤ t < t1 and its constant initial The output of G value on time interval t1 ≤ t < t2 is expressed by y˜1 (t) = L−1 {Y˜1 (s)},

ˆ 1 (s)L{e1 (t)}, Y˜1 (s) = G

yˆ1 (t) = L−1 {Y˜1 (s)}w1 (t). 31

ˆ 1 (s) block in front of sum block, the output signal Moving the initial value of G ˆ 2 (s) block on time interval t1 ≤ t < t2 is of G

  ˜2 (s) + G ˆ 1 (s)−1 Yˆ1 (s) G2 (s), Y˘2 (s) = U

Yˆ1 (s) = L−1 {ˆ y1 (t)}.

Finally, the solution of (42) on time interval t1 ≤ t < t2 leads us to y˘2 (t) = L−1 {Y˘2 (s)}.

y2 (t) = γ2 y˘2 (t),

Analogously as above, the solution of (42) on time interval t2 ≤ t < t3 can

ˆ 3 (s) integral block is added in series be found. At time instant t2 , the third G ˆ 2 (s) block and starts to integrate with 0 initial value. However, the first to G ˆ 1 (s) and second G ˆ 2 (s) blocks continue to integrate with the values where G these integrators finished in previous sub-interval. To obtain the initial values for single fractional order integrator we start with e2 (t) = (˜ u2 + a2 λ2 γ2 y˘2 (t)) w ˜2 (t). ˆ 1 (s) block on time interval t1 ≤ t < t2 and its constant initial The output of G value on time interval t2 ≤ t < t3 is expressed by y˜1 (t) = L−1 {Y˜1 (s)},

ˆ 1 (s)L{e2 (t)}, Y˜1 (s) = G

yˆ1 (t) = L−1 {Y˜1 (s)}w2 (t). ˆ 2 (s) block on time interval t2 ≤ t < t3 Similarly, the constant initial value of G is expressed by y˜2 (t) = L−1 {Y˜2 (s)},

ˆ 1 (s)G ˆ 2 (s)L{e2 (t)}, Y˜2 (s) = G

yˆ2 (t) = L−1 {Y˜2 (s)}w2 (t). ˆ 1 (s) and G ˆ 2 (s) integrators in front of sum Then, moving the initial value of G ˆ 3 (s) integrator on time interval t2 ≤ t < t3 is block, the output signal of G   ˜3 (s) + G ˆ 1 (s)−1 Yˆ1 (s) + G ˆ 1 (s)−1 G ˆ 2 (s)−1 Yˆ2 (s) G3 (s), Y˘3 (s) = U Yˆ2 (s) = L−1 {ˆ y2 (t)}.

32

Finally, the solution of (42) on time interval t2 ≤ t < t3 leads us to y3 (t) = γ3 y˘3 (t), 215

y˘3 (t) = L−1 {Y˘3 (s)}.

For the consecutive switches in time, we proceed analogously, obtaining, in general, iterated algorithm (44). Example 7. Using (44a), we calculate the solution of B-type differential equation (42), where u(t) = H(t), λ(t) = 1, and          −6 4 for t ∈ [0, 1), 1 for t ∈ [0, 1),                   −7 5 for t ∈ [1, 2), 2 for t ∈ [1, 2), a(t) = γ(t) = α(t) =          −8 6 for t ∈ [2, 3), 3 for t ∈ [2, 3),                −9 7 for t ∈ [3, 4), 4 for t ∈ [3, 4),

for t ∈ [0, 1), for t ∈ [1, 2), for t ∈ [2, 3), for t ∈ [3, 4),

and then

G1 (s) =

1 , s + 24

G2 (s) =

s2

1 , + 35

G3 (s) =

s3

1 , + 48

G4 (s) =

s4

1 , + 63

ˆ 1 (s) = G ˆ 2 (s) = G ˆ 3 (s) = G ˆ 4 (s) = 1 . G s The solution of (42) is 
 1 −24t  , 6 1 − e    √ √ √


2   35  35(t − 1) 1 − e−24 + 27 sin 21 35 (t − 1)  168 sin       0.125 − 0.044e7.28−3.64t         −e1.82t−3.64 (0.0005 cos (3.15t − 6.3) sin (3.15t − 6.3))      0.11 − cos (2t − 6) 0.111 cosh(2t − 6) + 0.018e6−2t     y(t) = −0.018e−6+2t + 0.0107e−2t+2.34 + 0.02e−2t+7.8       −0.01e2t−9.62 − 0.02e2t−4.16 + 0.0017e−21.7−2t     
  −0.0017e2t−33.6 + sin (2t − 6) 0.018e6−2t +         0.018e−6+2t − 0.01e−2t+2.34 − 0.009e7.8−2t        −0.01e−9.6+2t + 0.015e−4.16+2t − 0.0018e−21.7−2t       −0.0018e−33.6+2t + 0.00006 sinh (2t − 6)
33

for t ∈ [0, 1), for t ∈ [1, 2),

for t ∈ [2, 3)

for t ∈ [2, 3).

Plots of analytical and numerical solutions of differential equation (35) are depicted in Fig. 21, and the difference plot between these solutions is presented in Fig. 22, where the numerical step size is 0.001.

Figure 21: Solution plots (solid line—analytical, dotted line—numerical (for step size 0.001)) of differential equation (35) from Ex. 7.

Figure 22: Plot of difference between analytical and numerical solutions from Ex. 7.

Taking a(t) = 0 in (42), we get B α(t) x 0 Dt

= λ(t)u,

x(0) = 0,

y(t) = γ(t)x(t),

(45a) (45b)

which yields an integral equation y(t) = γ(t)E D−α(t) (λ(t)u) . 34

(46)

220

Corollary 4. The solution of (46), and thereby of (45), is given by (44a), where Gi ’s are taken with ai = 0. Example 8. Consider an integral system (46), where λ(t) = 1, γ(t) = 1, and     1 for t ∈ [0, 1),  1 for t ∈ [0, 1), 2 u(t) = α(t) =   1 for t ∈ [1, 2). 1 for t ∈ [1, 2), The solution is  √   2√πt ,    p y(t) = 2.25 · (0.32t − 0.32)      +0.32 · (−LerchPhi(1 − t, 1, 1.5)t + 2)(t − 1)3/2

for t ∈ [0, 1),

for t ∈ [1, 2).

Plots of analytical and numerical solutions of differential equation (46) are

depicted in Fig. 23, and the difference plot between these solutions is presented 225

in Fig. 24, where the numerical step size is 0.001.

Figure 23: Solution plots (solid line—analytical, dotted line—numerical (for step size 0.001)) of integral equation (46) from Ex. 8.

35

Figure 24: Plot of difference between analytical and numerical solutions from Ex. 8.

4. Conclusions In the paper, an original approach for analytical solving of linear fractional variable order equations was proposed. For different types of variable order derivative definitions, considerably divergent algorithms have been obtained. 230

The proposed methods are based on switching schemes approach and the duality properties between variable order derivatives. In order to confirm the proposed methodology, results of analytical solutions were compared with numerical implementation. Presented results fully confirm efficiency and correctness of the proposed methods. Differences between analytical and numerical solutions, pre-

235

sented in this paper, are mostly due to accuracy of used numerical algorithm for solving variable order differential equations rather than accuracy of the analytical method itself. Although the proposed methods provide analytical solutions, in general they may lead to non-trivial symbolical computations. These methods can be helpful, e.g., in analysis of properties of variable order systems and,

240

especially, to validate numerical algorithms for solving these types of equations. References [1] A. Dzieli´ nski, D. Sierociuk, Fractional order model of beam heating process and its experimental verification, in: D. Baleanu, Z. B. Guvenc, J. A. T. Machado (Eds.), New Trends in Nanotechnology and Fractional Calculus

245

Applications, Springer Netherlands, 2010, pp. 287–294. [2] D. Sierociuk, A. Dzieli´ nski, G. Sarwas, I. Petras, I. Podlubny, T. Skovranek, Modelling heat transfer in heterogeneous media using fractional calculus, 36

Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 371 (1990). 250

[3] P. A. Harris, R. Garra, Nonlinear heat conduction equations with memory: Physical meaning and analytical results, Journal of Mathematical Physics 58 (6) (2017) 063501. doi:doi.org/10.1063/1.4984583. [4] I. Colombaro, A. Giusti, F. Mainardi, A class of linear viscoelastic models based on bessel functions, Meccanica 52 (4) (2017) 825–832.

255

doi:10.1007/s11012-016-0456-5. [5] A. Giusti, On infinite order differential operators in fractional viscoelasticity, Fractional Calculus and Applied Analysis 20 (4) (2017) 854–867. doi:10.1515/fca-2017-0045. [6] M. R. Kumar, S. Ghosh, S. Das, Charge-discharge energy efficiency analysis

260

of ultracapacitor with fractional-order dynamics using hybrid optimization and its experimental validation, AEU - International Journal of Electronics and Communications 78 (2017) 274 – 280. [7] J.-D. Gabano, T. Poinot, H. Kanoun, Lpv continuous fractional modeling applied to ultracapacitor impedance identification, Control Engineering

265

Practice 45 (2015) 86 – 97. [8] A. Dzieli´ nski,

G. Sarwas,

D. Sierociuk,

Time domain validation

of ultracapacitor fractional order model, in:

Decision and Con-

trol (CDC), 2010 49th IEEE Conference on, 2010, pp. 3730–3735. doi:10.1109/CDC.2010.5717093. 270

[9] A. Dzieli´ nski, G. Sarwas, D. Sierociuk, Comparison and validation of integer and fractional order ultracapacitor models, Advances in Difference Equations 2011:11. doi:{10.1186/1687-1847-2011-11}. [10] R. Garrappa, F. Mainardi, G. Maione, Models of dielectric relaxation based on completely monotone functions, Fractional Calculus and Applied Anal-

275

ysis 19 (5) (2016) 1105–1160. doi:{10.1515/fca-2016-0060}. 37

[11] D. Sierociuk, G. Sarwas, M. Twardy, Resonance phenomena in circuits with ultracapacitors, in: Proceedings of 12th International Conference on Environment and Electrical Engineering (EEEIC), 2013, pp. 197–202. doi:10.1109/EEEIC.2013.6549616. 280

[12] P. Ostalczyk, T. Rybicki, Variable-fractional-order dead-beat control of an electromagnetic servo, Journal of Vibration and Control 14 (9-10) (2008) 1457–1471. [13] A. Dabiri, B. Moghaddam, J. T. Machado, Optimal variable-order fractional pid controllers for dynamical systems, Journal of Computational and

285

Applied Mathematics. [14] K. Miller, B. Ross, An Introduction to the Fractional Calculus and Fractional Differenctial Equations, John Wiley & Sons Inc., New York, 1993. [15] C. A. Monje, Y. Chen, B. M. Vinagre, D. Xue, V. Feliu, Fractional-order Systems and Controls, Springer, 2010.

290

[16] K. B. Oldham, J. Spanier, The Fractional Calculus, Academic Press, 1974. [17] I. Podlubny, Fractional Differential Equations, Academic Press, 1999. [18] S. Samko, A. Kilbas, O. Maritchev, Fractional Integrals and Derivative. Theory and Applications, Gordon & Breach Sci. Publishers, 1987. [19] H. Sheng, Y. Chen, T. Qiu, Signal Processing Fractional Processes and

295

Fractional-Order Signal Processing, Springer, London, 2012. [20] M. D. Ortigueira, J. T. Machado, What is a fractional derivative?, Journal of Computational Physics 293 (2015) 4 – 13, fractional PDEs. doi:https://doi.org/10.1016/j.jcp.2014.07.019. [21] R. Gorenflo, F. Mainardi, Fractional Calculus: Integral and Differential

300

Equations of Fractional Order, Springer Vienna, Vienna, 1997, pp. 223– 276. doi:10.1007/978-3-7091-2664-6_5.

38

[22] Z. Li, H. Wang, R. Xiao, S. Yang, A variable-order fractional differential equation model of shape memory polymers, Chaos, Solitons and Fractals 102 (2017) 473 – 485, future Directions in Fractional Calculus Research and 305

Applications. doi:https://doi.org/10.1016/j.chaos.2017.04.042. [23] D. Yin, P. Qu, Variable-order fractional MSD function to describe the evolution of protein lateral diffusion ability in cell membranes, Physica A: Statistical Mechanics and its Applications 492 (2018) 707 – 714. doi:https://doi.org/10.1016/j.physa.2017.10.030.

310

[24] H. Sheng, H. Sun, C. Coopmans, Y. Chen, G. W. Bohannan, Physical experimental study of variable-order fractional integrator and differentiator, in: Proceedings of The 4th IFAC Workshop Fractional Differentiation and its Applications FDA’10, 2010. [25] L. Ramirez, C. Coimbra, On the variable order dynamics of the nonlinear

315

wake caused by a sedimenting particle, Physica D-Nonlinear Phenomena 240 (13) (2011) 1111–1118. [26] P. Sakrajda, D. Sierociuk, Modeling heat transfer process in grid-holes structure changed in time using fractional variable order calculus, in: A. Babiarz, A. Czornik, J. Klamka, M. Niezabitowski (Eds.), Theory and

320

Applications of Non-integer Order Systems, Springer International Publishing, Cham, 2017, pp. 297–306. [27] C.

Lorenzo,

T.

Hartley,

fractional operators,

Variable

order

and

distributed

order

Nonlinear Dynamics 29 (1-4) (2002) 57–98.

doi:{10.1023/A:1016586905654}. 325

[28] D. Valerio, J. S. da Costa, Variable-order fractional derivatives and their numerical approximations, Signal Processing 91 (3, SI) (2011) 470–483. doi:{10.1016/j.sigpro.2010.04.006}. [29] L. E. S. Ramirez, C. F. M. Coimbra, On the selection and meaning of

39

variable order operators for dynamic modeling, International Journal of 330

Differential Equations. [30] D. Sierociuk, W. Malesza, M. Macias, Derivation, interpretation, and analog modelling of fractional variable order derivative definition, Applied Mathematical Modelling 39 (13) (2015) 3876–3888. [31] D. Sierociuk, W. Malesza, M. Macias, On the recursive fractional variable-

335

order derivative: Equivalent switching strategy, duality, and analog modeling, Circuits, Systems, and Signal Processing 34 (4) (2015) 1077–1113. doi:10.1007/s00034-014-9895-1. [32] M. Macias, D. Sierociuk, An alternative recursive fractional variable-order derivative definition and its analog validation, in: ICFDA’14 International

340

Conference on Fractional Differentiation and Its Applications 2014, 2014, pp. 1–6. doi:10.1109/ICFDA.2014.6967452. [33] M. Zaky,

E. Doha,

T. Taha,

D. Baleanu,

New recursive ap-

proximations for variable-order fractional operators with applications, Mathematical Modelling and Analysis 23 (2) (2018) 227–239. 345

doi:https://doi.org/10.3846/mma.2018.015. [34] D. Sierociuk, W. Malesza, M. Macias, Equivalent switching strategy and analog validation of the fractional variable order derivative definition, in: Proceedings of European Control Conference 2013, ECC’2013, Zurich, Switzerland, 2013, pp. 3464–3469.

350

[35] D. Sierociuk, M. Macias, W. Malesza, Analog modeling of fractional switched-order derivatives: Experimental approach, in: Advances in the Theory and Applications of Non-integer Order Systems, Springer International Publishing, 2013, pp. 271–280. [36] D. Sierociuk, M. Twardy, Duality of variable fractional order difference op-

355

erators and its application to identification, Bulletin of the Polish Academy of Sciences: Technical SciencesSubmitted to. 40

[37] Y. Xu, Z. He, Existence and uniqueness results for cauchy problem of variable-order fractional differential equations, Applied 360

Mathematics

and

Computing

43

(1-2)

Journal of

(2013)

295–306.

doi:10.1007/s12190-013-0664-2. [38] S. Zhang, Existence result of solutions to differential equations of variableorder with nonlinear boundary value conditions, Communications in Nonlinear Science and Numerical Simulation 18 (12) (2013) 3289 – 3297. [39] C. Wen, Z. Jianjun, Z. Jinyang, A variable-order time-fractional deriva-

365

tive model for chloride ions sub-diffusion in concrete structures, Fractional Calculus and Applied Analysis 16 (1) (2013) 76–92. [40] E. Demirci, N. Ozalp, A method for solving differential equations of fractional order, Journal of Computational and Applied Mathematics 236 (11) (2012) 2754 – 2762.

370

[41] W. Malesza, D. Sierociuk, M. Macias, Solution of fractional variable order differential equation, in: 2015 American Control Conference (ACC), 2015, pp. 1537–1542. doi:10.1109/ACC.2015.7170951. [42] A. Bhrawy, M. Zaky, Highly accurate numerical schemes for multidimensional space variable-order fractional Schr¨odinger equations, Com-

375

puters & Mathematics with Applications 73 (6) (2017) 1100 – 1117, advances in Fractional Differential Equations (IV): Time-fractional PDEs. doi:https://doi.org/10.1016/j.camwa.2016.11.019. [43] A. H. Bhrawy, M. A. Zaky, Numerical algorithm for the variable-order caputo fractional functional differential equation, Nonlinear Dynamics 85 (3)

380

(2016) 1815–1823. doi:10.1007/s11071-016-2797-y. [44] Y.-T. Jia, M.-Q. Xu, Y.-Z. Lin, A numerical solution for variable order fractional functional differential equation, Applied Mathematics Letters 64 (2017) 125 – 130. doi:https://doi.org/10.1016/j.aml.2016.08.018.

41

[45] A. H. Bhrawy, M. A. Zaky, Numerical simulation for two-dimensional 385

variable-order fractional nonlinear cable equation, Nonlinear Dynamics 80 (1) (2015) 101–116. doi:10.1007/s11071-014-1854-7. [46] M. D. Ortigueira, J. J. Trujillo, Generalized GL Fractional Derivative and its Laplace and Fourier Transform, in: Proceedings of ASME International Design Engineering Technical Conferences and Computers and Information

390

in Engineering Conference, Vol. 4, ASME, Design Engn Div; ASME, Comp & Informat Engn Div, 2010, pp. 1227–1231. [47] W. Malesza, D. Sierociuk, Analytical description and equivalence of additive-switching scheme for fractional variable-order continuous differ-integrals,

395

in:

ICFDA’14 International Conference on Frac-

tional Differentiation and Its Applications 2014, doi:10.1109/ICFDA.2014.6967453.

42

2014,

pp. 1–6.