# Approximate Analytical Solution for Nonlinear System of Fractional Differential Equations by BPs Operational Matrices

Authors: Mohsen Alipour; Dumitru Baleanu
Journal: Advances in Mathematical Physics
 (2013)

---

## Abstract

We present two methods for solving a nonlinear system of fractional differential equations within Caputo derivative. Firstly, we derive operational matrices for Caputo fractional derivative and for Riemann-Liouville fractional integral by using the Bernstein polynomials (BPs). In the first method, we use the operational matrix of Caputo fractional derivative (OMCFD), and in the second one, we apply the operational matrix of Riemann-Liouville fractional integral (OMRLFI). The obtained results are in good agreement with each other as well as with the analytical solutions. We show that the solutions approach to classical solutions as the order of the fractional derivatives approaches 1.


---

## Body

## 1. Introduction

Differential equations of fractional order have been subjected to many studies due to their frequent appearance in various applications in fluid mechanics, viscoelasticity, biology, physics, engineering, and so on. Recently, a large amount of literature was developed regarding the application of fractional differential equations in nonlinear dynamics (see, e.g., [1–11] and the references therein). Thus, a huge attention has been given to the solution of fractional ordinary differential equations, integral equations, and fractional partial differential equations of physical interest. As it is known, there exists no method that yields an exact solution for fractional differential equations. Various methods have been proposed in order to solve the fractional differential equations. These methods include the homotopy perturbation method [12–15], Adomian’s decomposition method [16–20], variation iteration method  [12–14, 21–23], homotopy analysis method [24], differential transform method [25], operational matrices [26–28], and nonstandard finite difference scheme [29].In this paper, we investigate the nonlinear system of fractional differential equations as(1)Dαixi(t)=gi(t,X(t)),i=1,…,n,0<t≤1,0<αi≤1,
and the initial condition
(2)X(0)=X0,
where X(t)=[x1(t),…,xn(t)]T and X0=[x0,1,…,x0,n]T. Also, gi:[0,1]×Rn→R are multivariable polynomial functions.The structure of the paper is given later. In Section2, we present some preliminaries and properties in fractional calculus and Bernstein polynomials. In Section 3, we make operational matrices for product, power, Caputo fractional derivative, and Riemann-Liouville fractional integral by BPs. In Section 4, we apply two methods for solving nonlinear system of fractional differential equations by BPs. In Section 5, numerical examples are simulated to demonstrate the high performance of the proposed method.  Conclusions are presented in Section 6.


## 2. Basic Tools

In this section, we recall some basic definitions and properties of the fractional calculus and Bernstein polynomials.Definition 1 (see [2, 7, 10]).
The Riemann-Liouville fractional integral operator of orderα≥0, of a function f∈Cμ,μ≥-1, is defined as
(3)Iαf(t)=1Γ(α)∫0t(t-x)α-1f(x)dx,α>0,t>0,I0f(t)=f(t),
and for n-1<α≤n,  n∈N,  t>0, f∈C-1n, the fractional derivative of f(t) in the Caputo sense is defined as
(4)Dαf(t)=1Γ(n-α)∫0t(t-x)n-α-1f(n)(x)dx,
where for n∈N,μ∈R we have
(5)Cμ={f(t)∣f(t)>0fort>0,f(t)=tpf1(t)wherep>μ,f1(t)∈C[0,∞)},Cμn={f(t)∣f(n)(t)∈Cμ}.
Also, if n-1<α≤n,n∈N, and f∈Cμn,μ≥-1, then
(6)(1)DαIαf(t)=f(t),(7)(2)IαDαf(t)=f(t)-∑k=0n-1f(k)(0+)xkk!,t>0.Definition 2 (see [30]).
The Bernstein polynomials (BPs) ofmth degree are defined on the interval [0,1] as follows:
(8)Bi,m(x)=(mi)xi(1-x)m-i,i=0,1,…,m.Lemma 3.
One can writeΦm(x)=ATm(x), where A is a matrix upper triangular,  Tm(x)=[1,x,…,xm]T, and Φm(x)=[B0,m(x),B1,m(x),…,Bm,m(x)]T.Proof.
(see  [26]).Definition 4.
We can define the dual matrixQ(m+1)×(m+1) on the basis of Bernstein polynomials of mth degree as follows:
(9)Q=∫01Φ(x)Φ(x)Tdx,
where
(10)(Q)i+1,j+1=∫01Bi,m(x)Bj,m(x)dx=(mi)(mj)(2m+1)(2mi+j),i,j=0,1,…,m.Lemma 5.
LetL2[0,1] be a Hilbert space with the inner product 〈f,g〉=∫01f(x)g(x)dx and y∈L2[0,1]. Then, we can find the unique vector c=[c0,c1,…,cm]T such that cTΦm(x)is the best approximation of y(x) from space Sm=Span{B0,m(x),B1,m(x),…,Bm,m(x)}. Moreover, one can get c=Q-1〈y,Φm〉, such that 〈y,Φm〉=[〈y,B0,m〉,〈y,B1,m〉,…,〈y,Bm,m〉]T.Proof .
(see [31]).Lemma 6.
Suppose that the functiony:[0,1]→R is m+1 times continuously differentiable (i.e.,y∈Cm+1([0,1])). If cTB is the best approximation y out of Sm, then
(11)∥y-cTB∥L2[0,1]≤K^(m+1)!2m+3,
where K^=maxx∈[0,1]|y(m+1)(x)|. Also, if y∈C∞([0,1]), then the error bound vanishes.Proof .
(see [32]).


