# An Efficient Series Solution for Nonlinear Multiterm Fractional Differential Equations

Authors: Moh’d Khier Al-Srihin; Mohammed Al-Refai
Journal: Discrete Dynamics in Nature and Society
 (2017)

---

## Abstract

In this paper, we introduce an efficient series solution for a class of nonlinear multiterm fractional differential equations of Caputo type. The approach is a generalization to our recent work for single fractional differential equations. We extend the idea of the Taylor series expansion method to multiterm fractional differential equations, where we overcome the difficulty of computing iterated fractional derivatives, which are difficult to be computed in general. The terms of the series are obtained sequentially using a closed formula, where only integer derivatives have to be computed. Several examples are presented to illustrate the efficiency of the new approach and comparison with the Adomian decomposition method is performed.


---

## Body

## 1. Introduction

Fractional differential equations (FDEs) are generalization to differential equations (DEs) for noninteger orders. In recent years, FDEs caught the attention of many researchers because of their appearance in modeling several phenomenon in the physical sciences. As many FDEs do not possess exact solutions on closed forms, analytical and numerical techniques have been implemented to study these equations. Iterative methods, such as the variational iteration method (VIM) in [1], the homotopy analysis method (HAM) in [2, 3], the Adomian decomposition method (ADM) in [4–9], and the fractional differential transform method in [10], have been implemented to solve various types of FDEs. These methods produce a solution in a series form whose terms are determined sequentially. We refer the reader to [11, 12] for a comprehensive study of series solutions of fractional differential equations. Recently, we have introduced a new series solution for single fractional differential equations [13]. The new approach is a modified form of the well-known Taylor series expansion, where we overcome the difficulty of computing iterative fractional derivatives. The efficiency of the new approach has been illustrated through several examples. In this paper we extend the idea to multiterm fractional differential equations. The presented work is a part of the Master thesis [14]. We consider the left Caputo fractional derivative D0+α, defined by [15, 16](1)D0+αft=1Γn-α∫0tt-τn-α-1fnτdτ,ifn-1<α<n,fnt,ifα=n∈Nprovided the integral exists, where Γ is the well-known Gamma function. The left Riemann-Liouville fractional integral, I0+α, of order α∈R+, is defined by(2)I0+αft=1Γα∫0tt-sα-1fsds,t>0,α>0.The left Caputo derivative is related to the left Riemann-Liouville fractional integral by(3)D0+αft=I0+n-αfnt,n-1<α<n.It is known that(4)I0+αD0+αft=ft-∑k=0m-1fk0k!tk,D0+αI0+αft=ft.Also, for n-1<α<n,n∈N,t>0, it holds that(5)D0+αtμ=Γμ+1Γμ-α+1tμ-α,ifμ>n-1,0,ifμ∈0,1,…,n-1.This paper is organized as follows. In Section 2, we present the series solution of nonlinear two-term fractional differential equations. We illustrate the efficiency of the presented technique through several examples. We also compare our results with the ones obtained by the Adomian decomposition method. In Section 3, we present and illustrate the efficiency of the new series solution for three-term fractional differential equations of several types. Finally, we conclude with some remarks in Section 4.


## 2. Two-Term Fractional Differential Equations

We start with the nonlinear two-term fractional initial value problems of the form(6)c1D0+α1ut+c2D0+α2ut=ft,ut,t>0with(7)u0=b,where 0<α2≤α1<1, and c1 and c2 are nonzero constants. We assume that f(t,u(t)) is continuous and smooth with respect to u(t). We also assume that α1 and α2 are rational numbers with α1=p1/q1 and α2=p2/q2, p1,p2,q1,q2∈N.


### 2.1. The Expansion Procedure

