# Solving Linear Coupled Fractional Differential Equations by Direct Operational Method and Some Applications

Authors: S. C. Lim; Chai Hok Eab; K. H. Mak; Ming Li; S. Y. Chen
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

A new direct operational inversion method is introduced for solving coupled linear systems of ordinary fractional differential equations. The solutions so-obtained can be expressed explicitly in terms of multivariate Mittag-Leffler functions. In the case where the multiorders are multiples of a common real positive number, the solutions can be reduced to linear combinations of Mittag-Leffler functions of a single variable. The solutions can be shown to be asymptotically oscillatory under certain conditions. This technique is illustrated in detail by two concrete examples, namely, the coupled harmonic oscillator and the fractional Wien bridge circuit. Stability conditions and simulations of the corresponding solutions are given.


---

## Body

## 1. Introduction

Fractional differential equations are well suited to model physical systems with memory or fractal attributes. This is particularly true in the fields of condensed matter physics, where fractional differential equations have been used to model various anomalous transport and relaxation phenomena [1–9]. Coupled fractional differential equations (CFDEs) of nonlinear type are widely used in studying various chaotic systems [10] such as the Lorentz system [11], fractional Chuah’s circuit [12], fractional Rössler system [13], and fractional Duffing system [14]. Since in most cases no analytic solutions for such nonlinear CFDEs exist, it is necessary to resort to numerical approximations and simulations [15–21]. Even for linear CFDEs with unequal multiorders, analytic solutions, if they exist, are difficult to obtain and very often numerical methods have to be used.In this paper we introduce a direct operational method to solve a system of linear inhomogeneous CFDEs. We will restrict our discussion to a system of linear nonhomogeneous ordinary differential equations of arbitrary fractional-orders. These equations based on two types of fractional derivatives will be considered, namely, the Caputo and Riemann-Liouville fractional derivatives. The main idea is to reexpress the coupled fractional equations by incorporating the initial conditions based on the definitions of these derivatives. The solutions obtained can be expressed in terms of multivariate Mittag-Leffler functions. When each order of the CFDEs is an integer multiple of a certain common real positive number, it is possible to further reduce the solutions to the single-variate Mittag-Leffler functions. For such cases, we study the conditions for the existence of asymptotic periodic solutions.In the next section we consider two types of coupled fractional differential equations of Caputo and Riemann-Liouville type, and a direct operational method is introduced to solve these equations. Subsequent sections deal with the applications of the coupled fractional differential equations to two physical systems, namely, the coupled fractional oscillator and the fractional Wien bridge circuit, as examples to illustrate the proposed method.


## 2. Linear-Coupled Fractional Differential Equations

