On Partial Differential Equations that Exhibit Fractional Behaviors

On Partial Differential Equations that Exhibit Fractional Behaviors J. Sabatier, H. C. Nguyen, X. Moreau, A. Oustaloup IMS Laboratory, LAPS/CRONE Group CNRS UMR 5218, Bordeaux 1 University 351 cours de la libération, 33405 Talence, France Tel: +33-5-40-00-66-07, Fax: +33-5-40-00-66-44 Email: {firstname.lastname}@ims-bordeaux.fr Abstract – The link between fractional differentiation and diffusion equations is reminded in the paper. The system infinite dimension along with the constant geometric characteristics is these equations are at the origin of the fractional behavior. This paper demonstrates that several other classes of differential equations also exhibit, on a frequency band, a fractional behavior. The fractional behavior is obtained with these equations on a space of finite dimension but with particular geometric characteristics. Keywords – Partial differential equations, fractional systems

1. INTRODUCTION The link between fractional systems, recursivity and diffusion equations is now well known. This link is reminded in the second section of the paper. This link should be used to implement fractional integration in software dedicated to numerical solving of differential equation such as for instance COMSOL Multiphysics software. However and as shown in section 2, the diffusion equation form of a fractional system requires the computation of an inverse Fourier transform that is in most case impossible to compute analytically. This is why this paper proposes alternative differential equations that exhibits fractional behaviours in a given frequency band. These differential equations can be easily implemented to simulate a fractional real or complex integrator. As also shown, they also permit to represent real mechanical systems.

ht

sinJS f J tx ³ x e dx .

S

(3)

0

Response of system (1) to an input u(t) is defined as the convolution product of the impulse response h(t) and the input u(t): y t

t

³ ht  W u W dW ,

(4)

0

and thus using relation (3) and through an integral permutation: y t

f sin

³0

JS x J §¨ t e  t W xu W dW ·¸dx . ¨³ ©0

S

(5)

¸ ¹

Let wt , x

t

³e

t W x

u W dW .

(6)

0

2. FRACTIONAL SYSTEMS, DIFFUSION AND RECURSIVITY

the following state space representation can be obtained for system (1):

The following fractional system (fractional integrator) is considered. H s s

J

,

(1)

with 0  J  1 . Its link with diffusion equation can be demonstrated through the system impulse response (Oustaloup, 1983). This impulse response is defined by the Mellin-Fourier integral of (1):

^ `

ht L1 s J

1 Z of 2 jS lim

c  jf st J ³c  jf e s ds

(2)

where c is greater than the abscissa of the singular points of H(s). Using poles definition that can be found in (Oustaloup, 1983) (Sabatier et al, 2011), this system do not generate poles and its impulse response is thus given by :

­ wwt , x  xwt , x  u t ° wt ® sin JS f J ° y t x wt , x dx S ³0 ¯

x  R .

(7)

Such a representation can be generalised to a large class of fractional systems as demonstrated in (Montseny, 1998) (Matignon, 1998). In these works, second relation in (7) is rewritten as: y t

f

³ P x wt , x dx ,

(8)

0

and representation representation.

(7)-(8)

is

denoted

diffusive

Initial conditions are defined for this system (7) by w0, x U x and thus permits give the exact expression of the system response with initial conditions (Sabatier et al, 2010): f

§

0

©

³ P x ¨¨ w0, x e

 xt

t ·  ³ e  x t W u W dW ¸ dx . ¸ 0 ¹

(9)

Through several changes of variables described in (Sabatier et al, 2010), system (1), (but also a large number of fractional systems) can be described by: ­ wI t , ] w I t , ]  u t G ] ° w] 2 ° wt . (10) ® f ° y t ³ m] I t , ] d] ° f ¯ Relation (10) shows that a fractional integrator can thus be seen as an infinite dimensional system described by a diffusion equation. This remark can be generalized to large number of fractional systems and thus demonstrates their link with diffusion equations.

1.208 Poles and zeros recursivity

y t

1.21

Zeros Poles

1.206

1.204

1.202

1.2

2

1

2

10

10 Frequency (rd/s

3

10

4

10

Fig. 2. Poles and zeros recursivity of Hd(s): Q = 0.5, 'z=0.182, N1 = 0 ( Z N1 0 rd / s ), N2=51 ( Z N2

10000 rd / s )

Note that this truncation can also be done on the second relation of (10) like that: y t

]

] 0

]1

³ m] I t , ] d]

] 2

^



.

(14)

`