## 3. Operational Matrices of Bernstein Polynomials

In Section3, we recall the operational matrices for product, power, Caputo fractional derivative and Riemann-Liouville fractional integral by BPs.Lemma 7.
Suppose thatc(m+1)×1 is an arbitrary vector. The operational matrix of product  C^(m+1)×(m+1) using BPs can be given as follows:
(12)cTΦm(x)Φm(x)T≈Φm(x)TC^.Proof .
(see [27]).Corollary 8.
Suppose thaty(t)≈cTΦm(t), x(t)≈dTΦm(t), and C^(m+1)×(m+1) is the operational matrix of product using BPs for vector c. One can get the approximate function for x(t)y(t) using BPs as follows:
(13)y(t)x(t)≈Φm(t)TC^d.Proof.
By using Lemma7, it is clear.Corollary 9.
Suppose thaty(t)≈cTΦm(t) and C^(m+1)×(m+1) is the operational matrix of product using BPs for vector c. One can get the approximate function for yk(t)(k∈N), using BPs as follows:
(14)yk(t)≈Φm(t)TC~k,
where C~k=C^k-1c.Proof .
(see [26]).Theorem 10.
One can get BPs operational matrixDα from order (m+1)×(m+1) for the Caputo fractional derivative as follows:
(15)DαΦm(t)=1Γ(n-α)∫0t(t-x)n-α-1Φm(n)(x)dx≈DαΦm(t).Proof.
See [26] for details.Theorem 11.
One can obtain the operational matrixFαfrom order (m+1)×(m+1) for the Riemann-Liouville fractional integral on the basis of BPs from order m as
(16)IαΦm(t)=1Γ(α)∫0t(t-x)α-1Φm(x)dx≈FαΦm(t).Proof.
See [28] for details.


## 4. Solving System of Fractional Differential Equations

In this section, we use two methods for solving system of fractional differential equations. In the first method, we use the operational matrix for Caputo fractional derivative (OMCFD), and in the second method, we apply the operational matrix for Riemann-Liouville fractional integral (OMRLFI).


### 4.1. Solving the Problem by OMCFD