Letq=lcm(q1,q2); we have q=sq1=rq2 for some s, r∈N.In the following we expand the solution of problem (6)-(7) in an infinite series of the form (8)ut=∑n=0∞antn/q,where the coefficients an: n⩾0 have to be determined sequentially in the following manner: From the initial condition (7) we have u(0)=b=a0. Since D0+α(c)=0, for c being constant, we have(9)D0+α1ut=∑n=1∞ansntn/q-p1/q1=∑n=1∞ansntn-sp1/q,D0+α2ut=∑n=1∞anrntn/q-p2/q2=∑n=1∞anrntn-rp2/q,where(10)sn=Γn/q+1Γn/q-α1+1,rn=Γn/q+1Γn/q-α2+1.By substituting (9) in (6) we have(11)c1∑n=1∞ansntn-sp1/q+c2∑n=1∞anrntn-rp2/q=ft,∑n=0∞antn/q.Applying the well-known Taylor series method to compute the coefficients {an;n≥1} will lead to computing iterated fractional derivatives, which are not easily computed in general. To avoid this difficulty, let t=wq; we have(12)c1∑n=1∞ansnwn-sp1+c2∑n=1∞anrnwn-rp2=fwq,∑n=0∞anwn.Shifting the index to zero yields(13)c1∑n=0∞an+1sn+1wn-sp1+1+c2∑n=0∞an+1rn+1wn-rp2+1=fwq,∑n=0∞anwn.To avoid the singularity at w=0, we multiply (13) by wsp1-1; we have(14)c1∑n=0∞an+1sn+1wn+c2∑n=0∞an+1rn+1wn-rp2+sp1=wsp1-1fwq,∑n=0∞anwn.Now, since α1=p1/q1=sp1/q>rp2/q=p2/q2=α2, thus sp1-rp2>0, and (14) has no singularity at w=0.Letk=sp1-rp2-1≥0, and g(w)=f(wq,∑n=0∞‍anwn); then (14) can be written as(15)c1∑n=0kan+1sn+1wn+∑n=k+1∞c1an+1sn+1+c2an-krn-kwn=wk+rp2gw.We first determine the coefficients an for n≤k. By performing the nth derivative of (15) with respect to w and substituting w=0, we have(16)c1n!an+1sn+1=dndwnwk+rp2gww=0which yields(17)an+1=1c1n!sn+1dndwnwk+rp2fwq,∑n=0∞anwnw=0.Since k+rp2≥n+1, and f is smooth, then for n≤k, we have (18)dndwnwk+rp2fwq,∑n=0∞anwnw=0=0,and hence an+1=0, for n≤k.We now determinean, for n≥k+1. By performing the nth derivative of (15) with respect to w and substituting w=0, we have(19)n!c1an+1sn+1+c2an-krn-k=dndwnwk+rp2gww=0=dndwnwk+rp2fwq,∑n=0∞anwnw=0.Using the well-known Leibniz rule for differentiating the products, we have (20)dndwnwk+rp2gw=∑j=0nnjdjdwjwk+rp2dn-jdwn-jgw.Since(21)djdwjwk+rp2w=0=0ifj≠k+rp2,j!ifj=k+rp2,we have(22)c1an+1sn+1+c2an-krn-k=1n-j!dn-jdwn-jgww=0,wherej=k+rp2.From the last equation we determine an:n≥k+1 and thus the solution (23)ut=∑n=0∞antn/q=a0+∑n=k+1∞antn/qis obtained.Remark 1.
In (6), assuming α1=α2=1, then sn=rn=n,k=-1 and k+rp2=0. Thus (17) is reduced to (24)an+1=1c1n+1!dndwnfw,∑n=0∞anwnw=0,which coincides with the coefficients obtained by the Taylor series expansion method. Same comment applies for the coefficient an+1 in (22).Remark 2.
The algorithm can be generalized for the two-term fractional differential equation (6) with 0<α2≤α1<N, for arbitrary integer N. But we have to take care of the existence of the fractional derivative D0+α1(tn/q) (see Eq. (5)), by choosing the coefficient an=0, for n=1,…,qN-1. This case is discussed in Section 3.


### 2.2. Numerical Results

