Advances in Mathematical
Sciences and Applications
Vol. 27, No. 1 (2018), pp. 169–180

GAKKOTOSHOTOKYO JAPANNUMERICAL EXPERIMENTS ON ANALYTICITY OF
SOLUTIONS TO FRACTIONAL DIFFERENTIAL
EQUATIONS
Hiroshi Fujiwara
Graduate School of Informatics, Kyoto University,
Yoshida-honmachi, Sakyo-ku, Kyoto, 606-8501, Japan(fujiwara@acs.i.kyoto-u.ac.jp)
Nobuyuki Higashimori
Center for the Promotion of Interdisciplinary Education and Research,
Kyoto University,
Yoshida-honmachi, Sakyo-ku, Kyoto, 606-8501, Japan(nobuyuki@acs.i.kyoto-u.ac.jp)
andHitoshi Imai
Faculty of Science and Engineering, Doshisha University,
1-3 Tatara Miyakodani, Kyotanabe, Kyoto, 610-0394, Japan(himai@mail.doshisha.ac.jp)

Abstract. In the paper, analyticity of solutions of fractional diﬀerential equations isinvestigated numerically. In numerical computation the Chebyshev spectral collocation
method is used for discretization to realize arbitrary order approximation. Multipleprecision arithmetic is also used for arbitrary reduction of rounding error. Some onedimensional fractional diﬀerential equations are solved numerically. Numerical results
clarify analyticity of solutions.

————————————————————
Communicated by Editors; Received April 17, 2018
This work was supported by JSPS KAKENHI Grant Numbers 16K13774 and 16H02155.
AMS Subject Classiﬁcation: 65-05, 65L05.
Keywords: fractional diﬀerential equation, analyticity, spectral method, multiple-precision arithmetic.

169

170

1

Introduction

Recently, there has been a tremendous increase in the use of fractional diﬀerential models
to simulate dynamics in many ﬁelds, e.g. physics, chemistry, biology, engineering and soon [1],[4],[6],[7],[17]. In the models fractional order derivatives are used. They have a long
history and they have been studied in the ﬁeld of fractional calculus. A list of prominent mathematicians contributing to fractional calculus is found in [20]. Besides them,
Caputo gave a convenient deﬁnition for initial value problems of diﬀerential equations asfollows [5]: For 0 < α non-integer,αc D0 u(t) =1
Γ(n − α)∫ t(t − s)n−α−10d n u(s)
ds,
dsn

(1)

where Γ is the gamma function, n = ⌈α⌉ is the least integer more than or equal to α, andu is assumed to be suﬃciently smooth. This deﬁnition is used in the paper because it isconvenient for our numerical computation.
It is hard to obtain exact solutions of fractional diﬀerential equations, or even if anexact solution is available it may be too complicated to be used in practice. Thus the roleof numerical computation becomes important. Various numerical methods to fractional
diﬀerential equations such as the Adams-type predictor-corrector method, linear multistep methods, are presented in [1] and [7]. Spectral methods have also been applied to
fractional diﬀerential equations [7],[9],[21]. As for spectral methods it is well known thatthey are superior in accuracy [3]. If a function u is m-times continuously diﬀerentiable,
then ∥u − uN ∥ = O(N −m ) where uN is the N -th order approximate function by the spectral methods, and ∥·∥ is a suitable norm. This property is called spectral accuracy or
inﬁnite-order accuracy. Moreover, if a function u is analytic in a strip of the complex
plane, then ∥u − uN ∥ = O(e−γN ) with some positive number γ. In [21] numerical results
for some fractional diﬀerential equations show spectral accuracy.
Spectral methods are roughly divided into the Galerkin type and the collocation type[2],[3]. The latter type is called the spectral collocation method, and it is more convenient
for nonlinear problems or higher dimensional problems. In the paper the Chebyshev spectral collocation method is adopted. Although in [7] and [21] spectral collocation methods
are already used, by using Chebyshev polynomials we can give a concrete expression ofcoeﬃcients in the linear combination of basis functions. In the paper multiple-precision
arithmetic is also used in numerical computation. To investigate mathematical properties
of solutions numerically it is necessary to eliminate the inﬂuence of rounding error. The
combination of spectral methods and multiple-precision arithmetic was presented in [16]to realize the arbitrary reduction of numerical error. Such reduction is indispensable fordirect numerical simulations of inverse problems [10],[13],[15]. By using this combination we also investigated numerically mathematical properties of functions or solutions,for instance, regularity or existence of solutions [8],[14],[19]. In the paper, by using the
combination we investigate analyticity of solutions of fractional diﬀerential equations.