Using Lemma5, we can approximate the functions xi(t) as follows:
(17)xi(t)≈CiTΦm(t),i=1,…,n,
where Ci∈R(m+1)×1.From (17) and (15)  we can write
(18)Dαixi(t)≈CiTDαiΦm(t),i=0,…,n.
Therefore, problem (1) and (2) reduces to the following problem:
(19)CiTDαiΦm(t)=gi(t,C1TΦm(t),…,CnTΦm(t)),i=1,…,n,
and the initial condition
(20)CiTΦm(0)=x0,i,i=1,…,n.
Now, using Lemma 5 we can approximate all of the known functions in the system (19). Then, by using Lemma 7 and Corollaries 8 and 9, since functions giare polynomial, we obtain the following approximations:
(21)gi(t,X(t))≈Gi(C1,…,Cn)Φm(t),i=1,…,n,
where Gi:R(m+1)×n→R1×(m+1).Also, for eachi(i=1,…,n), by using tau method [33] we can generate algebraic equations from (19) and (21) as follows
(22)G~i,j=∫01(CiTDαi-Gi(C1,…,Cn))×Φm(t)Bj,m(t)dt=0,j=0,…,m-1,
and from (23) we set G~i,m=CiTΦm(0)-x0,i.Finally, problem (1) and (2) has been reduced to the system of algebraic equations
(23)G~i,j(C1,…,Cn)=0,i=1,…,n,j=0,…,m.
The aforementioned system can be solved for Ci by Newton’s iterative method. Then, we get the approximate value of the functions xi(t) from (17).


### 4.2. Solving the Problem by OMRLFI

This method consists of two steps.Step 1.
Initial conditions are used to reduce a given initial-value problem to a      problem with zero initial conditions. Therefore we have a modified system, incorporating the initial values.Step 2.
The BPs operational matrix of Riemann-Liouville fractional integral is used to transform the problem into a system of algebraic equations.Now, from (2) we define
(24)xi(t)=x0,i+zi(t),i=1,2,…,n,
where zi(t),i=1,2,…,n, are the new unknown functions.Substituting (24) in (1) and (2), we have the following system:
(25)Dαizi(t)=fi(t,Z(t)),i=1,…,n,0<t≤1,0<αi≤1,
and the initial condition
(26)zi(0)=0,i=1,…,n,
where Z(t)=[z1(t),…,zn(t)]T and fi:[0,1]×Rn→R are multivariable polynomial functions. We use the following approximation:
(27)Dαizi(t)≈C~iTΦm(t),i=1,…,n,
where C~i∈R(m+1)×1 are unknown vectors. From (7), (27), and Theorem 11, we can write
(28)zi(t)=IαiDαizi(t)≈Iαi(C~iTΦm(t))=C~iTIαiΦm(t)≈C~iTFαiΦm(t).
So, by (27) and (28), problem (25) and (26) reduces to the following problem:
(29)C~iTΦm(t)=fi(t,C~1TFα1Φm(t),…,C~nTFαnΦm(t)),i=1,…,n.
As we saw in the previous section, we can obtain the following approximations:
(30)fi(t,X(t))≈Fi(C~1,…,C~n)Φm(t),i=1,…,n,
where Fi:R(m+1)×n→R1×(m+1). So, from (29) and (30) we have
(31)(C~iT-Fi(C~1,…,C~n))Φm(t)=0,i=1,…,n.
Therefore, we have reduced problem (1) and (2) to the system of algebraic equations as follows:
(32)C~iT-Fi(C~1,…,C~n)=0,
where this system can be solved for C~i by Newton’s iterative method. Finally we obtain the approximate of the functions xi(t) by
(33)xi(t)≈x0,i+C~iTFαiΦm(t),i=1,2,…,n.


## 4.1. Solving the Problem by OMCFD

Using Lemma5, we can approximate the functions xi(t) as follows:
(17)xi(t)≈CiTΦm(t),i=1,…,n,
where Ci∈R(m+1)×1.From (17) and (15)  we can write
(18)Dαixi(t)≈CiTDαiΦm(t),i=0,…,n.
Therefore, problem (1) and (2) reduces to the following problem:
(19)CiTDαiΦm(t)=gi(t,C1TΦm(t),…,CnTΦm(t)),i=1,…,n,
and the initial condition
(20)CiTΦm(0)=x0,i,i=1,…,n.
Now, using Lemma 5 we can approximate all of the known functions in the system (19). Then, by using Lemma 7 and Corollaries 8 and 9, since functions giare polynomial, we obtain the following approximations:
(21)gi(t,X(t))≈Gi(C1,…,Cn)Φm(t),i=1,…,n,
where Gi:R(m+1)×n→R1×(m+1).Also, for eachi(i=1,…,n), by using tau method [33] we can generate algebraic equations from (19) and (21) as follows
(22)G~i,j=∫01(CiTDαi-Gi(C1,…,Cn))×Φm(t)Bj,m(t)dt=0,j=0,…,m-1,
and from (23) we set G~i,m=CiTΦm(0)-x0,i.Finally, problem (1) and (2) has been reduced to the system of algebraic equations
(23)G~i,j(C1,…,Cn)=0,i=1,…,n,j=0,…,m.
The aforementioned system can be solved for Ci by Newton’s iterative method. Then, we get the approximate value of the functions xi(t) from (17).