Example 3.
Consider the nonlinear two-term fractional initial value problem(25)2D0+1/2ut+2Γ1310D0+1/5ut=Γ12u2+t3/10-t+1,t>0,with(26)u0=0.The exact solution of the problem is u(t)=t1/2.
Applying the current algorithm, we haveα1=1/2=p1/q1, α2=1/5=p2/q2, q=l.c.m(q1,q2)=10, s=5, and r=2. We expand the solution in an infinite series of the form(27)ut=∑n=0∞antn/10.The initial condition in (26) yields a0=0. We have t=wq=w10 and(28)gw=fw10,uw=Γ12∑n=0∞anwn2+w3-w10+1.Since g(w) is continuous and smooth with respect to w, we have (29)an+1=0,forn≤k=sp1-rp2-1=2.Thus a1=a2=a3=0. The function g(w) satisfies the assumption of the proposed algorithm, and it holds that(30)dmdwmgww=0=Γ121m=03!m=36!a32m=62×7!a3a4m=78!a42+2×8!a3a5m=82!a4a5+a3a3m=9-10!+10!a52+2×10!a4a6+a3a7m=10⋮⋮⁡The computation above is made using the software Mathematica version 9. For n≥3; substituting (30) in (22) yields(31)2an+1sn+1+2Γ1310an-2rn-2=1n-4!dn-4dwn-4gww=0,where(32)sn+1=Γn+1/10+1Γn+1/10+1/2,rn-2=Γn-2/10+1Γn-2/10+4/5.Applying (31) together with a1=a2=a3=0, we have (33)a4=0,a5=1,an=0,n≥6.
Thus,(34)ut=a5t5/10=t1/2,and the exact solution of problem (25)-(26) is obtained.In the following we compare our results with the Adomian decomposition method (ADM). Assume that the nonlinear functionf(u)=u2 and the solution u(t) can be expressed in the following series form: (35)ut=∑n=0∞unt,fu=∑n=0∞An,where An,n=0,1,2,… are the well-known Adomian polynomials that can be derived from [17] (36)An=1n!dndβnf∑j=0∞βjujβ=0,j≥0.For recent advances in the Adomian decomposition method we refer the reader to [18, 19]. For f(u)=u2, we have(37)A0=u02,A1=2u0u1,A2=2u0u2+u12,A3=2u0u3+2u1u2,A4=2u0u4+2u1u3+u22,⋮Applying the Riemann-Liouville fractional integral operator I0+1/2 to (25) and substituting (38)I0+1/2D0+1/5ut=I0+3/10ut,we have(39)ut=Γ1/22I0+1/21+t3/10-t-Γ1310I0+3/10ut+Γ1/22I0+1/2u2t=t1/2+0.853958t4/5-0.666667t3/2-Γ1310I0+3/10∑n=0∞unt+Γ1/22I0+1/2∑n=0∞An=∑n=0∞unt.Set(40)u0=t1/2+0.853958t4/5-0.666667t3/2and balancing (39) yields(41)un+1=-Γ1310I0+3/10un+Γ1/22I0+1/2An,n=0,1,2,…Evaluating the first 5 term of the power series solution, we obtain the approximate solution u(4)(t)=∑n=04‍un(t), where(42)u4t=t0.5+0.257998t2-2.798359t2.7-1.454522t3+8.855605t3.4+9.799202t3.7+2.317936t4-11.599906t4.1-22.429449t4.4-11.273123t4.7+6.685081t4.8-1.388180t5+23.109057t5.1+20.278401t5.4-1.404129t5.5+5.305238t5.7-11.070905t5.8+0.278605t6-17.334122t6.1-7.918245t6.4+2.012699t6.5-0.885674t6.7+7.156236t6.8+5.827855t7.1+1.140180t7.4-1.153180t7.5-2.127835t7.8-0.742611t8.1+0.309488t8.5+0.244589t8.8-0.032578t9.5, which is not the exact solution. Solving the problem by the ADM with u0=t0.5 yields U(4)(t)=∑n=04‍Un(t), where(43)U4t=t0.5-0.853958t0.8+0.682106t1.1-0.515733t1.4+0.666667t1.5+0.0372208t1.7-1.527774t1.8+2.25291t2.1-2.68363t2.4+0.711111t2.5-2.576649t2.8+5.473545t3.1+0.8533333t3.5-4.19714t3.8+1.078789t4.5,which is again not the exact solution. Figure 1 depicts the exact solution u(t) obtained by the current algorithm and the approximate solutions u(4)(t) and U(4)(t) obtained by the Adomian decomposition method.Figure 1
Comparison between the exact and approximate solutions obtained by the ADM for problem ((25)-(26)).Example 4.
Consider the nonlinear two-term fractional initial value problem(44)D0+4/5ut+D0+1/2ut=utsin⁡tu2t+1,t>0,with(45)u0=1.This example has been discussed in [20], where the problem is transformed to a fractional integral equation, and then the Adams-Bashforth-Moulton method is used with step size h=1/50 to approximate the solution.
Applying the current algorithm, we havep1=4,q1=5,p2=1,q2=2,q=l.c.m(5,2)=10,s=2,r=5,k=2, and(46)an+1sn+1+an-2rn-2=1n-7!dn-7dwn-7fw10,∑n=0∞anwnw=0,where (47)fw10,∑n=0∞anwn=sin⁡w10∑n=0∞anwn∑n=0∞anwn2+1.The initial condition u(0)=1 yields a0=1. We apply (46) to compute the first 13 nonzero terms of an, and we have (48)ut=1+0.298242t1.8-0.227519t2.1+0.167717t2.4-0.119885t2.7+0.0833333t3-0.0564631t3.3+0.0373656t3.6-0.0280303t3.8-0.0241927t3.9+0.0179008t4.1+0.0153477t4.2-0.0112111t4.4.Since the exact solution of problem (44)-(45) is not available in a closed form, we define the error EN(t) by(49)ENti=∫0tiPuNt2dt,t∈0,ti,where (50)Put=D0+4/5ut+D0+1/2ut-utsin⁡tu2t+1,uNt=∑k=0Naktk/10.Tables 1 and 2 present the error EN for different values of N. One can see that the error decreases with N and more accuracy can be achieved by considering more terms. Also, the error increases with t, as any other series solutions.Table 1
The error of Example4 for N=35 and different values of t.
t
i
E
35
(
t
i
)
0.1
1.16844
×
1
0
-
10
0.2
8.36081
×
1
0
-
9
0.3
9.93967
×
1
0
-
8
0.4
5.28261
×
1
0
-
7
0.5
1.93074
×
1
0
-
6
0.6
5.47664
×
1
0
-
6
0.7
0.0000130833
0.8
0.0000276624
0.9
0.0000534896
1.0
0.0000967893Table 2
The error of Example4 at t=1 and for different values of N.
N
EN(1)
5
0.0681689
10
0.0681689
15
0.0681689
20
0.067155
25
0.0248829
30
0.00882027
35
0.0000967893


## 2.1. The Expansion Procedure