We consider two types of fractional derivatives [22–26]:(2.1a)CaputoD*αf(t)=Im-αDmf(t),(2.1b)Riemann-LiouvilleD#αf(t)=DmIm-αf(t),where the fractional integral is defined for γ>0 as(2.2)Iγf(t)=1Γ(γ)∫0∞(t-τ)γ-1f(τ)dτ.
When referring to either definition, we simply use the notation Dα.Let us consider a linear-coupled system of inhomogeneous fractional differential equations of the form(2.3)DαX(t)=B(t)+AX(t),
where X=(x1,…,xn) and B=(b1,…,bn) are vectors of dimension n, A=(aij), i,j=1,…,n is an (n×n)-matrix, and Dα is the fractional differential operator given by the diagonal matrix operator:
(2.4)Dα=(Dα10⋯00Dα2⋯0⋮⋮⋱⋮00⋯Dαn).
The orders αi are real positive numbers with mi-1<αi<mi, where mi is a positive integer for each i=1,2,…,n. The boundary conditions for (2.3) are given by(2.5a)[D*kixi](0)=c*ikiforki=0,1,2,…,(mi-1),(2.5b)[D#αi-kixi](0)=c#ikiforki=1,2,…,mi,and i=1,2,…,n. We remark that the mi’s generally differ in value. Let its maximum and minimum values be denoted by m and mo, respectively.For a single fractional differential equation, its solution can be obtained by integral transform methods such as the Fourier, Laplace, and Mellin transforms (see, e.g., [24, Chapter 4]). However, in the case of a system of CFDEs, it is necessary to employ specific techniques appropriate to the given problem, that is, the form of matrix A and the type of fractional derivative involved. There exist several methods (see [24, Chapters 5 and 6]) for solving such problems. Here we want to develop a technique which is more direct, similar to Green’s function method.Let the operator be(2.6)L=Dα-A.
The solution of (2.3) can then be expressed as
(2.7)X(t)=L-1B.
Unfortunately, the inverse of L may not exist. However, the right-inverse G does exist for both the Riemann-Liouville and Caputo cases with
(2.8)LG=1≠GL.
The main task now is to determine the right-inverse of L associated with (2.3) for both Riemann-Liouville and Caputo fractional derivatives. We remark that our treatment is rather formal, aiming mainly to provide an alternative direct operational method to the usual Laplace transform technique in solving CFDEs. In particular, the existence of solutions will be assumed, and all operators considered are assumed to be well defined in a certain appropriate function space.


### 2.1. The Right-Inverse Operator

For the Caputo derivativeD*α with arbitrary α and m,(2.9a)IαD*αf(t)=IαIm-αDmf(t)=ImDmf(t)=f(t)-∑k=0m-1tkΓ(k+1)[Dkf](0),
and similarly for the Riemann-Liouville derivative,
(2.9b)IαD#αf(t)=IαDmIm-αf(t)=f(t)-∑k=1mtα-kΓ(α-k+1)[D#α-kf](0).Applying the fractional integral operator Iα to (2.3) and using the initial conditions (2.5a) and (2.5b) gives
(2.10)X=IaAX+IaB+W,
where W is given by
(2.11)w*i=∑k=0mi-1tkΓ(k+1)c*ik,w#i=∑k=1mitα-kΓ(α-k+1)c#ik.
Let
(2.12)Q=[1-IαA].
Now by rearranging (2.10) one gets
(2.13)QX=IαB+W.
The operator Q has an inverse K. One possible representation of K is given by
(2.14)K=∑0∞(IαA)p.
This form may not be a simple one, since in most cases Iα and A do not commute. However, the verification is straightforward:
(2.15)QK=KQ=∑p=0∞(IαA)p-∑p=0∞(IαA)p+1=1.
The other possible representation is
(2.16)K=ΨQ*,
where Q* is the adjoint of Q, that is, Q*Q=detQ, and the Ψ is the inverse operator of detQ such that ΨdetQf(t)=f(t). It is quite simple to verify that this representation is the inverse of Q, and this will be shown in Section 4.Now we define(2.17)G=KIα.
One can easily verify that G is the right inverse of L. Thus, the solution is given by
(2.18)X=GB+KW.
The detailed evaluation can be carried out by using different techniques, a few of which will be considered here.


### 2.2. System with Constant Inhomogeneous Terms

When the inhomogeneous term is a constant, we can absorb this term in the following way, though it may not be immediately obvious. For the Caputo case, let(2.19)X*=X̅*-A-1B,
with
(2.20)x̅*i(0)=x*i(0)+(a-1)ijbj.
The initial conditions have to be transformed accordingly and they become
(2.21)c̅*i0=c*i0+(a-1)ijbj.In the Riemann-Liouville case, we cannot absorb the termB into X̅ as in the Caputo case. However, if one compares the initial condition terms w#i and the term IαB, one can see that if we modify
(2.22)w̅#i=∑k=0mitα-kΓ(α-k+1)c#ik,
with c#i0=bi, then the solution can be written as if it is a solution of a homogeneous linear equation.Thus the inhomogeneous linear fractional differential equation with constant source termbi can be transformed into a homogeneous linear fractional differential equation. The solution of the transformed equation can then be written as
(2.23)X̅=KW̅.
When bi is time dependent with power up to mi, the above modification can still apply to the Caputo case. In the Riemann-Liouville case, if bi(·) are analytic functions, then the same modification as above holds if (2.22) is altered with the summation from -∞ to mi-1.In subsequent sections, we will consider the solution of (2.3) according to the above modifications.


## 2.1. The Right-Inverse Operator

For the Caputo derivativeD*α with arbitrary α and m,(2.9a)IαD*αf(t)=IαIm-αDmf(t)=ImDmf(t)=f(t)-∑k=0m-1tkΓ(k+1)[Dkf](0),
and similarly for the Riemann-Liouville derivative,
(2.9b)IαD#αf(t)=IαDmIm-αf(t)=f(t)-∑k=1mtα-kΓ(α-k+1)[D#α-kf](0).Applying the fractional integral operator Iα to (2.3) and using the initial conditions (2.5a) and (2.5b) gives
(2.10)X=IaAX+IaB+W,
where W is given by
(2.11)w*i=∑k=0mi-1tkΓ(k+1)c*ik,w#i=∑k=1mitα-kΓ(α-k+1)c#ik.
Let
(2.12)Q=[1-IαA].
Now by rearranging (2.10) one gets
(2.13)QX=IαB+W.
The operator Q has an inverse K. One possible representation of K is given by
(2.14)K=∑0∞(IαA)p.
This form may not be a simple one, since in most cases Iα and A do not commute. However, the verification is straightforward:
(2.15)QK=KQ=∑p=0∞(IαA)p-∑p=0∞(IαA)p+1=1.
The other possible representation is
(2.16)K=ΨQ*,
where Q* is the adjoint of Q, that is, Q*Q=detQ, and the Ψ is the inverse operator of detQ such that ΨdetQf(t)=f(t). It is quite simple to verify that this representation is the inverse of Q, and this will be shown in Section 4.Now we define(2.17)G=KIα.
One can easily verify that G is the right inverse of L. Thus, the solution is given by
(2.18)X=GB+KW.
The detailed evaluation can be carried out by using different techniques, a few of which will be considered here.


## 2.2. System with Constant Inhomogeneous Terms

When the inhomogeneous term is a constant, we can absorb this term in the following way, though it may not be immediately obvious. For the Caputo case, let(2.19)X*=X̅*-A-1B,
with
(2.20)x̅*i(0)=x*i(0)+(a-1)ijbj.
The initial conditions have to be transformed accordingly and they become
(2.21)c̅*i0=c*i0+(a-1)ijbj.In the Riemann-Liouville case, we cannot absorb the termB into X̅ as in the Caputo case. However, if one compares the initial condition terms w#i and the term IαB, one can see that if we modify
(2.22)w̅#i=∑k=0mitα-kΓ(α-k+1)c#ik,
with c#i0=bi, then the solution can be written as if it is a solution of a homogeneous linear equation.Thus the inhomogeneous linear fractional differential equation with constant source termbi can be transformed into a homogeneous linear fractional differential equation. The solution of the transformed equation can then be written as
(2.23)X̅=KW̅.
When bi is time dependent with power up to mi, the above modification can still apply to the Caputo case. In the Riemann-Liouville case, if bi(·) are analytic functions, then the same modification as above holds if (2.22) is altered with the summation from -∞ to mi-1.In subsequent sections, we will consider the solution of (2.3) according to the above modifications.


## 3. System with Equal Fractional Orders

In the case where allαi=α, then Iα=Iα1, and
(3.1)K=∑p=0∞AnInα.
It would be convenient if we introduce the Mittag-Leffler function with matrix argument [27–29]:
(3.2)Eα,β(Z)=∑n=0∞ZnΓ(nα+β).
Then
(3.3)Kδ(t)=t-1Eα,0(Atα).
Here we have used the following definition of the Dirac delta function:
(3.4)δ(t)=lim⁡ε→0t+ε-1Γ(ε).
The matrix A can always be decomposed into Jordan normal form. However, we consider only the case where it can be diagonalized:
(3.5)Λ=P-1AP,K̂=P-1KP,
with eigenvalues λj as the diagonal elements of Λ, and
(3.6)K̂ijδ(t)=δijt-1Eα,0(λjtα),(3.7)X*=∑k=0m-1PtkEα,k+1(Λtα)P-1C*k,X#=∑k=-1m-1Ptα-k-1Eα,α-k(Λtα)P-1C#k.
Thus, the solutions of the fractional differential equation under consideration based on both types of fractional derivatives are simply linear combinations of Mittag-Leffler functions.


### 3.1. Coupled Oscillator with Equal Fractional Orders

Here we demonstrate our method by considering a linear-coupled oscillator system given by(3.8)Dα1x1(t)=-ω2x1(t)+ω̅2(x2(t)-x1(t)),Dα2x2(t)=-ω2x2(t)+ω̅2(x1(t)-x2(t)),
where ω2 and ω̅2 are nonnegative real numbers. The initial conditions are
(3.9)x*j(0)=c*j0,Dx*j(0)=c*j1,Dαj-1x#j(0)=c#j0,Dαj-2x#j(0)=c#j1.
Thus
(3.10)A=-(ω2+ω̅2-ω̅2-ω̅2ω2+ω̅2).
For simplicity, we use the following notation:
(3.11)a11=a22=-(ω2+ω̅2)=-η,a12=a21=ω̅2=ε.
We first consider the simpler case with α1=α2=α in this section. Since A is symmetric, it can be diagonalized to give two eigenvalues:
(3.12)λ±=-η±ε={-ω2-ω2-ω̅2.
Both eigenvalues are nonpositive and
(3.13)P=12(11-11).
Just as we have shown in (3.7), the solutions are again linear combinations of Mittag-Leffler functions:(3.14a)x*1=(c*10+c*20)Eα,1(λ-tα)+(c*10-c*20)Eα,1(λ+tα)+(c*11+c*21)tEα,2(λ-tα)+(c*11-c*21)tEα,2(λ+tα),(3.14b)x#1=(c#11+c#21)tα-1Eα,α(λ-tα)+(c#11-c#21)tα-1Eα,α(λ+tα)+(c#12+c#22)tα-2Eα,α-1(λ-tα)+(c#12-c#22)tα-2Eα,α-1(λ+tα).Figure 1 shows simulations of the Caputo solution (3.14a) for orders 1<α≤2. An interesting feature of fractional oscillators in general is the presence of damping internal to the system, that is, an inherent decay in the amplitude which is not associated with any external friction. The variation in the amount of internal damping can be clearly seen as the order increases. In the limiting cases we obtain exponential decay (α=1) and undamped oscillation (α=2).Coupled oscillator with equal fractional orders. Parameters:ω=1, ω̅=0.5, x*1(0)=1.0, x*2(0)=0, Dx*1(0)=0, Dx*2(0)=0.1. Legend: x*1 (solid line), x*2 (dashed line).
(a)(b)(c)(d)


## 3.1. Coupled Oscillator with Equal Fractional Orders

Here we demonstrate our method by considering a linear-coupled oscillator system given by(3.8)Dα1x1(t)=-ω2x1(t)+ω̅2(x2(t)-x1(t)),Dα2x2(t)=-ω2x2(t)+ω̅2(x1(t)-x2(t)),
where ω2 and ω̅2 are nonnegative real numbers. The initial conditions are
(3.9)x*j(0)=c*j0,Dx*j(0)=c*j1,Dαj-1x#j(0)=c#j0,Dαj-2x#j(0)=c#j1.
Thus
(3.10)A=-(ω2+ω̅2-ω̅2-ω̅2ω2+ω̅2).
For simplicity, we use the following notation:
(3.11)a11=a22=-(ω2+ω̅2)=-η,a12=a21=ω̅2=ε.
We first consider the simpler case with α1=α2=α in this section. Since A is symmetric, it can be diagonalized to give two eigenvalues:
(3.12)λ±=-η±ε={-ω2-ω2-ω̅2.
Both eigenvalues are nonpositive and
(3.13)P=12(11-11).
Just as we have shown in (3.7), the solutions are again linear combinations of Mittag-Leffler functions:(3.14a)x*1=(c*10+c*20)Eα,1(λ-tα)+(c*10-c*20)Eα,1(λ+tα)+(c*11+c*21)tEα,2(λ-tα)+(c*11-c*21)tEα,2(λ+tα),(3.14b)x#1=(c#11+c#21)tα-1Eα,α(λ-tα)+(c#11-c#21)tα-1Eα,α(λ+tα)+(c#12+c#22)tα-2Eα,α-1(λ-tα)+(c#12-c#22)tα-2Eα,α-1(λ+tα).Figure 1 shows simulations of the Caputo solution (3.14a) for orders 1<α≤2. An interesting feature of fractional oscillators in general is the presence of damping internal to the system, that is, an inherent decay in the amplitude which is not associated with any external friction. The variation in the amount of internal damping can be clearly seen as the order increases. In the limiting cases we obtain exponential decay (α=1) and undamped oscillation (α=2).Coupled oscillator with equal fractional orders. Parameters:ω=1, ω̅=0.5, x*1(0)=1.0, x*2(0)=0, Dx*1(0)=0, Dx*2(0)=0.1. Legend: x*1 (solid line), x*2 (dashed line).
(a)(b)(c)(d)


## 4. The Adjoint Method

As we mentioned in Section2.1, (2.16), the main task now is to calculate the inverse of detQ. Note that (2.16) can be reexpressed as
(4.1)K=ψ[1-IαA]*,
so that
(4.2)G=ψ[1-IαA]*Iα.
Also recall that
(4.3)det⁡Q=det⁡(1-IαA).
We know that any finite dimension determinant can be evaluated; however, it can be easily obtained only for a few lower-dimensional cases. We will compute it explicitly for a two-dimensional system.


### 4.1. Two-Dimensional System

In a two-dimensional system the determinant is easy to calculate, and we have(4.4)det⁡Q=1-a11Iα1-a22Iα2+Iα1+α2det⁡A.
The inverse is given by
(4.5)ψ=∑r=0∞[a11Iα1+a22Iα2-Iα1+α2det⁡A]r(4.6)=∑r=0∞∑k1+k2+k3=rr!k1!k2!k3![a11Iα1]k1[a22Iα2]k2[-Iα1+α2det⁡A]k3.
Its kernel is the multivariate Mittag-Leffler function of the kind given by (B.3), Appendix B:
(4.7)ψδ(t)=t-1Ẽα1,α2,α1+α2,0(a11tα1,a22tα2,-det⁡Atα1+α2)=ε̃α1,α2,α1+α2,0(a11,a22,-det⁡A:t).
The adjoint of Q is
(4.8)Q*=(1-a22Iα2a12Iα1a21Iα21-a11Iα1).
In the following we evaluate explicitly only for the 11- and 12-elements, while the 22- and 21-elements can be obtained by just interchanging subscripts. The kernels of K are given by
(4.9)K11δ(t)=ε̃α1,α2,α1+α2,0(a11,a22,-det⁡A:t)-a22ε̃α1,α2,α1+α2,α2(a11,a22,-det⁡A:t),K12δ(t)=a12ε̃α1,α2,α1+α2,α1(a11,a22,-det⁡A:t).


### 4.2. The Solutions

The solutions based on the adjoint method are given by(4.10a)x*1=∑k=0m1-1c*1k[ε̃α1,α2,α1+α2,k+1(a11,a22,-det⁡A:t)-a22ε̃α1,α2,α1+α2,α2+k+1(a11,a22,-det⁡A:t)]+a12∑k=0m2-1c*2kε̃α1,α2,α1+α2,α1+k+1(a11,a22,-det⁡A:t),(4.10b)x#1=∑k=0m1-1c#1k[ε̃α1,α2,α1+α2,α1-k+1(a11,a22,-det⁡A:t)-a22ε̃α1,α2,α1+α2,α2+α1-k+1(a11,a22,-det⁡A:t)]+a12∑k=0m2-1c#2kε̃α1,α2,α1+α2,α1+α2-k+1(a11,a22,-det⁡A:t).


## 4.1. Two-Dimensional System

In a two-dimensional system the determinant is easy to calculate, and we have(4.4)det⁡Q=1-a11Iα1-a22Iα2+Iα1+α2det⁡A.
The inverse is given by
(4.5)ψ=∑r=0∞[a11Iα1+a22Iα2-Iα1+α2det⁡A]r(4.6)=∑r=0∞∑k1+k2+k3=rr!k1!k2!k3![a11Iα1]k1[a22Iα2]k2[-Iα1+α2det⁡A]k3.
Its kernel is the multivariate Mittag-Leffler function of the kind given by (B.3), Appendix B:
(4.7)ψδ(t)=t-1Ẽα1,α2,α1+α2,0(a11tα1,a22tα2,-det⁡Atα1+α2)=ε̃α1,α2,α1+α2,0(a11,a22,-det⁡A:t).
The adjoint of Q is
(4.8)Q*=(1-a22Iα2a12Iα1a21Iα21-a11Iα1).
In the following we evaluate explicitly only for the 11- and 12-elements, while the 22- and 21-elements can be obtained by just interchanging subscripts. The kernels of K are given by
(4.9)K11δ(t)=ε̃α1,α2,α1+α2,0(a11,a22,-det⁡A:t)-a22ε̃α1,α2,α1+α2,α2(a11,a22,-det⁡A:t),K12δ(t)=a12ε̃α1,α2,α1+α2,α1(a11,a22,-det⁡A:t).


## 4.2. The Solutions

The solutions based on the adjoint method are given by(4.10a)x*1=∑k=0m1-1c*1k[ε̃α1,α2,α1+α2,k+1(a11,a22,-det⁡A:t)-a22ε̃α1,α2,α1+α2,α2+k+1(a11,a22,-det⁡A:t)]+a12∑k=0m2-1c*2kε̃α1,α2,α1+α2,α1+k+1(a11,a22,-det⁡A:t),(4.10b)x#1=∑k=0m1-1c#1k[ε̃α1,α2,α1+α2,α1-k+1(a11,a22,-det⁡A:t)-a22ε̃α1,α2,α1+α2,α2+α1-k+1(a11,a22,-det⁡A:t)]+a12∑k=0m2-1c#2kε̃α1,α2,α1+α2,α1+α2-k+1(a11,a22,-det⁡A:t).


## 5. Laplace Transform Method

In this section we briefly discuss how the widely used Laplace transform technique can be employed to determine Green’s functionG and the operator K. There is no intention here to provide a detailed discussion of this method. Instead, it will be discussed as a complement to the adjoint method introduced above, so as to allow one to see the relation between the direct operational method presented here and usual Laplace transform method.Without loss of generality (see Section2.2) the system of CFDEs is assumed to be homogeneous. We begin by calculating the Laplace transform of detQδ(t) using (4.5):
(5.1)∫0∞det⁡Qδ(t)e-stdt=1-a11s-α1-a22s-α2+det⁡As-(α1+α2)
and the Laplace transform of the adjoint kernel Q*δ(t):
(5.2)∫0∞Q11*δ(t)e-stdt=1-a22s-α2,∫0∞Q12*δ(t)e-stdt=a12s-α1,∫0∞Q21*δ(t)e-stdt=a21s-α2,∫0∞Q22*δ(t)e-stdt=1-a11s-α1.
Next, the Laplace transforms of the part related to the initial conditions from (2.11) are
(5.3)∫0∞w*i(t)e-stdt=∑k=0mi-1s-k-1c*ik,∫0∞w#i(t)e-stdt=∑k=1mis-αi+k-1c#ik.
Thus the Laplace transforms of the solutions become(5.4a)x̂*1(s)=[1-a22s-α2]∑k=0m1-1s-k-1c*1k+a12s-α1∑k=0m2-1s-k-1c*2k1-a11s-α1-a22s-α2+det⁡As-(α1+α2),(5.4b)x̂#1(s)=[1-a22s-α2]∑k=1m1s-α1+k-1c#1k+a12s-α1∑k=1m2s-α2+k-1c#2k1-a11s-α1-a22s-α2+det⁡As-(α1+α2).x̂*2(s) and x̂#2(s) can be obtained just by interchanging 1↔2. From the complexity of the Laplace transforms, one sees that it is virtually impossible to obtain the analytic solutions by direct application of the inverse Laplace transform. To obtain the solution of this type one has to use the Laplace transform of the multivariate Mittag-Leffler function [27–29], which then gives the identity for getting the Laplace inversion of (5.4a) and (5.4b). This is one main advantage of the direct operational inversion method proposed here as it will give the solution directly.Clearly, it is important that both methods produce equivalent solutions. This is verified explicitly in AppendixC for a 2-dimensional system.


## 6. Multiple Fractional-Order System

In physics and engineering problems the fractional-orders can often be approximated by rational numbers, that is,αi=pi/qi, for some pi,qi∈ℕ. Thus one gets αi=μi/q, where q is the least common multiple of q1,q2,…,qn with some μi∈ℕ. However, we can also consider the more general case with αj=μjα0, for some μj∈ℕ. Here, α0∈ℝ+ can be either rational or irrational. In this section we show how such a system of CFDEs with these multiple fractional-orders can be solved.Referring to (4.4), if we assign the symbol ξ=Iα0, the expansion of the determinant will be a polynomial of order μ=μ1+μ2+μ3+⋯+μn in ξ. By the fundamental theorem of algebra, it must have in general μ complex roots, that is, ζj for j=1,2,3,…,μ. Note that since all coefficients of the polynomial are real, if any root ζp is complex, its complex conjugate ζ̅p is also a root. That means that the complex roots occur in pairs. For convenience we write λj=1/ζj, for j=1,2,3,…,μ. We can then factorize the polynomial(6.1)det⁡Q=∏j=1μ(1-λjξ).
If all roots are distinct, the inverse can be written as a partial fraction:(6.2)ψ=∑j=1μhj1-λjξ=∑j=1μhj∑k=0∞λjkIkα0.


### 6.1. Two-Dimensional System

We will explore this method further for two dimensions. Using (4.5) and the adjoint in (4.8), the solution is given by
(6.3)x*1=∑j=1μ1+μ2hj{∑k=0m1-1c*1k[tkEα0,k+1(λjtα0)-a22tμ2α0+kEα0,μ2α0+k+1(λjtα0)]+a12∑k=0m2-1c*2ktμ1α0+kEα0,μ1α0+k+1(λjtα0)},x#1=∑j=1μ1+μ2hj{∑k=1m1c#1k[tμ1α0-kEα0,μ1α0-k+1(λjtα0)-a22t(μ1+μ2)α0-kEα0,(μ1+μ2)α0-k+1(λjtα0)]+a12∑k=1m2c#2kt(μ1+μ2)α0-kEα0,(μ1+μ2)α0-k+1(λjtα0)}.To simplify the problem we consider the case where0<max(μ1α0,μ2α0)≤1, that is, m1=m2=1. The extension to the general case as above is straightforward:
(6.4)x*1=∑j=1μ1+μ2hj{c*10[Eα0,1(λjtα0)-a22tμ2α0Eα0,μ2α0+1(λjtα0)]+a12c*20tμ1α0Eα0,μ1α0+1(λjtα0)},x#1=∑j=1μ1+μ2hj{c#11tμ1α0-1Eα0,μ1α0(λjtα0)-(a22c#11-a12c#21)t(μ1+μ2)α0-1Eα0,(μ1+μ2)α0(λjtα0)}.
Using the following formula:
(6.5)zqEα,qα+γ(z)=Eα,γ(z)-∑p=0q-1zpΓ(pα+γ),(6.4) can be written as
(6.6)x*1=f*1+g*1,x#1=f#1+g#1,
where(6.7a)f*1=∑jμ1+μ2hjf*1j,g*1=∑jμ1+μ2hjg*1j,(6.7b)f#1=∑jμ1+μ2hjf#1j,g#1=∑jμ1+μ2hjg#1j.(6.8a)f*1j=[c*10(1-a22λj-μ2)+a12c*20λj-μ1]Eα0,1(λjtα0),(6.8b)f#1j=[c#11λj1-μ1-(a22c10-a12c#21)λj1-μ1-μ2]tα0-1Eα0,α0(λjtα0).(6.9a)g*1j=a22c*10λj-μ2∑p=0μ2-1λjptpα0Γ(pα0+1)-a12c*20λj-μ1∑p=0μ1-1λjptpα0Γ(pα0+1),(6.9b)g#1j=-c#11λj1-μ1∑p=0μ1-2λjptpα0+α0-1Γ(pα0+α0)+(a22c#11-a12c#21)λj1-μ1-μ2∑p=0μ1+μ2-2λjptpα0+α0-1Γ(pα0+α0).


### 6.2. Solutions with Asymptotic Oscillations

It is clear from the previous expansion of the Caputo terms that for the solutionx*1, as t→∞, any possible oscillation arises from the pair of complex roots (λj,λj¯), while all (negative) real roots must result in an asymptotic decay. Note that each term in (6.9a) with p>0 will grow asymptotically as the power law tpα0. However, when we combine the terms and reexpress the equation explicitly as
(6.10)g*1=∑j=1μ1+μ2hjg*1j=a22c*10∑p=0μ2-1[∑j=1μ1+μ2hjλjp-μ2]tpα0Γ(pα0+1)-a12c*20∑p=0μ1-1[∑j=1μ1+μ2hjλjp-μ1]tpα0Γ(pα0+1),
for the particular case with μ1=1,μ2=2, one has
(6.11)g*1=a22c*10∑p=01[∑j=13hjλjp-2]tpα0Γ(pα0+1)-a12c*20[∑j=13hjλj-1]1Γ(1)(6.12)=a22c*10{[∑j=13hjλj-2]1Γ(1)+[∑j=13hjλj-1]tα0Γ(α0+1)}-a12c*20[∑j=13hjλj-1]1Γ(1).
The second term is the only term that grows asymptotically as ~tα0, which does not contribute since [∑j=13hjλj-1] is zero. The verification of this result is given here for the general case of n-roots. Let us consider general partial fractions:
(6.13)1(1-λ1x)(1-λ2x)⋯(1-λnx)=h1(1-λ1x)+h2(1-λ2x)+⋯+hn(1-λnx).
We have
(6.14.0)h1+h2+⋯+hn=1,
(6.14.1)h1(λ2+λ3+⋯λn)+h2(λ3+λ4+⋯λn+λ1)+⋯+hn(λ1+λ2+⋯λn-1)=0,
(⋮)⋮
(6.14.n−1)h1λ2λ3⋯λn+h2λ3λ4⋯λnλ1+⋯+hnλ1λ2⋯λn-1=0.
We can rewrite (6.14.n−1) as(6.15)h1λ1+h2λ2+h3λ3+⋯+hnλn=0.
Using (6.15) with n=3 in (6.12), we have
(6.16)g*1=a22c*10[∑j=13hjλj-2],
which is a constant. Thus for this particular case with the Caputo derivative one can have asymptotic oscillations. In general, CFDEs based on the Caputo derivatives will not oscillate asymptotically, since one cannot find any rule for the power law terms to cancel out. However, this is possible for some special cases under suitable conditions on the elements of A.Similar consideration can be given to Riemann-Liouville system. However, nowg#1j approaches a constant or possible zero as t→∞ if 0<max(μ1α0,μ2α0)≤1, and min(μ1,μ2)≤2. Thus the possible asymptotically stable oscillations are due to the term f#1j with its corresponding complex conjugate. If one looks at (6.8b), one can write explicitly the asymptotic expansion:
(6.17)tα0-1Eα0,α0(λjtα0)=λj(1-α0)/α0α0eλj1/α0t+tα0-1O(1tα0)⟶λj(1-α0)/α0α0eλj1/α0t+O(1t)⟶λj(1-α0)/α0α0eλj1/α0t.
One has to impose the condition that λj1/α0 be a purely imaginary number. Let us denote the inversion of the root by λj=|λj|eiθ; we must have
(6.18)θ=±α0π2.
Furthermore, all the other roots that will not give rise to oscillation must have a negative real part so that their contributions will be asymptotically zero.Assume that there is only one pair of complex roots that satisfies (6.18). The coefficient of the exponential in (6.17) after substitution into (6.8b) gives the jth term of (6.7b) as
(6.19)hj[c#11λj1-μ1-(a22c#11-a12c#21)λj1-μ1-μ2]λj(1-α0)/α0α0=rjeiφ,
and we then have the asymptotic solution:(6.20a)x1⟶rjei|λj|t+iφ+rje-i|λj|t-iφ=2rjcos⁡(|λj|t+φ).x2 can be evaluated in a similar way; it has the same period but different modulus and phase:
(6.20b)x2⟶2rj′cos⁡(|λj|t+φ′).We omit the determination of the roots for each system, which can be computed without any difficulty.The oscillation condition (6.18) was first derived by Matignon [30] who also showed that an identical condition existed for the eigenfunction Eα,1(z) of the Caputo derivative. This will be elaborated in the context of a physical system in Section 7.


## 6.1. Two-Dimensional System

We will explore this method further for two dimensions. Using (4.5) and the adjoint in (4.8), the solution is given by
(6.3)x*1=∑j=1μ1+μ2hj{∑k=0m1-1c*1k[tkEα0,k+1(λjtα0)-a22tμ2α0+kEα0,μ2α0+k+1(λjtα0)]+a12∑k=0m2-1c*2ktμ1α0+kEα0,μ1α0+k+1(λjtα0)},x#1=∑j=1μ1+μ2hj{∑k=1m1c#1k[tμ1α0-kEα0,μ1α0-k+1(λjtα0)-a22t(μ1+μ2)α0-kEα0,(μ1+μ2)α0-k+1(λjtα0)]+a12∑k=1m2c#2kt(μ1+μ2)α0-kEα0,(μ1+μ2)α0-k+1(λjtα0)}.To simplify the problem we consider the case where0<max(μ1α0,μ2α0)≤1, that is, m1=m2=1. The extension to the general case as above is straightforward:
(6.4)x*1=∑j=1μ1+μ2hj{c*10[Eα0,1(λjtα0)-a22tμ2α0Eα0,μ2α0+1(λjtα0)]+a12c*20tμ1α0Eα0,μ1α0+1(λjtα0)},x#1=∑j=1μ1+μ2hj{c#11tμ1α0-1Eα0,μ1α0(λjtα0)-(a22c#11-a12c#21)t(μ1+μ2)α0-1Eα0,(μ1+μ2)α0(λjtα0)}.
Using the following formula:
(6.5)zqEα,qα+γ(z)=Eα,γ(z)-∑p=0q-1zpΓ(pα+γ),(6.4) can be written as
(6.6)x*1=f*1+g*1,x#1=f#1+g#1,
where(6.7a)f*1=∑jμ1+μ2hjf*1j,g*1=∑jμ1+μ2hjg*1j,(6.7b)f#1=∑jμ1+μ2hjf#1j,g#1=∑jμ1+μ2hjg#1j.(6.8a)f*1j=[c*10(1-a22λj-μ2)+a12c*20λj-μ1]Eα0,1(λjtα0),(6.8b)f#1j=[c#11λj1-μ1-(a22c10-a12c#21)λj1-μ1-μ2]tα0-1Eα0,α0(λjtα0).(6.9a)g*1j=a22c*10λj-μ2∑p=0μ2-1λjptpα0Γ(pα0+1)-a12c*20λj-μ1∑p=0μ1-1λjptpα0Γ(pα0+1),(6.9b)g#1j=-c#11λj1-μ1∑p=0μ1-2λjptpα0+α0-1Γ(pα0+α0)+(a22c#11-a12c#21)λj1-μ1-μ2∑p=0μ1+μ2-2λjptpα0+α0-1Γ(pα0+α0).


## 6.2. Solutions with Asymptotic Oscillations

It is clear from the previous expansion of the Caputo terms that for the solutionx*1, as t→∞, any possible oscillation arises from the pair of complex roots (λj,λj¯), while all (negative) real roots must result in an asymptotic decay. Note that each term in (6.9a) with p>0 will grow asymptotically as the power law tpα0. However, when we combine the terms and reexpress the equation explicitly as
(6.10)g*1=∑j=1μ1+μ2hjg*1j=a22c*10∑p=0μ2-1[∑j=1μ1+μ2hjλjp-μ2]tpα0Γ(pα0+1)-a12c*20∑p=0μ1-1[∑j=1μ1+μ2hjλjp-μ1]tpα0Γ(pα0+1),
for the particular case with μ1=1,μ2=2, one has
(6.11)g*1=a22c*10∑p=01[∑j=13hjλjp-2]tpα0Γ(pα0+1)-a12c*20[∑j=13hjλj-1]1Γ(1)(6.12)=a22c*10{[∑j=13hjλj-2]1Γ(1)+[∑j=13hjλj-1]tα0Γ(α0+1)}-a12c*20[∑j=13hjλj-1]1Γ(1).
The second term is the only term that grows asymptotically as ~tα0, which does not contribute since [∑j=13hjλj-1] is zero. The verification of this result is given here for the general case of n-roots. Let us consider general partial fractions:
(6.13)1(1-λ1x)(1-λ2x)⋯(1-λnx)=h1(1-λ1x)+h2(1-λ2x)+⋯+hn(1-λnx).
We have
(6.14.0)h1+h2+⋯+hn=1,
(6.14.1)h1(λ2+λ3+⋯λn)+h2(λ3+λ4+⋯λn+λ1)+⋯+hn(λ1+λ2+⋯λn-1)=0,
(⋮)⋮
(6.14.n−1)h1λ2λ3⋯λn+h2λ3λ4⋯λnλ1+⋯+hnλ1λ2⋯λn-1=0.
We can rewrite (6.14.n−1) as(6.15)h1λ1+h2λ2+h3λ3+⋯+hnλn=0.
Using (6.15) with n=3 in (6.12), we have
(6.16)g*1=a22c*10[∑j=13hjλj-2],
which is a constant. Thus for this particular case with the Caputo derivative one can have asymptotic oscillations. In general, CFDEs based on the Caputo derivatives will not oscillate asymptotically, since one cannot find any rule for the power law terms to cancel out. However, this is possible for some special cases under suitable conditions on the elements of A.Similar consideration can be given to Riemann-Liouville system. However, nowg#1j approaches a constant or possible zero as t→∞ if 0<max(μ1α0,μ2α0)≤1, and min(μ1,μ2)≤2. Thus the possible asymptotically stable oscillations are due to the term f#1j with its corresponding complex conjugate. If one looks at (6.8b), one can write explicitly the asymptotic expansion:
(6.17)tα0-1Eα0,α0(λjtα0)=λj(1-α0)/α0α0eλj1/α0t+tα0-1O(1tα0)⟶λj(1-α0)/α0α0eλj1/α0t+O(1t)⟶λj(1-α0)/α0α0eλj1/α0t.
One has to impose the condition that λj1/α0 be a purely imaginary number. Let us denote the inversion of the root by λj=|λj|eiθ; we must have
(6.18)θ=±α0π2.
Furthermore, all the other roots that will not give rise to oscillation must have a negative real part so that their contributions will be asymptotically zero.Assume that there is only one pair of complex roots that satisfies (6.18). The coefficient of the exponential in (6.17) after substitution into (6.8b) gives the jth term of (6.7b) as
(6.19)hj[c#11λj1-μ1-(a22c#11-a12c#21)λj1-μ1-μ2]λj(1-α0)/α0α0=rjeiφ,
and we then have the asymptotic solution:(6.20a)x1⟶rjei|λj|t+iφ+rje-i|λj|t-iφ=2rjcos⁡(|λj|t+φ).x2 can be evaluated in a similar way; it has the same period but different modulus and phase:
(6.20b)x2⟶2rj′cos⁡(|λj|t+φ′).We omit the determination of the roots for each system, which can be computed without any difficulty.The oscillation condition (6.18) was first derived by Matignon [30] who also showed that an identical condition existed for the eigenfunction Eα,1(z) of the Caputo derivative. This will be elaborated in the context of a physical system in Section 7.


## 7. Wien Bridge System

In this section we apply the solution methods discussed earlier to model a fractional-order Wien bridge oscillator (Figure2). The Wien bridge is a common electronic circuit that can generate a sinusoidal output signal without requiring an oscillatory input. The resistor-capacitor pairs form a frequency-selective network, hence allowing the selection of output sine wave frequency by varying the circuit parameters. Ahmad et al. [31] first proposed a generalization of this circuit using fractional-order capacitors. Since the authors did not obtain the solutions explicitly, we briefly show how analytic solutions for such a system can be obtained within our present framework. Also, we show solutions based on both the Caputo and the Riemann-Liouville derivatives. (Note that reference [31] does not mention the type of derivative used; we were informed by one of the authors, Professor Ahmad, that they used the Riemann-Liouville fractional derivative.)Wien bridge oscillator: (a) circuit schematic with operational amplifier and (b) simplified circuit diagram for voltage analysis.
(a)(b)It is well known that a fractional differential equation of order0<α<1 is usually used to describe relaxation phenomena [2]. In the case of Wien bridge system, however, oscillation is achieved via the active elements and feedback provided in the circuit (see Section 7.3).In the following we use normalized voltagesxi=vi/Vsat where Vsat is the amplifier saturation voltage and time axes (normalized with respect to time constant τ=RC). Using basic circuit analysis, it can be shown that the capacitor voltages are related via a 2-dimensional CFDE:
(7.1)DαX=AX+B,
where(7.2a)A=(a-2-1a-1-1),B=(bb),(7.2b)(a,b)={(0,1),Kx1≥1,(K,0),-1<Kx1<1,(0,1),Kx1≤1. Here K is the amplifier gain (i.e., vo=Kv1). In the linear region of the amplifier, -1<Kx1<1 and (7.2a) simplifies to
(7.3)A=(K-2-1K-1-1),B=(00),
Thus we have to solve the homogeneous linear fractional-order differential system with
(7.4)a11=K-2,a21=K-1,a12=-1,a22=-1.
For fractional capacitors, the real orders are restricted to
(7.5)0<α1≤1,0<α2≤1.
We remark that the boundary conditions associated with the Caputo derivative seem more “physical” as they can be verified by experiments, whereas for the Riemann-Liouville case, the fractional derivative boundary conditions cannot be measured. However, the Riemann-Liouville operators are popular with mathematicians and theoretical physicists. We can write the initial conditions explicitly as
(7.6)x1(0)=c*10,x2(0)=c*20,D#α1-1x1(0)=c#11,D#α2-2x2(0)=c#21.


### 7.1. Solution Using the Adjoint Method

Substituting (7.4) into the solution (4.10a) we get for the Caputo case
(7.7)x*1=c*10[ε̃α1,α2,α1+α2,1(K-2,-1,-1:t)+ε̃α1,α2,α1+α2,α2(K-2,-1,-1:t)]-c*20ε̃α1,α2,α1+α2,α1(K-2,-1,-1:t),x*2=c*20[ε̃α1,α2,α1+α2,1(K-2,-1,-1:t)-(K-2)ε̃α1,α2,α1+α2,α1(K-2,-1,-1:t)]+(K-1)c*10ε̃α1,α2,α1+α2,α2(K-2,-1,-1:t).
Similarly, for the Riemann-Liouville case, substituting (7.4) into the solution (4.10b), we get
(7.8)x#1=c#11[ε̃α1,α2,α1+α2,α1(K-2,-1,-1:t)+ε̃α1,α2,α1+α2,α2+α1-1(K-2,-1,-1:t)]-c#21ε̃α1,α2,α1+α2,α1+α2-1(K-2,-1,-1:t),x#2=c#21[ε̃α1,α2,α1+α2,α2(K-2,-1,-1:t)-(K-2)ε̃α1,α2,α1+α2,α2+α1-1(K-2,-1,-1:t)]-c#11(K-1)ε̃α1,α2,α1+α2,α1+α2-1(K-2,-1,-1:t).


### 7.2. Solution Using the Laplace Transform

Substituting (7.4) into the Laplace transform solution (5.4a), we obtain(7.9a)x̂*1(s)=[1+s-α2]s-1c*10-s-α1-1c*201-(K-2)s-α1+s-α2+s-(α1+α2),(7.9b)x̂*2(s)=[1-(K-2)s-α1]s-1c*20+(K-1)s-α2-1c*101-(K-2)s-α1+s-α2+s-(α1+α2).Similarly, for the Laplace transform of the solution (5.4b),
(7.10)x̂#1(s)=[1+s-α2]s-α1c#11-s-(α1+α2)c#211-(K-2)s-α1+s-α2+s-(α1+α2),x̂#2(s)=[1-(K-2)s-α1]s-α2c#21+(K-1)s-(α1+α2)c#111-(K-2)s-α1+s-α2+s-(α1+α2).
In the following subsections, we study the conditions under which asymptotically stable oscillations are possible for a fractional Wien bridge oscillator and also present numerical simulations of the capacitor voltages.


### 7.3. Equal-Order Fractional Wien Bridge

The classical Wien bridge oscillator produces a stable sinusoidal output when its amplifier gainK=3. The amplitude and frequency of the sinusoid are a function of the initial capacitor voltages and the circuit time constant. For a fractional capacitor, the current-voltage relationship is dependent on both capacitor value and order; hence an additional degree of freedom is introduced into the Wien bridge circuit. We consider first the simple case with equal real fractional orders α1=α2=α. Equation (7.7) then simplifies to
(7.11)x*1=c*10ε̃α,α,2α,1(K-2,-1,-1:t)+(c*10-c*20)ε̃α,α,2α,α(K-2,-1,-1:t),x*2=c*20ε̃α,α,2α,1(K-2,-1,-1:t)+[c*10(K-1)-c*20(K-2)]ε̃α,α,2α,α(K-2,-1,-1:t).
To gain further insight into the system’s behaviour, it is advantageous to express the solution in a simpler form using only 1-parameter Mittag-Leffler functions Eα,1(λitα). From (7.9a),
(7.12)x̂*1(s)=sα-1c*10sα+(c*10-c*20)s2α-(K-3)sα+1=σ1sα-1sα-λ1+σ2sα-1sα-λ2,
where the σi∈ℂ are constants to be determined from partial fraction decomposition (see equivalent method in Section 6). The inverse transform yields
(7.13)x*1(t)=σ1Eα,1(λ1tα)+σ2Eα,1(λ2tα).
The solution for x*2(t) can be found in a similar manner. For brevity, we present only numerical simulations of the Caputo solutions (one plot of the Riemann-Liouville solution is presented for comparison). In order for the Wien bridge to produce sustained oscillations, we need to impose condition (6.18) on the complex roots λi. For the current system, this translates to the following expression for K:
(7.14)K=3+2cos⁡(απ2).
Hence, the amplifier gain K is no longer a constant as in the case of the classical Wien bridge but a function of capacitor order α. Simulations of (7.13) were plotted using Mathematica.Figure3 shows plots of x*1 and x*2 for α=0.3,0.5,0.7, and 1.0. There is a clear dependence of waveform amplitude on the fractional-order. The plot for α=1 corresponds to the classical Wien bridge (with ordinary capacitors) and is included for comparison. It is important to keep in mind that the time axes are normalized and oscillation frequency ωα actually varies with order as
(7.15)ωα=(RC)-1/α.
This can be seen in Figure 4(a). It is interesting to note that the values of resistance and capacitance have no effect on the output waveform amplitude (Figure 4(b)). Hence, the frequency of oscillation can be controlled by both the value C and order α of the capacitors. As noted in [31], a clear advantage of this is that high frequencies can be obtained by reducing the order of the capacitors rather than their value, which can remain sufficiently large.Caputo model of fractional Wien bridge with equal orders—comparison of phase and amplitude for capacitor voltages. Parameters:x*1(0)=x*2(0)=0.03. Legend: x*1 (solid line); x*2 (dashed line).
(a)(b)(c)(d)Variation of waveform characteristics forx1(t). Parameters: x*1(0)=x*2(0)=0.05, (a) Capacitor order affects both frequency and amplitude. Legend: α=0.3 (solid), 0.5 (dashed), 0.7 (dash-dotted), 1.0 (dotted) and (b) Time constant affects frequency while amplitude remains constant. Legend: RC=0.8 (solid), 0.9 (dashed), and 1.0 (dash-dotted), and 1.1 (dotted).
(a)(b)In Figure5 we see that the Riemann-Liouville solutions (7.8) are very similar to the Caputo solutions in terms of frequency and amplitude but differ in phase due to the second parameter β of the Mittag-Leffler function. Of particular concern is the fact that the former tends to diverge at the origin since the fractional initial conditions (2.5b) do not correspond to measurable physical quantities. This is an important distinction between the two definitions and has to be taken into consideration when modeling physical systems.Figure 5
Comparison ofx1(t) waveform for both derivative models. Solution for Riemann-Liouville model diverges as t→0. Parameters: x*1(0)=x*2(0)=0.03. Legend: Caputo (solid line), and Riemann-Liouville (dashed line).


### 7.4. Multiorder Fractional Wien Bridge

The Wien bridge circuit can be further generalized by allowing the orders to assume valuesαj=μjα as detailed in Section 6. We use the case α2=2α1 as a starting point. Substituting α1=α and α2=2α into (7.7), (7.9a) and (7.9b), we have
(7.16)x*1=c*10[ε̃α,2α,3α,1(K-2,-1,-1:t)+ε̃α,2α,3α,2α(K-2,-1,-1:t)]-c*20ε̃α,2α,3α,α(K-2,-1,-1:t),x*2=c*20[ε̃α,2α,3α,1(K-2,-1,-1:t)-(K-2)ε̃α,2α,3α,α(K-2,-1,-1:t)]+(K-1)c*10ε̃α,2α,3α,2α(K-2,-1,-1:t).
As with the equal-order bridge, we express the solution in Laplace domain and use partial fractions to obtain a more tractable form:
(7.17)x̂*1(s)=sα-1(s2α+1)c*10-sαc*20s3α-(K-2)s2α+sα+1=∑k=13σksα-1sα-λk,(7.18)x*1(t)=∑k=13σkEα,1(λktα).Only solutions with one negative real and two complex-conjugate roots will be of concern to our present discussion. To justify this, we note that the alternative case of three real roots is of no physical interest as it does not produce oscillatory solutions. With the exception of the exponentially decaying term (due to the negative real root, see asymptotic analysis in Section6.2), the solution is hence similar to the case with equal capacitor orders; that is, the output of the fractional Wien bridge can be expressed as a linear combination of Mittag-Leffler functions. Unfortunately, the relationship between K and α is not as simple as in the equal-order case (7.14). Although K is still a function of α, its form is sufficiently complex that a more convenient alternative is to define K implicitly, that is, find α=ϕ(K), and use polynomial curve-fitting as shown in Figure 6.Determination of amplifier gainK: (a) restrictions on possible values of K and α; (b) third-order least-squares approximation of K(α) for x∈(0,0.5).
(a)(b)Using the method of least-squares, we obtain a third-order approximation forK:
(7.19)K≈4.611-4.821α2+0.008α3.Two restrictions apply to the usable range ofK and α. The first is the requirement that the real root be negative. Plotting the denominator of (7.17) as a function of sα for various amplifier gains, we obtain the relationship in Figure 7. For 0<K<K0≈4.611, we have λ1∈ℝ- and λ2=λ̅3∈ℂ as required. Therefore, this serves as an upper limit to the amplifier gain. The lower limit can be determined by recalling that 0<α2=2α1<1 so that 0<α<0.5. The result of these restrictions is also shown in Figure 6.Figure 7
Properties of roots of (7.19) (shown as horizontal intercepts) for parameter K.Within the stipulated range, the values of gain calculated from the polynomial curve (7.19) are sufficiently accurate to create oscillatory solutions, as demonstrated in the following simulations (Figure 8).Caputo model of fractional Wien bridge withα2=2α1. The effect of the nonoscillatory term in (7.20) can be observed as an initial offset that decays asymptotically as t→∞. Parameters: x*1(0)=x*2(0)=0.03, (a) Comparison of phase and amplitude for capacitor voltages. (b) The limiting case of one ordinary capacitor and one semicapacitor (order 1/2). Legend: and x*1 (solid line), x*2 (dashed line).
(a)(b)As mentioned earlier, it is possible to extend this procedure to any case where one order is an integer multiple of the other. Consider the solution ofx*1 for α1=α, α2=υα:
(7.20)x̂*1(s)=sα-1(sυα+1)c*10-s(υ-1)αc*20s(υ+1)α-(K-2)sυα+sα+1=∑k=1υ+1σksα-1sα-λk,(7.21)x*1(t)=∑k=1μ+1σkEα,1(λktα).
Indeed, unless we are concerned about obtaining the actual Laplace solution in partial fraction form, the value of K that produces asymptotic oscillations can be found by simply studying the roots of the denominator in (7.20) (a polynomial in sα of order υ+1) and imposing suitable conditions as previously shown so that at least one pair of Mittag-Leffler terms has an eigenvalue that satisfies (6.18). For example, when υ=3, a possible solution contains 4 roots in 2 complex-conjugate pairs. One can adjust the amplifier gain such that the roots satisfy |θ|=α0π/2 for the first pair and |θ|>α0π/2 for the second pair, hence resulting in sustained oscillation and asymptotically decaying oscillation, respectively.


## 7.1. Solution Using the Adjoint Method

Substituting (7.4) into the solution (4.10a) we get for the Caputo case
(7.7)x*1=c*10[ε̃α1,α2,α1+α2,1(K-2,-1,-1:t)+ε̃α1,α2,α1+α2,α2(K-2,-1,-1:t)]-c*20ε̃α1,α2,α1+α2,α1(K-2,-1,-1:t),x*2=c*20[ε̃α1,α2,α1+α2,1(K-2,-1,-1:t)-(K-2)ε̃α1,α2,α1+α2,α1(K-2,-1,-1:t)]+(K-1)c*10ε̃α1,α2,α1+α2,α2(K-2,-1,-1:t).
Similarly, for the Riemann-Liouville case, substituting (7.4) into the solution (4.10b), we get
(7.8)x#1=c#11[ε̃α1,α2,α1+α2,α1(K-2,-1,-1:t)+ε̃α1,α2,α1+α2,α2+α1-1(K-2,-1,-1:t)]-c#21ε̃α1,α2,α1+α2,α1+α2-1(K-2,-1,-1:t),x#2=c#21[ε̃α1,α2,α1+α2,α2(K-2,-1,-1:t)-(K-2)ε̃α1,α2,α1+α2,α2+α1-1(K-2,-1,-1:t)]-c#11(K-1)ε̃α1,α2,α1+α2,α1+α2-1(K-2,-1,-1:t).


## 7.2. Solution Using the Laplace Transform

Substituting (7.4) into the Laplace transform solution (5.4a), we obtain(7.9a)x̂*1(s)=[1+s-α2]s-1c*10-s-α1-1c*201-(K-2)s-α1+s-α2+s-(α1+α2),(7.9b)x̂*2(s)=[1-(K-2)s-α1]s-1c*20+(K-1)s-α2-1c*101-(K-2)s-α1+s-α2+s-(α1+α2).Similarly, for the Laplace transform of the solution (5.4b),
(7.10)x̂#1(s)=[1+s-α2]s-α1c#11-s-(α1+α2)c#211-(K-2)s-α1+s-α2+s-(α1+α2),x̂#2(s)=[1-(K-2)s-α1]s-α2c#21+(K-1)s-(α1+α2)c#111-(K-2)s-α1+s-α2+s-(α1+α2).
In the following subsections, we study the conditions under which asymptotically stable oscillations are possible for a fractional Wien bridge oscillator and also present numerical simulations of the capacitor voltages.


## 7.3. Equal-Order Fractional Wien Bridge

The classical Wien bridge oscillator produces a stable sinusoidal output when its amplifier gainK=3. The amplitude and frequency of the sinusoid are a function of the initial capacitor voltages and the circuit time constant. For a fractional capacitor, the current-voltage relationship is dependent on both capacitor value and order; hence an additional degree of freedom is introduced into the Wien bridge circuit. We consider first the simple case with equal real fractional orders α1=α2=α. Equation (7.7) then simplifies to
(7.11)x*1=c*10ε̃α,α,2α,1(K-2,-1,-1:t)+(c*10-c*20)ε̃α,α,2α,α(K-2,-1,-1:t),x*2=c*20ε̃α,α,2α,1(K-2,-1,-1:t)+[c*10(K-1)-c*20(K-2)]ε̃α,α,2α,α(K-2,-1,-1:t).
To gain further insight into the system’s behaviour, it is advantageous to express the solution in a simpler form using only 1-parameter Mittag-Leffler functions Eα,1(λitα). From (7.9a),
(7.12)x̂*1(s)=sα-1c*10sα+(c*10-c*20)s2α-(K-3)sα+1=σ1sα-1sα-λ1+σ2sα-1sα-λ2,
where the σi∈ℂ are constants to be determined from partial fraction decomposition (see equivalent method in Section 6). The inverse transform yields
(7.13)x*1(t)=σ1Eα,1(λ1tα)+σ2Eα,1(λ2tα).
The solution for x*2(t) can be found in a similar manner. For brevity, we present only numerical simulations of the Caputo solutions (one plot of the Riemann-Liouville solution is presented for comparison). In order for the Wien bridge to produce sustained oscillations, we need to impose condition (6.18) on the complex roots λi. For the current system, this translates to the following expression for K:
(7.14)K=3+2cos⁡(απ2).
Hence, the amplifier gain K is no longer a constant as in the case of the classical Wien bridge but a function of capacitor order α. Simulations of (7.13) were plotted using Mathematica.Figure3 shows plots of x*1 and x*2 for α=0.3,0.5,0.7, and 1.0. There is a clear dependence of waveform amplitude on the fractional-order. The plot for α=1 corresponds to the classical Wien bridge (with ordinary capacitors) and is included for comparison. It is important to keep in mind that the time axes are normalized and oscillation frequency ωα actually varies with order as
(7.15)ωα=(RC)-1/α.
This can be seen in Figure 4(a). It is interesting to note that the values of resistance and capacitance have no effect on the output waveform amplitude (Figure 4(b)). Hence, the frequency of oscillation can be controlled by both the value C and order α of the capacitors. As noted in [31], a clear advantage of this is that high frequencies can be obtained by reducing the order of the capacitors rather than their value, which can remain sufficiently large.Caputo model of fractional Wien bridge with equal orders—comparison of phase and amplitude for capacitor voltages. Parameters:x*1(0)=x*2(0)=0.03. Legend: x*1 (solid line); x*2 (dashed line).
(a)(b)(c)(d)Variation of waveform characteristics forx1(t). Parameters: x*1(0)=x*2(0)=0.05, (a) Capacitor order affects both frequency and amplitude. Legend: α=0.3 (solid), 0.5 (dashed), 0.7 (dash-dotted), 1.0 (dotted) and (b) Time constant affects frequency while amplitude remains constant. Legend: RC=0.8 (solid), 0.9 (dashed), and 1.0 (dash-dotted), and 1.1 (dotted).
(a)(b)In Figure5 we see that the Riemann-Liouville solutions (7.8) are very similar to the Caputo solutions in terms of frequency and amplitude but differ in phase due to the second parameter β of the Mittag-Leffler function. Of particular concern is the fact that the former tends to diverge at the origin since the fractional initial conditions (2.5b) do not correspond to measurable physical quantities. This is an important distinction between the two definitions and has to be taken into consideration when modeling physical systems.Figure 5
Comparison ofx1(t) waveform for both derivative models. Solution for Riemann-Liouville model diverges as t→0. Parameters: x*1(0)=x*2(0)=0.03. Legend: Caputo (solid line), and Riemann-Liouville (dashed line).


## 7.4. Multiorder Fractional Wien Bridge

The Wien bridge circuit can be further generalized by allowing the orders to assume valuesαj=μjα as detailed in Section 6. We use the case α2=2α1 as a starting point. Substituting α1=α and α2=2α into (7.7), (7.9a) and (7.9b), we have
(7.16)x*1=c*10[ε̃α,2α,3α,1(K-2,-1,-1:t)+ε̃α,2α,3α,2α(K-2,-1,-1:t)]-c*20ε̃α,2α,3α,α(K-2,-1,-1:t),x*2=c*20[ε̃α,2α,3α,1(K-2,-1,-1:t)-(K-2)ε̃α,2α,3α,α(K-2,-1,-1:t)]+(K-1)c*10ε̃α,2α,3α,2α(K-2,-1,-1:t).
As with the equal-order bridge, we express the solution in Laplace domain and use partial fractions to obtain a more tractable form:
(7.17)x̂*1(s)=sα-1(s2α+1)c*10-sαc*20s3α-(K-2)s2α+sα+1=∑k=13σksα-1sα-λk,(7.18)x*1(t)=∑k=13σkEα,1(λktα).Only solutions with one negative real and two complex-conjugate roots will be of concern to our present discussion. To justify this, we note that the alternative case of three real roots is of no physical interest as it does not produce oscillatory solutions. With the exception of the exponentially decaying term (due to the negative real root, see asymptotic analysis in Section6.2), the solution is hence similar to the case with equal capacitor orders; that is, the output of the fractional Wien bridge can be expressed as a linear combination of Mittag-Leffler functions. Unfortunately, the relationship between K and α is not as simple as in the equal-order case (7.14). Although K is still a function of α, its form is sufficiently complex that a more convenient alternative is to define K implicitly, that is, find α=ϕ(K), and use polynomial curve-fitting as shown in Figure 6.Determination of amplifier gainK: (a) restrictions on possible values of K and α; (b) third-order least-squares approximation of K(α) for x∈(0,0.5).
(a)(b)Using the method of least-squares, we obtain a third-order approximation forK:
(7.19)K≈4.611-4.821α2+0.008α3.Two restrictions apply to the usable range ofK and α. The first is the requirement that the real root be negative. Plotting the denominator of (7.17) as a function of sα for various amplifier gains, we obtain the relationship in Figure 7. For 0<K<K0≈4.611, we have λ1∈ℝ- and λ2=λ̅3∈ℂ as required. Therefore, this serves as an upper limit to the amplifier gain. The lower limit can be determined by recalling that 0<α2=2α1<1 so that 0<α<0.5. The result of these restrictions is also shown in Figure 6.Figure 7
Properties of roots of (7.19) (shown as horizontal intercepts) for parameter K.Within the stipulated range, the values of gain calculated from the polynomial curve (7.19) are sufficiently accurate to create oscillatory solutions, as demonstrated in the following simulations (Figure 8).Caputo model of fractional Wien bridge withα2=2α1. The effect of the nonoscillatory term in (7.20) can be observed as an initial offset that decays asymptotically as t→∞. Parameters: x*1(0)=x*2(0)=0.03, (a) Comparison of phase and amplitude for capacitor voltages. (b) The limiting case of one ordinary capacitor and one semicapacitor (order 1/2). Legend: and x*1 (solid line), x*2 (dashed line).
(a)(b)As mentioned earlier, it is possible to extend this procedure to any case where one order is an integer multiple of the other. Consider the solution ofx*1 for α1=α, α2=υα:
(7.20)x̂*1(s)=sα-1(sυα+1)c*10-s(υ-1)αc*20s(υ+1)α-(K-2)sυα+sα+1=∑k=1υ+1σksα-1sα-λk,(7.21)x*1(t)=∑k=1μ+1σkEα,1(λktα).
Indeed, unless we are concerned about obtaining the actual Laplace solution in partial fraction form, the value of K that produces asymptotic oscillations can be found by simply studying the roots of the denominator in (7.20) (a polynomial in sα of order υ+1) and imposing suitable conditions as previously shown so that at least one pair of Mittag-Leffler terms has an eigenvalue that satisfies (6.18). For example, when υ=3, a possible solution contains 4 roots in 2 complex-conjugate pairs. One can adjust the amplifier gain such that the roots satisfy |θ|=α0π/2 for the first pair and |θ|>α0π/2 for the second pair, hence resulting in sustained oscillation and asymptotically decaying oscillation, respectively.


## 8. Concluding Remarks

We have proposed a new direct operational method for solving coupled linear fractional differential equations of multiorders. This technique provides an alternative way for solving some linear CFDEs, and the solutions so obtained can be expressed in terms of multivariate Mittag-Leffler functions. For the special cases where each of the multiorders is an integer multiple of a real positive number, the solutions can be further reduced to linear combinations of Mittag-Leffler functions of a single variable. Conditions for asymptotically oscillatory solutions are considered. Two examples, namely, the coupled fractional harmonic oscillator and the fractional Wien bridge circuit, are given to illustrate our method. Simulations of solutions and stability conditions are given. Note that to obtain the solution based on our method requires the use of the Laplace transform of the multivariate Mittag-Leffler function, which then gives the identity for getting the Laplace inversion for the solution. This is one main advantage of the direct operational inversion method proposed here as it will give the solution directly. We remark that our method does not actually simplify the computational aspect of obtaining solutions, though intuitively it allows one to obtain the solution in explicit form.Here we would like to remark that there were attempts recently to transform CFDEs with different multiorders into an equivalent system of CFDEs of a single order [32, 33]. Such a method again does not reduce the amount of computation necessary to obtain the solutions; instead, due to the increase in the number of the auxiliary equations in the latter system, it is actually more tedious to obtain the full solutions. Our view on CFDEs is that, in general, one still has to use numerical methods to obtain approximate solutions. The point is to find a method that provides a more efficient way of doing so. We hope to look into this aspect in a future work. Finally, it will be interesting to consider whether the above method can be extended to nonlinear CFDEs. One expects that such a generalization will not be straightforward.



---

*Source: 653939-2011-12-29.xml*