f

where F 1 In (14) m] F 1 4S 2 x P 4S 2 x 2 denotes the inverse Fourier transform. This relation is in practice not easy to compute analytically. Alternative differential equations are now proposed for real and complex fractional integrators.

f

3. DIFFERENTIAL EQUATIONS CONSIDERED

ut

I t , ]

y t

1.198 0 10

³ m] I t , ] d]

Fig. 1. Representation of system (10) and (1) Computing Laplace transform of relation (3) and using P(x) function previously introduced, the following transfer function for a fractional system is obtained: f

H s

P x

dx , ³ 0 px

or using the change of variable x H s

f

³

(11)

 dz .

(12)

p f 1 z e

¦



wxz, t w 2 x z , t  J z wz wtwz

p



e  k'z

.

x0, t  xZ , t .

(16)

E 1 z s z and J z g 1 z

(13)

Let Z k e  k'z . Relation (13) highlights the recursivity of Hd(s) poles. Due to truncation with a finite N1 and N2 required for its implementation, figure 3 shows that the zeros have also a recursive distribution in a frequency band located between Z  N 2 ...Z N1 . This thus highlights the link of fractional systems with recursivity.

(15)

In (17), coefficients E(z) and Jz) are supposed to be exponential laws of the form :

1

@

u t .

Also, let the system output y(t) be given by : s t

P e  k'z 'z

N1 of k  N 2 N 2 of

>

dimension ( z  >0..Z @ , Z  R  ) and of the time variable t. This function satisfies the class of partial differential equations:

E z

P e z

N1

lim

Let x(z, t) be a function of the space variable z of finite

e z

The discrete form of H(s) is H d s

3.1 Differential equations definition

with s z s 0 e  Az

and g z

g 0 e Bz

(17)

Transfer function that link the system input and output is defined by:

S  p E p

E 1 z ³ 1  J z E 1 z s dz

Z

s k 1 sk

0

s 0 e  Az dz ³ s0  A B z 0 1 s e g 0

Z

.

(18)

Using s0 g 0

B , G A B

Q

and c

e  A B Z , (19)

D 'z K 'z ,

(26)

1.21

integral H1(s) can be expressed as a sum of two hypergeometric functions (Erdelyi, 1955) (Oustaloup et al, 1995):

H 1 s

Z k 1 1 and Zk K 'z

thus showing the recursivity of H1d(s) poles. The recursivity of the zeros is represented in figure 3. This figure highlights an equivalent recursivity with figure 2. This recursivity, used in (Oustaloup 1983) for the synthesis of fractional integrator operators, is another proof of the link of the considered differential equation with fractional differentiation.

ª sG · º Q 1 § ¸ » «sG  1 F ¨1 Q ,1 Q ,2 Q , sG  1 ¹ » s 0 « © . A « sG · » Q 1 § ¸» « sG  c F ¨1 Q ,1 Q ,2 Q , sG  c ¹ ¼ © ¬ (20)

1.208

Poles and zeros recursivity

H 1 s

Zeros Poles

1.206

1.204

1.202

1.2

1.198

Given that sG 1 sG  1

and

sG 1 , sG  c

s  C , (21)

2

3

10 Fréquence (rd/s)

10

4

10