Letq=lcm(q1,q2); we have q=sq1=rq2 for some s, r∈N.In the following we expand the solution of problem (6)-(7) in an infinite series of the form (8)ut=∑n=0∞antn/q,where the coefficients an: n⩾0 have to be determined sequentially in the following manner: From the initial condition (7) we have u(0)=b=a0. Since D0+α(c)=0, for c being constant, we have(9)D0+α1ut=∑n=1∞ansntn/q-p1/q1=∑n=1∞ansntn-sp1/q,D0+α2ut=∑n=1∞anrntn/q-p2/q2=∑n=1∞anrntn-rp2/q,where(10)sn=Γn/q+1Γn/q-α1+1,rn=Γn/q+1Γn/q-α2+1.By substituting (9) in (6) we have(11)c1∑n=1∞ansntn-sp1/q+c2∑n=1∞anrntn-rp2/q=ft,∑n=0∞antn/q.Applying the well-known Taylor series method to compute the coefficients {an;n≥1} will lead to computing iterated fractional derivatives, which are not easily computed in general. To avoid this difficulty, let t=wq; we have(12)c1∑n=1∞ansnwn-sp1+c2∑n=1∞anrnwn-rp2=fwq,∑n=0∞anwn.Shifting the index to zero yields(13)c1∑n=0∞an+1sn+1wn-sp1+1+c2∑n=0∞an+1rn+1wn-rp2+1=fwq,∑n=0∞anwn.To avoid the singularity at w=0, we multiply (13) by wsp1-1; we have(14)c1∑n=0∞an+1sn+1wn+c2∑n=0∞an+1rn+1wn-rp2+sp1=wsp1-1fwq,∑n=0∞anwn.Now, since α1=p1/q1=sp1/q>rp2/q=p2/q2=α2, thus sp1-rp2>0, and (14) has no singularity at w=0.Letk=sp1-rp2-1≥0, and g(w)=f(wq,∑n=0∞‍anwn); then (14) can be written as(15)c1∑n=0kan+1sn+1wn+∑n=k+1∞c1an+1sn+1+c2an-krn-kwn=wk+rp2gw.We first determine the coefficients an for n≤k. By performing the nth derivative of (15) with respect to w and substituting w=0, we have(16)c1n!an+1sn+1=dndwnwk+rp2gww=0which yields(17)an+1=1c1n!sn+1dndwnwk+rp2fwq,∑n=0∞anwnw=0.Since k+rp2≥n+1, and f is smooth, then for n≤k, we have (18)dndwnwk+rp2fwq,∑n=0∞anwnw=0=0,and hence an+1=0, for n≤k.We now determinean, for n≥k+1. By performing the nth derivative of (15) with respect to w and substituting w=0, we have(19)n!c1an+1sn+1+c2an-krn-k=dndwnwk+rp2gww=0=dndwnwk+rp2fwq,∑n=0∞anwnw=0.Using the well-known Leibniz rule for differentiating the products, we have (20)dndwnwk+rp2gw=∑j=0nnjdjdwjwk+rp2dn-jdwn-jgw.Since(21)djdwjwk+rp2w=0=0ifj≠k+rp2,j!ifj=k+rp2,we have(22)c1an+1sn+1+c2an-krn-k=1n-j!dn-jdwn-jgww=0,wherej=k+rp2.From the last equation we determine an:n≥k+1 and thus the solution (23)ut=∑n=0∞antn/q=a0+∑n=k+1∞antn/qis obtained.Remark 1.
In (6), assuming α1=α2=1, then sn=rn=n,k=-1 and k+rp2=0. Thus (17) is reduced to (24)an+1=1c1n+1!dndwnfw,∑n=0∞anwnw=0,which coincides with the coefficients obtained by the Taylor series expansion method. Same comment applies for the coefficient an+1 in (22).Remark 2.
The algorithm can be generalized for the two-term fractional differential equation (6) with 0<α2≤α1<N, for arbitrary integer N. But we have to take care of the existence of the fractional derivative D0+α1(tn/q) (see Eq. (5)), by choosing the coefficient an=0, for n=1,…,qN-1. This case is discussed in Section 3.


## 2.2. Numerical Results