## 4.2. Solving the Problem by OMRLFI

This method consists of two steps.Step 1.
Initial conditions are used to reduce a given initial-value problem to a      problem with zero initial conditions. Therefore we have a modified system, incorporating the initial values.Step 2.
The BPs operational matrix of Riemann-Liouville fractional integral is used to transform the problem into a system of algebraic equations.Now, from (2) we define
(24)xi(t)=x0,i+zi(t),i=1,2,…,n,
where zi(t),i=1,2,…,n, are the new unknown functions.Substituting (24) in (1) and (2), we have the following system:
(25)Dαizi(t)=fi(t,Z(t)),i=1,…,n,0<t≤1,0<αi≤1,
and the initial condition
(26)zi(0)=0,i=1,…,n,
where Z(t)=[z1(t),…,zn(t)]T and fi:[0,1]×Rn→R are multivariable polynomial functions. We use the following approximation:
(27)Dαizi(t)≈C~iTΦm(t),i=1,…,n,
where C~i∈R(m+1)×1 are unknown vectors. From (7), (27), and Theorem 11, we can write
(28)zi(t)=IαiDαizi(t)≈Iαi(C~iTΦm(t))=C~iTIαiΦm(t)≈C~iTFαiΦm(t).
So, by (27) and (28), problem (25) and (26) reduces to the following problem:
(29)C~iTΦm(t)=fi(t,C~1TFα1Φm(t),…,C~nTFαnΦm(t)),i=1,…,n.
As we saw in the previous section, we can obtain the following approximations:
(30)fi(t,X(t))≈Fi(C~1,…,C~n)Φm(t),i=1,…,n,
where Fi:R(m+1)×n→R1×(m+1). So, from (29) and (30) we have
(31)(C~iT-Fi(C~1,…,C~n))Φm(t)=0,i=1,…,n.
Therefore, we have reduced problem (1) and (2) to the system of algebraic equations as follows:
(32)C~iT-Fi(C~1,…,C~n)=0,
where this system can be solved for C~i by Newton’s iterative method. Finally we obtain the approximate of the functions xi(t) by
(33)xi(t)≈x0,i+C~iTFαiΦm(t),i=1,2,…,n.


## 5. Examples