D('z)K('z)=1.2, K('z)=1,095; Z1=1 rd/s; ZN=10000 rd/s

As constants A and B are considered positive, parameter Q is fractional and ranges from 0 to 1. When the transitional frequency 1/G is far enough from the transitional frequency c/G, H1(s) behaves as K0sQ within interval [1/G, c/G], The gain decreases at a rate of 6(1-Q) dB/octave and the phase is constant and equal to (Q)S/2. This result reveals the link between the class of partial differential equations considered and fractional differentiation. Another way to highlight the link with fractional differentiation is the analysis of the poles and zeros of the discrete form of H1(s) defined by:

s

1

10

Fig. 3. Poles and zeros recursivity of H1d(s): Q = 0.5,

relation (20) converges on C.

H 1d

1.196 0 10

s0 e  Ak'z 'z , Z s¦ s 0  A B k'z k 11  s e g 0 N

N'z . (22)

3.2. An example of system As an example, the mechanical behaviour of the connection of two shaft represented by figure 4 is now studied. This connection shaft is composite material (carbon-epoxy) made. As shown in figure 4, shaft section variation imposes a variation of the winding angle Tdefined as the angle formed by the fibers in relation to the longitudinal axis. T is a linear function of z, namely:

T z J z  Z ' ,

where Z' is greater than Z and denotes the fillet length defined by :

Introducing the recursive factors D('z) and K('z) such that s  z  'z 1  A'z sz

1 , K 'z

g  z  'z 1  B'z g z

D 'z

1

,

(27)

Z' Z 

T z . J

(28)

Carbon fiber

(23)

(24)

0

T

z S1

relation (22) becomes N

H 1d s s ¦

sk

sk

k 11 

s

Zk

,

Zk

In (25), it is obvious to check that

s0 e  Ak'z 'z g 0  A B k'z . (25) e s0

S0 Fig. 4. Carbon-epoxy shaft connection Angle T varies on the interval T = [10°, 30°]. On this interval, Young modulus E(T(z)) is an exponential function defined by:

E T z E 0 e O1T  z ,

(29)

4.1. Class 2 Let x(z, t) be a function of the space variable z of finite

with: E0

O1

and

180.5 MPa

0.0569 .

(30)

On this interval, damping coefficient <(T(z)) is also an exponential function given by: < T z <0 e

O2T  z

,

(31)

dimension ( z  >0..Z @ , Z  R  ) and of the time variable t. This function satisfies the class of partial differential equations:

G 2 z

w 3 xz , t wt 2 wz

 G 1 z

with: <0

4.15%

O1

and

0.082 .

(32)

System output is supposed defined by:

Moreover, section variation is also an exponential function of the form: S z S 0 e O3T  z ,

(33)

s t

O3

· ¸¸ . ¹

(34)

H 2 s

1 , E T z S z

1 E0 S 0 e

(35)

O1JL '

e O1J O3 z

.

G , < T z E T z S z

(37)

G

<0 E 0 S 0 e O2 O1 JZ ' e O1J O2J O3 z

.

(43)

D 1 z

(44)

s 0

e E0 S 0

g 0

O1J  O3

and

B

r r

2

r

2

2

r i

i

2

2

(45) (46)

G 2 z E z 2 br 2  bi 2

(47)



then S  p E p

H 2 s

E z dz ,  0 1 s E z J z

Z

³

(38)

(48)

using

E z E z br  ibi

(49)

J z J z g r  ig i .

(50)

E z E 0 e  Az br  ibi

(51)

J z J 0 e Bz g r  ig i .

(52)

Ge O1 O2 JZ ' <0 E 0 S 0

Let now (39)

A

2

G 1 z 2J z E z br g r  bi g i

Using O1JZ '

 E z J z g b  g b z J z g  g 

where G is a constant and thus given (29), (31) and (33): g z

Z

D 0 z J z 2 E Z br g r 2  br g i 2

G0

(36)

Moreover, linear viscous sliding coefficient is given by: g z

sD 1  z  D 0  z ³ s 2G z  sG z  G z dz . (42) 0 2 1 0

S  p E p

G0(z), G1(z) et G2(z) are defined by:

or given relations (29) and (33) : sz

(41)

Moreover, it is supposed that coefficients D0(z), D1(z),

Linear flexibility is defined by: sz

x0, t  xZ , t .

Transfer function that links the system input and output is defined by:

with: §S 1 log¨¨ 0 Z © S1

w 2 x z , t wxz , t  G 0 z wtwz wz . (40) wu t D 1 z  D 0 z u t wz

O3  O1J  O 2 J ,

linear flexibility and viscous sliding coefficients expressions are similar to relations (17). The shaft behaviour can then be modelled using transfer function (18). 4. OTHER CLASSES OF SYSTEMS This section proposes others classes of differential equation which behaves as a fractional real or complex integrator.

Coefficients G(z), G(z), G(z), D(z) et Dz) (relations 43 to 47) thus become:

 D z E 0 J 0 g b G z J 0 g