Example 3.
Consider the nonlinear two-term fractional initial value problem(25)2D0+1/2ut+2Γ1310D0+1/5ut=Γ12u2+t3/10-t+1,t>0,with(26)u0=0.The exact solution of the problem is u(t)=t1/2.
Applying the current algorithm, we haveα1=1/2=p1/q1, α2=1/5=p2/q2, q=l.c.m(q1,q2)=10, s=5, and r=2. We expand the solution in an infinite series of the form(27)ut=∑n=0∞antn/10.The initial condition in (26) yields a0=0. We have t=wq=w10 and(28)gw=fw10,uw=Γ12∑n=0∞anwn2+w3-w10+1.Since g(w) is continuous and smooth with respect to w, we have (29)an+1=0,forn≤k=sp1-rp2-1=2.Thus a1=a2=a3=0. The function g(w) satisfies the assumption of the proposed algorithm, and it holds that(30)dmdwmgww=0=Γ121m=03!m=36!a32m=62×7!a3a4m=78!a42+2×8!a3a5m=82!a4a5+a3a3m=9-10!+10!a52+2×10!a4a6+a3a7m=10⋮⋮⁡The computation above is made using the software Mathematica version 9. For n≥3; substituting (30) in (22) yields(31)2an+1sn+1+2Γ1310an-2rn-2=1n-4!dn-4dwn-4gww=0,where(32)sn+1=Γn+1/10+1Γn+1/10+1/2,rn-2=Γn-2/10+1Γn-2/10+4/5.Applying (31) together with a1=a2=a3=0, we have (33)a4=0,a5=1,an=0,n≥6.
Thus,(34)ut=a5t5/10=t1/2,and the exact solution of problem (25)-(26) is obtained.In the following we compare our results with the Adomian decomposition method (ADM). Assume that the nonlinear functionf(u)=u2 and the solution u(t) can be expressed in the following series form: (35)ut=∑n=0∞unt,fu=∑n=0∞An,where An,n=0,1,2,… are the well-known Adomian polynomials that can be derived from [17] (36)An=1n!dndβnf∑j=0∞βjujβ=0,j≥0.For recent advances in the Adomian decomposition method we refer the reader to [18, 19]. For f(u)=u2, we have(37)A0=u02,A1=2u0u1,A2=2u0u2+u12,A3=2u0u3+2u1u2,A4=2u0u4+2u1u3+u22,⋮Applying the Riemann-Liouville fractional integral operator I0+1/2 to (25) and substituting (38)I0+1/2D0+1/5ut=I0+3/10ut,we have(39)ut=Γ1/22I0+1/21+t3/10-t-Γ1310I0+3/10ut+Γ1/22I0+1/2u2t=t1/2+0.853958t4/5-0.666667t3/2-Γ1310I0+3/10∑n=0∞unt+Γ1/22I0+1/2∑n=0∞An=∑n=0∞unt.Set(40)u0=t1/2+0.853958t4/5-0.666667t3/2and balancing (39) yields(41)un+1=-Γ1310I0+3/10un+Γ1/22I0+1/2An,n=0,1,2,…Evaluating the first 5 term of the power series solution, we obtain the approximate solution u(4)(t)=∑n=04‍un(t), where(42)u4t=t0.5+0.257998t2-2.798359t2.7-1.454522t3+8.855605t3.4+9.799202t3.7+2.317936t4-11.599906t4.1-22.429449t4.4-11.273123t4.7+6.685081t4.8-1.388180t5+23.109057t5.1+20.278401t5.4-1.404129t5.5+5.305238t5.7-11.070905t5.8+0.278605t6-17.334122t6.1-7.918245t6.4+2.012699t6.5-0.885674t6.7+7.156236t6.8+5.827855t7.1+1.140180t7.4-1.153180t7.5-2.127835t7.8-0.742611t8.1+0.309488t8.5+0.244589t8.8-0.032578t9.5, which is not the exact solution. Solving the problem by the ADM with u0=t0.5 yields U(4)(t)=∑n=04‍Un(t), where(43)U4t=t0.5-0.853958t0.8+0.682106t1.1-0.515733t1.4+0.666667t1.5+0.0372208t1.7-1.527774t1.8+2.25291t2.1-2.68363t2.4+0.711111t2.5-2.576649t2.8+5.473545t3.1+0.8533333t3.5-4.19714t3.8+1.078789t4.5,which is again not the exact solution. Figure 1 depicts the exact solution u(t) obtained by the current algorithm and the approximate solutions u(4)(t) and U(4)(t) obtained by the Adomian decomposition method.Figure 1
Comparison between the exact and approximate solutions obtained by the ADM for problem ((25)-(26)).Example 4.
Consider the nonlinear two-term fractional initial value problem(44)D0+4/5ut+D0+1/2ut=utsin⁡tu2t+1,t>0,with(45)u0=1.This example has been discussed in [20], where the problem is transformed to a fractional integral equation, and then the Adams-Bashforth-Moulton method is used with step size h=1/50 to approximate the solution.
Applying the current algorithm, we havep1=4,q1=5,p2=1,q2=2,q=l.c.m(5,2)=10,s=2,r=5,k=2, and(46)an+1sn+1+an-2rn-2=1n-7!dn-7dwn-7fw10,∑n=0∞anwnw=0,where (47)fw10,∑n=0∞anwn=sin⁡w10∑n=0∞anwn∑n=0∞anwn2+1.The initial condition u(0)=1 yields a0=1. We apply (46) to compute the first 13 nonzero terms of an, and we have (48)ut=1+0.298242t1.8-0.227519t2.1+0.167717t2.4-0.119885t2.7+0.0833333t3-0.0564631t3.3+0.0373656t3.6-0.0280303t3.8-0.0241927t3.9+0.0179008t4.1+0.0153477t4.2-0.0112111t4.4.Since the exact solution of problem (44)-(45) is not available in a closed form, we define the error EN(t) by(49)ENti=∫0tiPuNt2dt,t∈0,ti,where (50)Put=D0+4/5ut+D0+1/2ut-utsin⁡tu2t+1,uNt=∑k=0Naktk/10.Tables 1 and 2 present the error EN for different values of N. One can see that the error decreases with N and more accuracy can be achieved by considering more terms. Also, the error increases with t, as any other series solutions.Table 1
The error of Example4 for N=35 and different values of t.
t
i
E
35
(
t
i
)
0.1
1.16844
×
1
0
-
10
0.2
8.36081
×
1
0
-
9
0.3
9.93967
×
1
0
-
8
0.4
5.28261
×
1
0
-
7
0.5
1.93074
×
1
0
-
6
0.6
5.47664
×
1
0
-
6
0.7
0.0000130833
0.8
0.0000276624
0.9
0.0000534896
1.0
0.0000967893Table 2
The error of Example4 at t=1 and for different values of N.
N
EN(1)
5
0.0681689
10
0.0681689
15
0.0681689
20
0.067155
25
0.0248829
30
0.00882027
35
0.0000967893