To demonstrate the applicability and to validate the numerical scheme, we apply the present method for the following examples.Example 12.
Consider the following linear system of fractional differential equations [24, 25]:
(34)Dα1x1(t)=x1(t)+x2(t),Dα2x2(t)=-x1(t)+x2(t),0<α1,α2≤1,
with initial condition
(35)x1(0)=0,x2(0)=1.
For this problem we have the exact solution in the case of  α1=α2=1 as
(36)x1(t)=etsin(t),x2(t)=etcos(t).
We solved this problem by OMCFD and OMRLFI. Figures 1 and 2 show the approximate solutions of x1(t) and x2(t), respectively, as a function of time for m=10, for different values of α1, α2. The results show that numerical solutions are in good agreement with each other, in both methods. Also, these figures show that as α1, α2 approach close to 1, the numerical solutions approach to the solutions for α1=α2=1 as expected. In Figures 3 and 4, we see the absolute error of both methods, for m=10, α1=α2=1. In these figures, we can see that obtained results using the presented methods agree well with the analytical solutions for α1=α2=1.Figure 1
Approximate solutions ofx1(t) for m=10 and different values of α1, α2 by OMCFD and OMRLFI in Example 12.Figure 2
Approximate solutions ofx2(t) for m=10 and different values of α1, α2 by OMCFD and OMRLFI in Example 12.Plot of absolute error functionx1(t) for α1=α2=1 and m=10 by OMCFD and OMRLFI in Example 12.
(a)
(b)Plot of absolute error functionx2(t) for α1=α2=1 and m=10 by OMCFD and OMRLFI in Example 12.
(a)
(b)Example 13.
Let us consider the following nonlinear fractional system [24] as follows:
(37)Dα1x1(t)=x1(t)2,Dα2x2(t)=x12(t)+x2(t),0<α1,α2≤1,
such that
(38)x1(0)=1,x2(0)=0.The exact solution of this system, whenα1=α2=1, is
(39)x1(t)=et/2,x2(t)=tet.
Figures 5 and 6 show the approximate solutions of x1(t) and x2(t), respectively, for different values of α1, α2 by OMCFD and OMRLFI. We conclude that as α1, α2 approach close to 1, the numerical solutions approach solutions for α1=α2=1 as expected. Furthermore, in both methods, the results agree well with each other. Figures 7 and 8 show that, the absolute error of obtained results for m=10 and α1=α2=1 using OMCFD and OMRLFI is in good agreement with the exact solution.Figure 5
Approximate solutions ofx1(t) for m=10 and different values of α1, α2 by OMCFD and OMRLFI in Example 13.Figure 6
Approximate solutions ofx2(t) for m=10 and different values of α1, α2 by OMCFD and OMRLFI in Example 13.Plot of absolute error functionx1(t) for α1=α2=1 and m=10 by OMCFD and OMRLFI in Example 13.
(a)
(b)Plot of absolute error functionx2(t) for α1=α2=1 and m=10 by OMCFD and OMRLFI in Example 13.
(a)
(b)Example 14.
Consider the nonlinear system of fractional differential equations [24]:
(40)Dα1x1(t)=x1(t),Dα2x2(t)=2x12(t),Dα3x3(t)=3x1(t)x2(t),0<α1,α2,α3≤1,
with the initial conditions given by
(41)x1(0)=1,x2(0)=1,x3(0)=0.
The exact solution of this system, when α1=α2=α3=1, becomes
(42)x1(t)=et,x2(t)=e2t,x3(t)=e3t-1.
We can see the approximate solutions of x1(t),x2(t) and x3(t), by OMCFD and OMRLFI for m=10 and different values of α1, α2 and α3, in Figures 9, 10, and 11. These figures show that, when α1, α2, and α3 approach close to 1, the numerical solutions approach the solutions for α1=α2=α3=1 as expected. In Figures 9–11, we observe that results of OMCFD and OMRLFI overlap. In Figures 12, 13, and 14, we see the absolute error of the obtained results for m=10 and α1=α2=α3=1 in both methods.Figure 9
Approximate solutions ofx1(t) for m=10 and different values of α1, α2, α3 by OMCFD and OMRLFI in Example 14.Figure 10
Approximate solutions ofx2(t) for m=10 and different values of α1, α2, α3 by OMCFD and OMRLFI in Example 14.Figure 11
Approximate solutions ofx3(t) for m=10 and different values of  α1, α2, α3 by OMCFD and OMRLFI in Example 14.Plot of absolute error functionx1(t) for α1=α2=α3=1 and m=10 by OMCFD and OMRLFI in Example 14.
(a)
(b)Plot of absolute error functionx2(t) for α1=α2=α3=1 and m=10 by OMCFD and OMRLFI in Example 14.
(a)
(b)Plot of absolute error functionx3(t) for α1=α2=α3=1 and m=10 by OMCFD and OMRLFI in Example 14.
(a)
(b)


## 6. Conclusion

In this paper, we get operational matrices of the product, Caputo fractional derivative, and Riemann-Liouville fractional integral by Bernstein polynomials. Then by using these matrices, we proposed two methods that reduced the nonlinear systems of fractional differential equations to the two system of algebraic equations that can be solved easily. Finally, numerical examples are simulated to demonstrate the high performance of the proposed method. We saw that the results of both methods were in good agreement with each other, and the classical solutions were recovered when the order of the fractional derivative goes to 1.



---

*Source: 954015-2013-04-04.xml*