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

#### 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]:where the fractional integral is defined for as When referring to either definition, we simply use the notation .

Let us consider a linear-coupled system of inhomogeneous fractional differential equations of the form where and are vectors of dimension , , is an -matrix, and is the fractional differential operator given by the diagonal matrix operator: The orders are real positive numbers with , where is a positive integer for each . The boundary conditions for (2.3) are given byand . We remark that the ’s generally differ in value. Let its maximum and minimum values be denoted by and , 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 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 The solution of (2.3) can then be expressed as Unfortunately, the inverse of may not exist. However, the right-inverse does exist for both the Riemann-Liouville and Caputo cases with The main task now is to determine the right-inverse of 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 derivative with arbitrary and , and similarly for the Riemann-Liouville derivative, Applying the fractional integral operator to (2.3) and using the initial conditions (2.5a) and (2.5b) gives where is given by Let Now by rearranging (2.10) one gets The operator has an inverse . One possible representation of is given by This form may not be a simple one, since in most cases and do not commute. However, the verification is straightforward: The other possible representation is where is the adjoint of , that is, , and the is the inverse operator of such that . It is quite simple to verify that this representation is the inverse of , and this will be shown in Section 4.

Now we define One can easily verify that is the right inverse of . Thus, the solution is given by 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 with The initial conditions have to be transformed accordingly and they become

In the Riemann-Liouville case, we cannot absorb the term into as in the Caputo case. However, if one compares the initial condition terms and the term , one can see that if we modify with , 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 can be transformed into a homogeneous linear fractional differential equation. The solution of the transformed equation can then be written as When is time dependent with power up to , the above modification can still apply to the Caputo case. In the Riemann-Liouville case, if are analytic functions, then the same modification as above holds if (2.22) is altered with the summation from to .

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 , then , and It would be convenient if we introduce the Mittag-Leffler function with matrix argument [27–29]: Then Here we have used the following definition of the Dirac delta function: The matrix can always be decomposed into Jordan normal form. However, we consider only the case where it can be diagonalized: with eigenvalues as the diagonal elements of , and 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 where and are nonnegative real numbers. The initial conditions are Thus For simplicity, we use the following notation: We first consider the simpler case with in this section. Since is symmetric, it can be diagonalized to give two eigenvalues: Both eigenvalues are nonpositive and Just as we have shown in (3.7), the solutions are again linear combinations of Mittag-Leffler functions:Figure 1 shows simulations of the Caputo solution (3.14a) for orders . 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 () and undamped oscillation ().

**(a)**

**(b)**

**(c)**

**(d)**

#### 4. The Adjoint Method

As we mentioned in Section 2.1, (2.16), the main task now is to calculate the inverse of . Note that (2.16) can be reexpressed as so that Also recall that 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 The inverse is given by Its kernel is the multivariate Mittag-Leffler function of the kind given by (B.3), Appendix B: The adjoint of 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 are given by

##### 4.2. The Solutions

The solutions based on the adjoint method are given by

#### 5. Laplace Transform Method

In this section we briefly discuss how the widely used Laplace transform technique can be employed to determine Green’s function and the operator . 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 using (4.5): and the Laplace transform of the adjoint kernel : Next, the Laplace transforms of the part related to the initial conditions from (2.11) are Thus the Laplace transforms of the solutions become and can be obtained just by interchanging . 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.

#### 6. Multiple Fractional-Order System

In physics and engineering problems the fractional-orders can often be approximated by rational numbers, that is, , for some . Thus one gets , where is the least common multiple of with some . However, we can also consider the more general case with , for some . Here, 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 , the expansion of the determinant will be a polynomial of order in . By the fundamental theorem of algebra, it must have in general complex roots, that is, for . Note that since all coefficients of the polynomial are real, if any root is complex, its complex conjugate is also a root. That means that the complex roots occur in pairs. For convenience we write , for . We can then factorize the polynomial If all roots are distinct, the inverse can be written as a partial fraction:

##### 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

To simplify the problem we consider the case where , that is, . The extension to the general case as above is straightforward: Using the following formula: (6.4) can be written as where

##### 6.2. Solutions with Asymptotic Oscillations