## 3. Three-Term Fractional Differential Equations

We consider the nonlinear three-term fractional initial value problem of the form(51)c1D0+α1ut+c2D0+α2ut+c3D0+α3ut=ft,yt,t>0,with(52)ui0=bi,i=0,1,…,N1-1,where(53)0<α3≤α2≤α1<N1,N1-1<α1<N1,N2-1<α2≤N2≤N1,0<N3-1<α3≤N3≤N2,and c1, c2, and c3 are nonzero constants. We assume that f(t,y(t)) is continuous and smooth with respect to u(t). We also assume that α1,α2, and α3 are rational numbers with α1=p1/q1,α2=p2/q2, and α3=p3/q3, where p1,p2,p3,q1,q2,q3∈N. Let q=lcm(q1,q2,q3); we have q=sq1=rq2=vq3 for s, r and v∈N.Applying analogous steps for the case of the two-term fractional equations, we seek a solution of the form(54)ut=∑n=0∞antn/q=∑n=0N1-1qantn/q+∑n=N1-1q+1∞antn/q=∑n=0N1-1anqtn+∑n=N1-1q+1∞antn/q. From the initial condition (52) we have u(i)(0)=bi=i!aiq, and thus, aiq=bi/i!,i=0,1,…,N1-1.Using the fact thatD0+α(tn)=0, for α>n∈N, we have(55)D0+α1ut=∑n=N1-1q+1∞ansntn-sp1/q,D0+α2ut=∑n=N2N1-1anqρnqtn-α2+∑n=N1-1q+1∞anrntn-rp2/q,D0+α3ut=∑n=N3N1-1anqνnqtn-α3+∑n=N1-1q+1∞anvntn-vp3/q,where sn=Γ(n/q+1)/Γ(n/q-α1+1),rn=Γ(n/q+1)/Γ(n/q-α2+1),ρn=Γ(n+1)/Γ(n+1-α2), vn=Γ(n/q+1)/Γ(n/q-α3+1), and νn=Γ(n+1)/Γ(n+1-α3).Substituting (55) in (51), we have(56)c1∑n=N1-1q+1∞ansntn-sp1/q+c2∑n=N2N1-1anqρnqtn-α2+c2∑n=N1-1q+1∞anrntn-rp2/q+c3∑n=N3N1-1anqνnqtn-α3+c3∑n=N1-1q+1∞anvntn-vp3/q=ft,∑n=0∞antn/q.Since the coefficients anq,n=0,…,N1-1, are known, we write the last equation in the form(57)c1∑n=N1-1q+1∞ansntn-sp1/q+c2∑n=N1-1q+1∞anrntn-rp2/q+c3∑n=N1-1q+1∞anvntn-vp3/q=ft,∑n=0∞antn/q-c2∑n=N2N1-1anqρnqtn-α2-c3∑n=N3N1-1anqνnqtn-α3.Let w=t1/q; then (57) reduces to(58)∑n=N1-1q+1∞c1ansnwn-sp1+∑n=N1-1q+1∞c2anrnwn-rp2+∑n=N1-1q+1∞c3anvnwn-vp3=fwq,∑n=0∞anwn-∑n=N2N1-1c2anqρnqwqn-α2-∑n=N3N1-1c3anqνnqwqn-α3.Substituting z=(N1-1)q+1 and(59)gw=fwq,∑n=0∞anwn-∑n=N2N1-1c2anqρnqwqn-α2-∑n=N3N1-1c3anqνnqwqn-α3and shifting the indices to zero, we have(60)∑n=0∞c1an+zsn+zwn-sp1+z+∑n=0∞c2an+zrn+zwn-rp2+z+∑n=0∞c3an+zvn+zwn-vp3+z=gw.To avoid singularity we multiply both sides in (60) by wsp1-z; we have(61)∑n=0∞c1an+zsn+zwn+∑n=0∞c2an+zrn+zwn-rp2+sp1+∑n=0∞c3an+zvn+zwn-vp3+sp1=wsp1-zgw.Equation (61) can be written as(62)∑n=0k1-1c1an+zsn+zwn+∑n=k1k2-1c1an+zsn+z+c2an+z-k1rn+z-k1wn+∑n=k2∞c1an+zsn+z+c2an+z-k1rn+z-k1+c3an+z-k2vn+z-k2wn=wsp1-zgw,where k1=sp1-rp2 and k2=sp1-vp3.Applying analogs steps for deriving (22), we have(63)c1an+zsn+z+c2an+z-k1rn+z-k1+c3an+z-k2vn+z-k2=1n-j!dn-jdwn-jgww=0forn≥k2,where j=sp1-z. From the last recursion we determine an+z,n≥k2.Remark 5.
Our derivation is based on the facts that the nonlinear function is smooth and the fractional differential equations is of constant coefficients. In case if one of these conditions does not hold, a modified treatment will be considered as we will see in Example7.