e 

D 0 z J 0 2 E 0 br g r 2  br g i 2 e 2 B  A z (53) 2

1

0

r r

2

r

2

 g r bi 2

2

 g i 2 e 2 Bz



B  2 A z

(54) (55)

G 1 z 2J 0 E 0 br g r  bi g i e B  A z

(56)

G 2 z E 0 2 br 2  bi 2 e 2 Az

(57)





Discrete form of H2(s) is defined by: H 2 d s

N

E 0 e  Ak'z br  ibi 'z , E 0 br  ibi  A B k'z 11  s e J 0 g r  ig i

¦

k

Z

N 'z .

4.2. Class 3 In equation (40), we suppose now that coefficients

D0(z), D1(z), G0(z), G1(z) et G2(z) are defined by: (58)

Introducing the recursive factors D('z) and K('z) such that 1 E z  'z , (59) 1  A'z K  'z E z

J z  'z 1  B'z J z

1

D 'z

,

(60)

H 2 d s

¦

Ek

k 11 

E 0 br  ibi e  Ak'z 'z J 0 g r  ig i  A B k'z . e E 0 br  ibi

Ek ,

s

Zk

Zk

(61) In (61), it is obvious to check that

E k 1 Ek

Z k 1 1 and Zk K 'z

D 'z K 'z ,

D 1 z J r

(63)

2

i

G 0  z J r  z 2  J i  z 2

(64)

G 1 z 2E r z J r z  E i z J i z

(65)

G 2  z E r  z 2  E i  z 2 .

(66)

E  z E r  z  iE i  z

(67)

J  z J r  z  iJ i  z ,

(68)

H 3 s

S  p E p

ª º «Z » E z Re / i « ³ dz » . « E z » «0 1 s » J z ¼ ¬

(69)

If

Z i+1

Zi

Z i-1

E z E 0 e  A z

E 0 e  Ar z e iTz

E 0 e  Ar z e

Re



i Tz T E 0



> 





E 0 e  Ar z cos Tz  T E 0  i sin Tz  T E 0

DK

DK

J z J 0 e B z

Fig. 5. H2d(s) poles recursivity The zeros recursivity is represented in figure 6. This figure highlights an equivalent recursivity with figure 2. 1.21

J 0 e Br z e iTz

> 

J 0 e Br z e







i Tz TJ 0

J 0 e Br z cos Tz  T J 0  i sin Tz  T J 0

@

@

, (71)

 E 0 2 J 0 cosTz  T J e B  2 A z

D 0 z E 0 J 0 2 cos Tz  T E 0 e  A 2 B z D 1 z

Zeros Poles

1.206

,(70)

coefficients G(z), G(z), G(z), D(z) et Dz) (relations (62) to (66)) become:

1.208

Poles and zeros recursivity

 z 2

then transfer function that links the system input and output is defined by:

thus showing the recursivity of H2d(s) poles. This recursivity is illustrated in the complex plane by figure 5.

1.204

0

G 0 z J 0 2 e 2 Bz



1.202

1.2

1.198

1.196 0 10

r

Let

(2)

Im

 E z

(62)

relation (58) becomes N

 z E

D 0  z E r  z J r  z 2  J i  z 2

(72) (73) (74)



G 1 z 2J 0 E 0 cos T E 0  T J 0 e 2 Bz

(75)

G 2 z E 0 2 e 2 Az .

(76)

Discrete form of transfer function (69) is given by 1

10

2

10 Fréquence (rd/s)