171

2

The Chebyshev Spectral Collocation Method

We derive a convenient form of fractional derivatives in terms of the Chebyshev spectral
method, which plays an essential role in our numerical studies.
First we brieﬂy recall the Chebyshev spectral collocation method [2]. Let ŨN (x) be apolynomial function of degree N deﬁned on the interval [−1, 1]. The Chebyshev expansionof ŨN (x) isŨN (x) =N∑āk Tk (x),

(2)

k=0

where Tk (x) = cos (k arccos x) is the Chebyshev polynomial of order k. We also introduce
Chebyshev-Gauss-Lobatto collocation points as follows:x̄j = cosjπN(j = 0, 1, 2, · · · , N ).

(3)

Letting Ũj = ŨN (x̄j ), from (2) we have
Ũj =

N
∑

āk Tk (x̄j ).

k=0

The coeﬃcients ak are given by the following inverse formula
2 ∑ 1
āk =
Tk (x̄j )Ũj ,N ck j=0 cj
N(k = 0, 1, 2, . . . , N )

(4)

where
{
2, j = 0, N,
cj =1, otherwise.
Substituting (4) to āk in (2) we obtain
)
(
N
N
∑
1 ∑ 2
ŨN (x) =Tk (x̄j )Tk (x) Ũj .
cj k=0 N ck
j=0This representation motivates an N -th order approximation UN (x) to a function U (x)on [−1, 1] as
)
(
N
N
∑
1 ∑ 2
UN (x) =Tk (x̄j )Tk (x) Uj ,
(5)
c
N
c
j
k
j=0
k=0
where Uj = U (x̄j ). We should remark that if U (x) is a polynomial of degree m(≤ N ) thenU (x) = UN (x). The truncation number N in (5) denotes the approximation order. On theother hand, N in (3) denotes the number of collocation points; to be exact, the number
is (N + 1). This means that in spectral collocation methods the approximation ordercan be easily controlled by the number of collocation points. As stated in Introduction,

172
spectral collocation methods are applicable to nonlinear problems or higher dimensional
problems. These features are quite convenient to investigate mathematical properties ofsolutions to various problems.
In spectral collocation methods, diﬀerential equations are satisﬁed at the collocation
points. Thus they require the derivative of the approximate function at each collocation
point. For instance, the ﬁrst order derivative of UN (x) at the collocation point xi is given
by diﬀerentiating (5):
(
)
N
N
N
∑
∑
d
d
1 ∑ 2
Tk (x̄j ) Tk (x̄i ) Uj ,UN (x̄i ) =(Dx )i,j Uj =
(6)
dx
c
N
c
dx
j
k
j=0
j=0
k=0
where Dx is called the ﬁrst order derivative matrix.
T
Let u(t) be a function deﬁned on [0, T ]. By substituting t = (x + 1), u
2
is a function on [−1, 1] and is approximated by the Chebyshev expansion
 a()  .0 
UN (x) = T0 (x) · · · TN (x)  ..  .
aN( )By using a square matrix T̃ = T̃kj of order N + 1 with entries T̃kj =(T(x + 1)
2)

2
Tk (xj ), theN ck cj

inverse formula (4) in the matrix form is
 
 
a0
U0
 .. 
 .. 
 .  = T̃  .  ,aN(where Uj = uUN)T(xj + 1) . Let C be a square matrix of order N + 1 which satisﬁes
2
 


1
T0 (x)
 T (x + 1)   T1 (x) 
 
 2

C
(7)
..
 =  ..  .
  . .
( T )NTN (x)(x + 1)N
2

We remark that T̃ and C are constant matrices. They lead to the expression
 
U0
)
(
(
)


N
T
C T̃  ...  ,UN (x) = 1 T2 (x + 1) · · · T2 (x + 1)
UNT
where C T is the transpose of C. Let tj = (xj + 1), uj = u(tj ), and uN (t) = UN (x) with2T
t = (x + 1). Then uj = Uj and the function uN is an approximation to u and is written
2

173
as




u0
(
)
 
uN (t) = 1 t · · · tN C T T̃  ...  .
uNThis immediately leads to the fractional derivative of uN at ti asu0) T  . 
( α
α
α N
α
c D0 uN (ti ) = c D0 [1](ti ) c D0 [t](ti ) · · · c D0 [t ](ti ) C T̃  ..  ,uNwhich is a linear combination of uj , j = 0, . . . , N , and is similar to (6). It should beemphasized that the Caputo derivatives of the unknown function at each collocation
points c D0α [tk ](ti ) are easily calculated since the Caputo derivatives of monomials tk for anonnegative integer k are given by0,
0 ≤ k < α,α kc D0 [t ] = Γ(k+1) tk−α , α ≤ k.
Γ(k+1−α)It is noted in [21] that the entries of C are calculated by the use of inverse of the
Vandermonde matrix which is obtained by substituting N + 1 collocation points xj , j =0, 1, . . . , N , into (7). However, the entries of C are (calculated
explicitly)( ) as follows. For
the purpose, we introduce two square matrices F = fkj and L = ℓkj of order N + 1,
which satisfy
 
 




1
1
1
T0 (x)xx
 T (x + 1)  .. 
 
 

 2
..
.
 .  = F  ..  and  ..  = L 
 . 
 . 


.
( T )NTN (x)
NNNxx(x + 1)
2
Firstly, the entries fkj are calculated from the recurrence relation T0 (x) = 1, T1 (x) = x
and Tk+1 (x) = 2xTk (x) − Tk−1 (x). Secondly, it follows from the binomial expansion thatx = (x + 1 − 1) =
k

k

k
∑
j=0

(−1)

k−j

( )
k
(x + 1)j .
jSince (x + 1)j , j = 0, . . . , N , are linearly independent, ℓkj is given by( )j ( )k(−1)k−j 2, j ≤ k,
T
ℓkj =
j


0,
k < j.
Finally the product of two lower triangular matrices gives C = F L.

174

3

Numerical Evaluation of Fractional Derivatives

This section is devoted to investigate accuracy of numerical integration appearing in (1)and computational environments.
The integrand in (1) has a singularity at boundary s = t, and thus special attentions
should be paid for the sake of reliable and accurate computation. The use of double
exponential rule [18] is known as a novel method to overcome the diﬃculty. It is also
convenient due to returning numerical results in given tolerance automatically.
We demonstrate an example with 0 <(α < 1 andu(t) = t2 . Using the standard double)πexponential transformation ϕ(y) = tanhsinh y , the evaluation of integral at t = 1 is2
∫ ∞1
ϕ(y) + 1 ϕ′ (y)
α
(
)α
D
u(1)
=
dy
c 0
1−ϕ(y)
Γ(1 − α) −∞2
2

by introducing the discretization parameter h > 0 and yk = kh, k ∈ Z, “the composite
trapezoidal rule” gives an approximate series as≈∑1Γ(1 − α) k∈Zϕ(y ) + 1
( k )α1−ϕ(yk )2

ϕ′ (yk )h2which is truncated to a ﬁnite sum as≈∑1Γ(1 − α) K ≤k≤K
0

ϕ(y ) + 1
( k )α1

1−ϕ(yk )2

ϕ′ (yk )h2with some integers K0 and K1 .
As |k| → ∞, ϕ′ (yk ) converges to the zero more rapidly than growth of the integrandcaused by the singularity, and their products and accumulations are negligible for suﬃciently large |k|. Therefore choice of h, K0 , and K1 is crucial for accuracy of numericalintegration. Table 1 and Table 2 show numerical evaluation of∑1Ih =Γ(1 − α) K ≤k≤K
0

ϕ(y ) + 1
( k )α1−ϕ(yk )2

1

ϕ′ (yk )h2with α = 0.5, 0.9 respectively. In the tables, we show discretization and truncation
parameters as L = K0 h, U = K1 h, and K = K1 − K0 . Values of the exact solution[ ]α 2
c D0 t=t=12
Γ(3 − α)are also shown in tables.
From Table 1, both double precision arithmetic and 100 decimal digit computation
enable us to ﬁnd approximate values with (L, U ) = (−5, 3.16) for α = 0.5. On theother hand, for α = 0.9, parameters (L, U ) = (−5, 5) and K ≥ 16 in 100 decimal digit
return approximations while the double precision arithmetic returns inﬁnity if U is greater

175
Table 1: Numerical integration Ih with α = 0.5.
LU

K

double precision

100 digits−5.00 3.16

8
16
32

1.50652114572
1.50451382597
1.50450554499

1.50652118788
1.50451384705
1.50450555553−5.00 3.17

8
16
32

+∞
+∞
+∞

1.50798169193
1.50451519723
1.50450555566

exact value

1.50450555613 . . .

Table 2: Numerical integration Ih with α = 0.9.
LU

K

double precision

100 digits−5.00 3.16

8
16
32

1.87072828925
1.87900443987
1.87188368373

1.90636895125
1.89682476417
1.88079405669−5.00 3.17

8
16
32

+∞
+∞
+∞

1.90645242453
1.89759689138
1.88210055004−5.00 5.00

8
16
32

+∞
+∞
+∞

1.90117136452
1.91113179851
1.91115819291

exact value

1.91115819293 . . .

than or equal to 3.17. This means that double precision arithmetic is not enough fornumerical evaluation of fractional derivatives without special treatments of singularities
in Ih , especially when the power −α is close to −1. In our numerical computation,
exflib [11] is used for multiple-precision arithmetic.
Figure 1 shows proﬁles of fractional derivatives c D0α u(t) for u(t) = t2 obtained bynumerically integrating (1). Figure 1(a) and (b) are processed with double precision
arithmetic and 100 decimal digits arithmetic respectively. We set (L, U ) = (−5, 3.16) indouble precision arithmetic, and (L, U ) = (−5, 5) in 100 decimal digits arithmetic, with
suﬃciently large K. Obviously, results for α = 0.9 do not coincide. This also shows that
double precision is not enough for reliable computations. More detail discussions andother numerical examples will be found in [12].

176
2

2

2

u(t)=t2
u'(t)=2talpha=0.2
alpha=0.4
alpha=0.6
alpha=0.8
alpha=0.9u(t)=tu’(t)=2talpha=0.2
alpha=0.4
alpha=0.6
alpha=0.8
alpha=0.9
1.5

1.5

1

1

0.5

0.5

0

0
0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

t

t

(a) Results by double precision(b) Results by 100 decimal digits

Figure 1: Fractional derivatives c D0α u(t) of u(t) = t2 .

4

Numerical Results for Fractional Diﬀerential Equations

In this section we shall demonstrate some examples which shows eﬃciency of combinationof the Chebyshev spectral method and multiple-precision arithmetic to achieve arbitrary
accurate numerical solutions to fractional order diﬀerential equations.
Example 1.
u(1/3) (t) = f (t),u(0) = 0.

0 < t < 3,

5!
10 !
15 !
Three cases are considered: f (t) = Γ(17/3)
t14/3 , Γ(32/3)
t29/3 , and Γ(47/3)
t44/3 . The exact
solutions are u(t) = t5 , t10 , and t15 respectively. Since each solution is a monomial, it isexactly treated by the Chebyshev spectral method, and hence only rounding errors appear
as numerical errors.
In Figure 2, the horizontal axis is the spectral degree N and the vertical axis shows
errors in numerical results:
max |u(xj ) − uN (xj )|0≤j≤Nin the logarithmic scale. In the ﬁgure, purple, green, and blue graph represent numericalerrors for u(t) = t5 , t10 , and t15 respectively. In each color, plus (+) sign and cross (×) areerror in 50 decimal digits and 100 decimal digits respectively. In each setting, behaviourof errors suddenly changes at a certain order N0 . More precisely, when N becomes larger
than or equal to N0 , the numerical errors are less than 10−50 or 10−100 which correspond to
the rounding errors in each computation. This indicates that the solution is a polynomialof degree almost N0 .

177
1e+40
f(x)=x**4.666 100 digit
f(x)=x**4.666 50 digit
f(x)=x**9.666 100 digit
f(x)=x**9.666 50 digit
f(x)=x**14.666 100 digit
f(x)=x**14.666 50 digit

1e+20

Maximum Error (logarithmic scale)

1e+00

1e-20

1e-40

1e-60

1e-80

1e-1001e-1200

5

10
N : Spectral Degree

15

20

Figure 2: Error decay for Example 1, the exact solutions are monomials.
Example 2. [9]u(5/2) (t) − 3u(2/3) (t) = f (t),0 < t ≤ 1,′′′
u(0) = 1, u (0) = β, u (0) = β 2 .

1e+20

1e+20

1e+00

1e+00

Maximum Error (Logarithmic Scale)

Maximum Error (Logarithmic Scale)

We set the exact solution u(t) = eβt , and calculate f (t) as the left-hand side of theequation by using (1). Figure 3 shows numerical results by double precision (blue), 50
decimal digits (green), and 100 decimal digits (purple). The error decreases exponentially
until the level of rounding error is almost reached. This behavior of errors is expectedfrom the property of the spectral methods.

1e-20

1e-40

1e-60

1e-80

1e-100

1e-20

1e-40

1e-60

1e-80

1e-100
gamma=1, 100 digits
gamma=1, 50 digits
gamma=1, double

gamma=6, 100 digits
gamma=6, 50 digits
gamma=6, double

1e-1201e-120
0

10

20

30

40
Spectral Degree N(a) β = 1

50

60

70

80

0

10

20

30

40

50

Spectral Degree N(b) β = 6

Figure 3: Error decay for Example 2, the exact solution is eβt .

60

70

80

178
Example 3. [9]u′′ (t) + u(1.5) (t) + u(t) = f (t),u(0) = 0, u′ (0) = ω.

0 < t ≤ 1,

1e+20

1e+20

1e+00

1e+00

Maximum Error (Logarithmic Scale)

Maximum Error (Logarithmic Scale)

The equation is known as the Bagley-Torvik equation. We set the exact solution u(t) =sin ωt. Then f (t) is calculated by using (1). Numerical results in several precisions areshown in Figure 4. The error decreases similarly in the preceding example, since the exact
solution is also analytic in this example.

1e-20

1e-40

1e-60

1e-80

1e-100

1e-20

1e-40

1e-60

1e-80

1e-100alpha=1, 100 digits
alpha=1, 50 digits
alpha=1, doublealpha=4 pi, 100 digits
alpha=4 pi, 50 digits
alpha=4 pi, double

1e-1201e-120
0

10

20

30

40
Spectral Degree N(a) ω = 1

50

60

70

80

0

10

20

30

40

50

60

70

80

Spectral Degree N(b) ω = 4π

Figure 4: Error decay for Example 3, the exact solution is sin ωt.

Example 4. In the same equation and initial condition as Example 1, we adopt f (t) =Γ(7/3)t, and the exact solution is u(t) = t4/3 . Figure 5 shows numerical results, where
both the horizontal and vertical axes are in the logarithmic scale. In this setting numericalerrors decay as O(N −2.66 ), which is not exponential decay but polynomial decay withrespect to N . This suggests that the exact solution is not analytic.

5

Conclusion

In the paper, by using the Chebyshev spectral collocation method a concrete expressionof coeﬃcients in the linear combination of basis functions is derived. It is used for solving
one-dimensional fractional diﬀerential equations numerically. We also investigate analyticity of solutions numerically. To do so, in numerical computation the multiple-precision
arithmetic is used to estimate the inﬂuence of rounding error. Numerical results are verysatisfactory. Spectral accuracy is seen and it shows analyticity of solutions. Moreover,when exact solutions are monomials, numerical results are succeeded to determine thedegree of monomials.

179
1e+01

Maximum Error (Logarithmic Scale)

1e+00

1e-01

1e-02

1e-03

1e-04

1e-05
O(N**(-2.66))
D**(1/3)u=t, 100 digits
D**(1/3)u=t, double
1e-06
1

10
Spectral Degree N (logarithmic scale)100

Figure 5: Error decay for Example 4 in double precision arithmetic (Blue) and 100 decimal
digits (Green), the exact solution is t4/3 .

References[1] D. Baleanu et al., Fractional Calculus: Models and Numerical Methods, World
Scientiﬁc, Singapore, 2012.
[2] C. Canuto et al., Spectral Methods: Fundamentals in Single Domains, Springer,
Berlin, 2006.
[3] C. Canuto et al., Spectral Methods: Evolution to Complex Geometries and Applications to Fluid Dynamics, Springer, Berlin, 2007.
[4] R. Caponetto et al., Fractional Order Systems: Modeling and Control Applications,
World Scientiﬁc, Singapore, 2010.
[5] M. Caputo, Linear Model of Dissipation Whose Q is Almost Frequency Independent— II, Geophys J. R. Astron. Soc., 13(1967), pp. 529–539.
[6] C. Cattani et al., Fractional Dynamics, De Gruyter Open, Berlin, 2015.
[7] A. Dabiri and E. A. Butcher, Eﬃcient Modiﬁed Chebyshev Diﬀerentiation Matrices
for Fractional Diﬀerential Equations, Commun. Nonlinear Sci. Numer. Simulat.,50(2017), pp. 284–310.
[8] K. C. Datta, H. Imai and H. Sakaguchi, On numerical Challenge for the Distinction between Smooth Functions and Analytic Functions, Theo. Appl. Mech. Jpn.,
62(2013), pp. 107–117.
[9] E. H. Doha, A. H. Bhrawy and S. S. Ezz-Eldien, Eﬃcient Chebyshev Spectral
Methods for Solving Multi-term Fractional Orders Diﬀerential Equations, Appl.
Math. Model., 35(2011), pp. 5662–5672.

180[10] H. Fujiwara et al., High-precision Numerical Computation of Integral Equations of
the First Kind, Transactions of JSIAM, 15(3) (2005), pp. 419–434 (in Japanese).
[11] H. Fujiwara and Y. Iso, Design of a Multiple-precision Arithmetic Package for a
64-bit Computing Environment and Its Application to Numerical Computation of
Ill-posed Problems, IPSJ J., 44(3) (2003), pp. 925–931.
[12] H. Fujiwara, Multiple-precision Arithmetic Environment in MATLAB and Its Application to Reliable Computation of Fractional Order Derivatives, Math for Industry
Kenkyu, 8 (2018) pp. 83–120[13] H. Imai et al., A Direct Approach to an Inverse Problem, GAKUTO Int. Ser. Math.
Sci. Appl., 12 (1999), pp. 223–232.
[14] H. Imai, S. Sakaguchi and Y. Iso, On Numerical Computation of the Toricomi
Equation, Theo. Appl. Mech. Jpn., 59 (2010), pp. 359–372.
[15] H. Imai and T. Takeuchi, Some Advanced Applications of the Spectral Collocation
Method, GAKUTO Int. Ser. Math. Sci. Appl., 17 (2001), pp. 323–335.
[16] H. Imai, T. Takeuchi and M. Kushida, On Numerical Simulation of Partial Diﬀerential Equations in Inﬁnite Precision, Adv. Math. Sci. Appl., 9(2) (1999), pp. 1007–
1016.
[17] R. L. Magin, Fractional Calculus in Bioengineering, Begell House, Connecticut,
2006.
[18] H. Takahashi and M. Mori, Double Exponential Formulas for Numerical Integration,
Publ. Res. Inst. Math. Sci., Kyoto Univ., 9 (1974), pp. 721–741.
[19] T. Takeuchi and H. Imai, Direct Numerical Simulations of Cauchy Problems for the
Laplace Operator, Adv. Math. Sci. Appl., 13(2) (2003), pp. 587–609.
[20] J. A. Tenreiro Machado, V. Kiryakova, and F. Mainardi, A Poster about the Recent
History of Fractional Calculus, Fract. Calc. Appl. Anal., 13(3)(2010), pp. 329–334.
[21] W. Y. Tian, W. Deng and Y. Wu, Polynomial Spectral Collocation Method forSpace Fractional Advection Diﬀusion Equation, Numer. Methods Partial Diﬀer.
Equ., 30(2)(2014), pp. 514–535.