### 3.1. Numerical Results

Example 6.
Consider the Bagely-Torvik initial value problem(64)D0+5/2ut+D0+2ut-2πD0+1/2ut+4ut=gt,t∈0,1,(65)u0=u′0=u′′0=0,where(66)gt=4t9-13107212155t17/2+72t7+49152143πt13/2.
This example has been discussed in [21] using a Chebyshev spectral method, where the solution has been approximated by the shifted Chebyshev polynomials with different degrees. Then the exact solution u(t)=t9 was obtained by considering the shifted Chebyshev polynomial of degree 9. We mention here that there are simple typos in presenting the example in (64) and in (66) and we correct them here. Applying the current algorithm we have (67)p1=5,p2=4,p3=1,q1=q2=q3=2,q=lcmq1,q2,q3=2,s=1,r=1,v=1,k1=1,k2=4.We expand the solution in infinite series of the form u(t)=∑n=0∞‍antn/2. The initial condition in (65) yields a0=a2=a4=0. Let t=w2; applying (59) we have (68)gw=4w18-13107212155w17+72w14+49152143πw13-4∑n=0∞anwn, and it holds that(69)dmdwmgww=0=-4a1ifm=1-4×2!a2ifm=2⋮⋮13!49152143π-4a13ifm=1314!72-4a14ifm=14-4×15!a15ifm=15-4×16!a16ifm=1617!-13107212155-4a17ifm=1718!4-4a18ifm=18-4!a19ifm=19⋮⋮Using (63) and (69), we have(70)an+5sn+5+an+4rn+4-2πan+1vn+1=1n!dndwngww=0,n≥4,where (71)sn+5=Γn+5/2+1Γn+5/2-3/2,rn+4=n+4n+24,vn+1=Γn+1/2+1Γn+1/2+1/2.Since g(w) is smooth, then an+5=0, for n<4. We now apply the last recursion together with a0=⋯=a8=0, to compute an+5, for n≥4. For n=4, we have (72)a9s9+a8r8-2πa5v5=g40=0,and thus a9=0.
Applying analogous arguments yieldsa10=⋯=a17=0. For n=13, we have (73)a18s18+a17r17-2πa14v14=49152143π-4a13,which, together with a17=a14=a13=0, implies that a18=1. For n=14, we have (74)a19s19+a18r18-2πa15v15=72-4a14,but a18=1 and a15=a14=0; thus a19=0. Applying the same steps yields a20=a21=0. For n=17, we have (75)a22s22+a21r21-2πa18v18=-13107212155-4a17,and thus a22=0.
Proceeding in the same manner, we havean+1=0 for n≥18. Thus, (76)ut=a18t18/2=t9,and the exact solution of problem (64)-(65) is obtained.In the following example we show that the current algorithm can be applied to more general multiterm fractional differential equations which are not necessary of the form in (51).Example 7.
Consider the nonlinear three-term fractional initial value problem(77)tD0+4/3ut+D0+1/3ut+t19/6D0+1/6ut=334Γ8/3t8/3+1296935Γ5/6u2t,with(78)u0=u′0=0.The exact solution for this problem is u(t)=t3.Applying the current algorithm we haveq=6. We then expand the solution in an infinite series of the form u(t)=∑n=0∞‍antn/6. From the initial condition (78), we have a0=a6=0. To guarantee the existence of the fractional derivative D0+4/3u(t) we choose a1=⋯=a5=0. Let w=t1/6; then (77) reduces to(79)a1r1+∑n=018an+2sn+2+an+2rn+2wn+∑n=19∞an+2sn+2+an+2rn+2+an-18vn-18wn=33w164Γ8/3+1296935Γ5/6∑n=0∞anwn2,where sn=Γn/6+1/Γn/6-1/3, rn=Γ(n/6+1)/Γ(n/6+2/3), and vn=Γ(n/6+1)/Γ(n/6+5/6).Performing thenth derivatives one can see that an=0 for n<16. For n=16, we have (80)a18s18+a18r18=334Γ8/3,which yields a18=1.Also,a19=⋯=a35=0.Forn=36, we have (81)a38s38+a38r38+a18v18=1296935Γ5/6,which yields a38=0.Following analogous steps, we havean=0 for n>36. Thus (82)ut=a18t18/6=t3,which is the exact solution.


## 3.1. Numerical Results