It is clear from the previous expansion of the Caputo terms that for the solution , as , any possible oscillation arises from the pair of complex roots , while all (negative) real roots must result in an asymptotic decay. Note that each term in (6.9a) with will grow asymptotically as the power law . However, when we combine the terms and reexpress the equation explicitly as
for the particular case with , one has
The second term is the only term that grows asymptotically as , 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 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 .

Similar consideration can be given to Riemann-Liouville system. However, now approaches a constant or possible zero as if , and . Thus the possible asymptotically stable oscillations are due to the term with its corresponding complex conjugate. If one looks at (6.8b), one can write explicitly the asymptotic expansion: One has to impose the condition that be a purely imaginary number. Let us denote the inversion of the root by ; 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 th term of (6.7b) as and we then have the asymptotic solution: can be evaluated in a similar way; it has the same period but different modulus and phase: 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 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 (Figure 2). 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.)

**(a)**

**(b)**

It is well known that a fractional differential equation of order 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 where is the amplifier saturation voltage and time axes (normalized with respect to time constant ). Using basic circuit analysis, it can be shown that the capacitor voltages are related via a 2-dimensional CFDE: where Here is the amplifier gain (i.e., ). In the linear region of the amplifier, 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.1. Solution Using the Adjoint Method

Substituting (7.4) into the solution (4.10a) we get for the Caputo case Similarly, for the Riemann-Liouville case, substituting (7.4) into the solution (4.10b), we get

##### 7.2. Solution Using the Laplace Transform

Substituting (7.4) into the Laplace transform solution (5.4a), we obtainSimilarly, for the Laplace transform of the solution (5.4b), 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 gain . 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 . Equation (7.7) then simplifies to 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 . From (7.9a), where the are constants to be determined from partial fraction decomposition (see equivalent method in Section 6). The inverse transform yields The solution for 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 . For the current system, this translates to the following expression for : Hence, the amplifier gain 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 and for , and 1.0. There is a clear dependence of waveform amplitude on the fractional-order. The plot for 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 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 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.

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

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.

##### 7.4. Multiorder Fractional Wien Bridge

The Wien bridge circuit can be further generalized by allowing the orders to assume values as detailed in Section 6. We use the case as a starting point. Substituting and into (7.7), (7.9a) and (7.9b), we have 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 and is not as simple as in the equal-order case (7.14). Although is still a function of , its form is sufficiently complex that a more convenient alternative is to define implicitly, that is, find , and use polynomial curve-fitting as shown in Figure 6.

**(a)**

**(b)**

Using the method of least-squares, we obtain a third-order approximation for :

Two restrictions apply to the usable range of and . The first is the requirement that the real root be negative. Plotting the denominator of (7.17) as a function of for various amplifier gains, we obtain the relationship in Figure 7. For , we have and as required. Therefore, this serves as an upper limit to the amplifier gain. The lower limit can be determined by recalling that so that . 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).

**(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 of for , : Indeed, unless we are concerned about obtaining the actual Laplace solution in partial fraction form, the value of that produces asymptotic oscillations can be found by simply studying the roots of the denominator in (7.20) (a polynomial in of order ) 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 , a possible solution contains 4 roots in 2 complex-conjugate pairs. One can adjust the amplifier gain such that the roots satisfy for the first pair and 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.

#### Appendices

#### A. Mittag-Leffler Function and Related Functions

The Mittag-Leffler function [26, 28] and its generalizations are defined as follows: where

Note that Thus we have For convenience we define the following functions:

##### A.1. Asymptotic Expansion of Mittag-Leffler Function [8, 24, 26]

For , Similarly one has

#### B. Multivariate Mittag-Leffler Functions [27–29]

Let us adopt the following notations: and the binomial coefficient is generalized to 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 with . The transform of (4.10a) and (4.10b) is then which agrees precisely with (5.4a) and (5.4b).

#### Acknowledgments

The authors would like to thank Professor W. Ahmad for correspondence. S. C. Lim would like to thank the Malaysian Ministry of Science, Technology and Innovation for the support under its Brain Gain Malaysia (Back to Lab) Program. Li acknowledges the 973 plan under the project no. 2011CB302802, and the NSFC under the project Grant nos. 60873264, 61070214. S. Y. Chen acknowledges the NSFC under the project Grant no. 60870002 and Zhejiang Provincial Natural Science Foundation (R1110679).