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

1 Faculty of Engineering, Multimedia University, Selangor Darul Ehsan, 63100 Cyberjaya, Malaysia 2 Department of Chemistry, Faculty of Science, Chulalongkorn University, Bangkok 10330, Thailand 3 School of Information Science & Technology, East China Normal University, No. 500, Dong-Chuan Road, Shanghai 200241, China 4 College of Computer Science, Zhejiang University of Technology, Hangzhou 310023, China


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

Linear-Coupled Fractional Differential Equations
We consider two types of fractional derivatives 22-26 : where the fractional integral is defined for γ > 0 as 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 where X x 1 , . . ., x n and B b 1 , . . ., b n are vectors of dimension n, A a ij , i, j 1, . . ., n is an n×n -matrix, and D α is the fractional differential operator given by the diagonal matrix operator: The solution of 2.3 can then be expressed as Unfortunately, the inverse of L may not exist.However, the right-inverse G does exist for both the Riemann-Liouville and Caputo cases with LG 1 / GL.

2.8
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.

The Right-Inverse Operator
For the Caputo derivative D α * with arbitrary α and m, and similarly for the Riemann-Liouville derivative,

2.9b
Applying the fractional integral operator I α to 2.3 and using the initial conditions 2.5a and 2.5b gives where W is given by

2.11
Let Now by rearranging 2.10 one gets QX I α B W.

2.13
The operator Q has an inverse K.One possible representation of K is given by

2.14
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
The other possible representation is where Q * is the adjoint of Q, that is, Q * Q det Q, and the Ψ is the inverse operator of det Q such that Ψ det Qf 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 G KI α .

2.17
One can easily verify that G is the right inverse of L. Thus, the solution is given by X GB KW.

2.18
The detailed evaluation can be carried out by using different techniques, a few of which will be considered here.

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 with

2.20
The initial conditions have to be transformed accordingly and they become

2.21
In the Riemann-Liouville case, we cannot absorb the term B 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 with c 0 #i b i , 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 term b i can be transformed into a homogeneous linear fractional differential equation.The solution of the transformed equation can then be written as X KW.

2.23
When b i is time dependent with power up to m i , the above modification can still apply to the Caputo case.In the Riemann-Liouville case, if b i • are analytic functions, then the same modification as above holds if 2.22 is altered with the summation from −∞ to m i − 1.
In subsequent sections, we will consider the solution of 2.3 according to the above modifications.

System with Equal Fractional Orders
In the case where all α i α, then I α I α 1, and 3.1 It would be convenient if we introduce the Mittag-Leffler function with matrix argument 27-29 :

3.3
Here we have used the following definition of the Dirac delta function:

3.4
The matrix A can always be decomposed into Jordan normal form.However, we consider only the case where it can be diagonalized: with eigenvalues λ j as the diagonal elements of Λ, and

3.7
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.

Coupled Oscillator with Equal Fractional Orders
Here we demonstrate our method by considering a linear-coupled oscillator system given by where ω 2 and ω 2 are nonnegative real numbers.The initial conditions are

3.10
For simplicity, we use the following notation:

3.11
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
Both eigenvalues are nonpositive and Just as we have shown in 3.7 , the solutions are again linear combinations of Mittag-Leffler functions: x * 1 c 0 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 .

The Adjoint Method
As we mentioned in Section 2.1, 2.16 , the main task now is to calculate the inverse of det Q.Note that 2.16 can be reexpressed as 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 twodimensional system.

Two-Dimensional System
In a two-dimensional system the determinant is easy to calculate, and we have The inverse is given by Its kernel is the multivariate Mittag-Leffler function of the kind given by B.3 , Appendix B:

4.7
The adjoint of Q is 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

The Solutions
The solutions based on the adjoint method are given by x * 1 4.10b

Laplace Transform Method
In this section we briefly discuss how the widely used Laplace transform technique can be employed to determine Green's function G 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 Section 2.2 the system of CFDEs is assumed to be homogeneous.We begin by calculating the Laplace transform of det Qδ t using 4.5 : and the Laplace transform of the adjoint kernel

5.2
Next, the Laplace transforms of the part related to the initial conditions from 2.11 are

5.3
Thus the Laplace transforms of the solutions become 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 Appendix C for a 2-dimensional system.

Multiple Fractional-Order System
In physics and engineering problems the fractional-orders can often be approximated by rational numbers, that is, α i p i /q i , for some p i , q i ∈ AE.Thus one gets α i μ i /q, where q is the least common multiple of q 1 , q 2 , . . . ,q n with some μ i ∈ AE.However, we can also consider the more general case with α j μ j α 0 , for some μ j ∈ AE.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 If all roots are distinct, the inverse can be written as a partial fraction: 6.2

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
To simplify the problem we consider the case where 0 < max μ 1 α 0 , μ 2 α 0 ≤ 1, that is, m 1 m 2 1.The extension to the general case as above is straightforward:

6.4
Using the following formula: 6.4 can be written as where

Solutions with Asymptotic Oscillations
It is clear from the previous expansion of the Caputo terms that for the solution x * 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 t pα 0 .However, when we combine the terms and reexpress the equation explicitly as for the particular case with μ 1 1, μ 2 2, one has 11

6.12
The second term is the only term that grows asymptotically as ∼ t α 0 , which does not contribute since is zero.The verification of this result is given here for the general case of n-roots.Let us consider general partial fractions: We have . . . . . .
We can rewrite 6.14.n−1 as Using 6.15 with n 3 in 6.12 , we have 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, now g j #1 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 j #1 with its corresponding complex conjugate.If one looks at 6.8b , one can write explicitly the asymptotic expansion: 6.17 One has to impose the condition that λ 1/α 0 j be a purely imaginary number.Let us denote the inversion of the root by λ j |λ j |e iθ ; we must have 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 r j e iϕ , 6.19 and we then have the asymptotic solution: x 1 −→ r j e i|λ j |t iϕ r j e −i|λ j |t−iϕ 2r j cos λ j t ϕ .

6.20a
x 2 can be evaluated in a similar way; it has the same period but different modulus and phase: x 2 −→ 2r j cos λ j t ϕ .

6.20b
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.

Wien Bridge System
In this section we apply the solution methods discussed earlier to model a fractional-order Wien bridge oscillator Figure 2 .The Wien bridge is a common electronic circuit that can generate a sinusoidal output signal without requiring an oscillatory input.The resistorcapacitor 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.
It is well known that a fractional differential equation of order 0 < α < 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 voltages x i v i /V sat where V sat 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 2dimensional CFDE: where

7.2b
Here K is the amplifier gain i.e., v o Kv 1 .In the linear region of the amplifier, −1 < Kx 1 < 1 and 7.2a simplifies to Thus we have to solve the homogeneous linear fractional-order differential system with For fractional capacitors, the real orders are restricted to 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

Solution Using the Adjoint Method
Substituting 7.4 into the solution 4.10a we get for the Caputo case

7.7
Similarly, for the Riemann-Liouville case, substituting 7.4 into the solution 4.10b , we get 7.8

Solution Using the Laplace Transform
Substituting 7.4 into the Laplace transform solution 5.4a , we obtain Similarly, for the Laplace transform of the solution 5.4b ,

7.10
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.

Equal-Order Fractional Wien Bridge
The classical Wien bridge oscillator produces a stable sinusoidal output when its amplifier gain K 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
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 λ i t α .From 7.9a , where the σ i ∈ are constants to be determined from partial fraction decomposition see equivalent method in Section 6 .The inverse transform yields 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: 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.Figure 3 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 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 .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.
In Figure 5 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.

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
As with the equal-order bridge, we express the solution in Laplace domain and use partial fractions to obtain a more tractable form: 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 Section 6.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.
Using the method of least-squares, we obtain a third-order approximation for K: K ≈ 4.611 − 4.821α 2 0.008α 3 .

7.19
Two restrictions apply to the usable range of K 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 < K 0 ≈ 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.
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 .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 of x * 1 for α 1 α, α 2 υα: 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.

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.

A. Mittag-Leffler Function and Related Functions
The Mittag-Leffler function 26, 28 and its generalizations are defined as follows: where Note that γ 0 1, 0 n 0 for n 0, 0 0 1.A.3 Thus we have For convenience we define the following functions:

B. Multivariate Mittag-Leffler Functions [27-29]
Let us adopt the following notations: The multivariate Mittag-Leffler functions is defined as:

C. Equivalence of Adjoint Method and Laplace Transform Solutions
We demonstrate the equivalence of the time-domain and frequency-domain Laplace solutions for the 2-dimensional system as given by 4.10a and 4.10b and 5.4a and 5.4b , respectively.The generalization to systems of higher dimension is straightforward and will be omitted for brevity.The Laplace transform of the multivariate Mittag-Leffler function in A.5a is L ε α 1 ,...,α n ,β a 1 , . . ., a n : t s −β 1 − n i 1 a i s −α i , C.1 with β > 0. The transform of 4.10a and 4.10b is then which agrees precisely with 5.4a and 5.4b .

Figure 2 :
Figure 2: Wien bridge oscillator: a circuit schematic with operational amplifier and b simplified circuit diagram for voltage analysis.

Figure 6 :Figure 7 :
Figure 6: Determination of amplifier gain K: a restrictions on possible values of K and α; b third-order least-squares approximation of K α for x ∈ 0, 0.5 .

α 1 = 5 bFigure 8 :
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 .