Example 6.
Consider the Bagely-Torvik initial value problem(64)D0+5/2ut+D0+2ut-2πD0+1/2ut+4ut=gt,t∈0,1,(65)u0=u′0=u′′0=0,where(66)gt=4t9-13107212155t17/2+72t7+49152143πt13/2.
This example has been discussed in [21] using a Chebyshev spectral method, where the solution has been approximated by the shifted Chebyshev polynomials with different degrees. Then the exact solution u(t)=t9 was obtained by considering the shifted Chebyshev polynomial of degree 9. We mention here that there are simple typos in presenting the example in (64) and in (66) and we correct them here. Applying the current algorithm we have (67)p1=5,p2=4,p3=1,q1=q2=q3=2,q=lcmq1,q2,q3=2,s=1,r=1,v=1,k1=1,k2=4.We expand the solution in infinite series of the form u(t)=∑n=0∞‍antn/2. The initial condition in (65) yields a0=a2=a4=0. Let t=w2; applying (59) we have (68)gw=4w18-13107212155w17+72w14+49152143πw13-4∑n=0∞anwn, and it holds that(69)dmdwmgww=0=-4a1ifm=1-4×2!a2ifm=2⋮⋮13!49152143π-4a13ifm=1314!72-4a14ifm=14-4×15!a15ifm=15-4×16!a16ifm=1617!-13107212155-4a17ifm=1718!4-4a18ifm=18-4!a19ifm=19⋮⋮Using (63) and (69), we have(70)an+5sn+5+an+4rn+4-2πan+1vn+1=1n!dndwngww=0,n≥4,where (71)sn+5=Γn+5/2+1Γn+5/2-3/2,rn+4=n+4n+24,vn+1=Γn+1/2+1Γn+1/2+1/2.Since g(w) is smooth, then an+5=0, for n<4. We now apply the last recursion together with a0=⋯=a8=0, to compute an+5, for n≥4. For n=4, we have (72)a9s9+a8r8-2πa5v5=g40=0,and thus a9=0.
Applying analogous arguments yieldsa10=⋯=a17=0. For n=13, we have (73)a18s18+a17r17-2πa14v14=49152143π-4a13,which, together with a17=a14=a13=0, implies that a18=1. For n=14, we have (74)a19s19+a18r18-2πa15v15=72-4a14,but a18=1 and a15=a14=0; thus a19=0. Applying the same steps yields a20=a21=0. For n=17, we have (75)a22s22+a21r21-2πa18v18=-13107212155-4a17,and thus a22=0.
Proceeding in the same manner, we havean+1=0 for n≥18. Thus, (76)ut=a18t18/2=t9,and the exact solution of problem (64)-(65) is obtained.In the following example we show that the current algorithm can be applied to more general multiterm fractional differential equations which are not necessary of the form in (51).Example 7.
Consider the nonlinear three-term fractional initial value problem(77)tD0+4/3ut+D0+1/3ut+t19/6D0+1/6ut=334Γ8/3t8/3+1296935Γ5/6u2t,with(78)u0=u′0=0.The exact solution for this problem is u(t)=t3.Applying the current algorithm we haveq=6. We then expand the solution in an infinite series of the form u(t)=∑n=0∞‍antn/6. From the initial condition (78), we have a0=a6=0. To guarantee the existence of the fractional derivative D0+4/3u(t) we choose a1=⋯=a5=0. Let w=t1/6; then (77) reduces to(79)a1r1+∑n=018an+2sn+2+an+2rn+2wn+∑n=19∞an+2sn+2+an+2rn+2+an-18vn-18wn=33w164Γ8/3+1296935Γ5/6∑n=0∞anwn2,where sn=Γn/6+1/Γn/6-1/3, rn=Γ(n/6+1)/Γ(n/6+2/3), and vn=Γ(n/6+1)/Γ(n/6+5/6).Performing thenth derivatives one can see that an=0 for n<16. For n=16, we have (80)a18s18+a18r18=334Γ8/3,which yields a18=1.Also,a19=⋯=a35=0.Forn=36, we have (81)a38s38+a38r38+a18v18=1296935Γ5/6,which yields a38=0.Following analogous steps, we havean=0 for n>36. Thus (82)ut=a18t18/6=t3,which is the exact solution.


## 4. Conclusion

For fractional differential equations of ordern-1<α<n, it is common to obtain a series solution in the form ∑n=0∞‍antα(n). The question is how to obtain the coefficients an,n=0,1,…. Naturally, if the problem is of fractional order, the differentiation is also of fractional order. In this paper, we presented a new algorithm for obtaining a series solution for nonlinear multiterm fractional differential equations of Caputo type, where we overcome the use of fractional differentiation. We employed a transformation that allows us to use ordinary differentiation rather than fractional differentiation to recursively compute the coefficient of the series expansion. Then the terms of the series, an, are computed sequentially using a closed form formula. We applied the new algorithm to several types of multiterm fractional differential equations, where accurate solutions as well as exact solutions in closed forms have been obtained. For one example it is noted that the current algorithm is more efficient than the ADM as it is more easier to apply and it produces the exact solution while the ADM does not. We have developed the new algorithm for two- and three-term fractional equations, while the idea can be extended to multiterm fractional equations of arbitrary order; obtaining a general formula in this case is not an easy task. The current algorithm can be modified to deal with the fractional multiterm time-diffusion equations.



---

*Source: 5234151-2017-03-08.xml*