3

10

4

10

H 3 d s

Fig. 6. Ratio between two consecutive zeros for :

D('z)K('z) = 1,2, K('z) = 1,095; Z1 = 1+i rd/s; Z N = 10000+10000i rd/s with

ª «N E k Re / i « ¦ s «k 1 « 1 Z k ¬

º » », » » ¼

(77)

Ek Zk

> 



@



E 0 e  Ar k'z cos Tk'z  T E 0  i sin Tk'z  T E 0 'z J 0 cos Tk'z  T J 0  i sin Tk'z  T J 0  A B k'z .

>   @ e E 0 >cosTk'z  T E  i sin Tk'z  T E @ 0

0

(78) As shown by Nichols diagram of figure 7, transfer function transmittance H3d(s) behaves, on a bounded frequency interval, as complex fractional integrator. Its phase diagram is a linear function of the frequency. 15 10 5

Gain (dB)

0 -5 -10 -15 -20 -25 -30 -1200

-1000

-800

-600 Phase (°)

-400

-200

0

Fig. 7. Nichols diagram of H3d(s) for : Z1 =1+5i, Z N =1000+5000i, D('z)K ('z)=1.2,

K 'z 1.0955e iS 6 , D 'z 1.0955e iS

6

One can note that the poles and zeros recursivity obtained in H3d(s), is the one, used in (Oustaloup et al, 2000) for the synthesis of a complex fractional integrator, and is a proof of the link of the considered differential equation with fractional differentiation. 5. CONCLUSION This paper first reminds the link of fractional systems with diffusion equations and recursivity. It is shown that the diffusion equation representation is obtained through the system impulse response computation and several transformations involving inverse Fourier transform. The output of the system represented by the diffusion equation is then the weighed sum of the system state of infinite dimension (system spatial dimension is infinite). The implementation of a fractional system in a software dedicated to differential equation simulation (like Comsol Multiphisics) should be done using this representation along with a convenient truncation in space. However, the diffusion equation representation requires the computation of an inverse Fourier transform involving fractional terms, which is in practice impossible to obtain analytically in most cases. To overcome this problem, the paper proposes others differential equations whose behaviour is fractional within a frequency domain. In the present state of our study, the fractional behavior obtained is those of a fractional integrator, real and complex. This is proved through the introduction of hypergeometric function of through the analysis of poles and zeros recursivity.

In future works, we intend now to extend the number of fractional differentiations classes to cover a large number of fractional systems. REFERENCES Erdelyi A. (1955), Higher transcendental functions", Vol. 1, California Institute of Technology, Mac Graw-Hill. Levron F., Sabatier J., Oustloup A., Habsieger L. (1999), From partial differential equations of propagative recursive systems to non integer differentiation Fractional Calculus and Applied Analysis (FCAA): an international journal for theory and applications, Vol. 2, N° 3, pp. 246-264. Matignon D. (1998), Stability properties for generalized fractional differential systems, ESAIM: Proc. 5, pp 145-158. Montseny G. (1998), Diffusive representation of pseudodifferential time-operators, ESAIM: Proc. vol 5, pp 159-175. Oustaloup A. (1983), Systèmes asservis linéaires d’ordre fractionnaire. Editions Masson, Paris. Oustaloup A., Nouillant M., (1991), Distributed Parameter Systems", European Control Conference, ECC-91, Grenoble, France, July 2-5. Oustaloup A, Sabatier J. (1995), Recursive Distributed Parameter Systems and Non Integer Derivation European Conference on Circuit Theory and Design ECCTD'95, 27-31 August 1995, Istanbul, Turkey. Oustaloup, A. Levron, F. Mathieu, B. Nanot, F. (2000), Frequency-band complex noninteger differentiator: characterization and synthesis, IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications, Vol 47, N° 1, pp 25-39. Sabatier J., Merveillaut M., Malti R., Oustaloup A. (2010), How to Impose Physically Coherent Initial Conditions to a Fractional System ?, Communications in Nonlinear Science and Numerical Simulation, Volume 15, Issue 5.