On the numerical solution of linear fractional differential equations

Introduction

The numerical analysis of fractional differential equations and operators has been widely developed in the last decades, see, for example,^{1,2,4–17,19,20}. The aim of the present paper is to report some results regarding collocation-type approximations for solving initial value problems for linear Caputo-type fractional differential equations. For simplicity, we restrict ourselves to differential equations with a single fractional derivative of order α ∈ (0, 1).

Let R = ( − ∞ , ∞ ) , N = { 1,2 , … } . Let C ^k (Ω) be the set of all k times (k ≥ 0) continuously differentiable functions on Ω ⊂ R , C⁰(Ω) = C(Ω). In particular, by C[0, b] we denote the Banach space of continuous functions u on [0, b] with the typical norm ‖u‖ _∞ ≔ max_0≤t≤b|u(t)|. Let α ∈ (0, 1). We consider the following initial value problem:

( D Cap α y ) ( t ) + a ( t ) y ( t ) = f ( t ) , 0 ≤ t ≤ b , b > 0 ,

(1)

y ( 0 ) = y 0 ,

(2)

( D Cap α y ) ( t ) = 1 Γ ( 1 − α ) d d t ∫ 0 t ( t − s ) − α ( y ( s ) − y ( 0 ) ) d s , 0 < t ≤ b .

Here by Γ we denote the Euler gamma function: Γ ( x ) = ∫ 0 ∞ e − s s x − 1 d s , x > 0. Note that for any solution y ∈ C[0, b] of (1)–(2) we have that D Cap α y ∈ C [ 0 , b ] . To study the existence, uniqueness and regularity of the exact solution y of (1)–(2), we introduce the following class of weighted functions C^m,ν(0, b], an adaptation of a more general class of functions first studied by Vainikko in ¹⁸ .

For given b ∈ R , b > 0 , m ∈ N and 0 < ν < 1, by C^m,ν(0, b] we denote the set of continuous functions u : [ 0 , b ] → R which are m times continuously differentiable in (0, b] such that

| u ( i ) ( t ) | ≤ c t 1 − ν − i , t ∈ ( 0 , b ] , i = 1 , … , m ,

‖ u ‖ C m , ν ( 0 , b ] ≔ ‖ u ‖ ∞ + ∑ i = 1 m sup 0 < t ≤ b t i − 1 + ν u ( i ) ( t ) , u ∈ C m , ν ( 0 , b ] ,

C q [ 0 , b ] ⊂ C q , ν ( 0 , b ] ⊂ C m , μ ( 0 , b ] ⊂ C [ 0 , b ] , q ≥ m ≥ 1 , 0 < ν ≤ μ < 1 .

We have the following existence, uniqueness and regularity result, see, for example,^14,19.

Theorem 1

Let α ∈ (0, 1) and a, f ∈ C[0, b]. Then (1)–(2) has a unique solution y ∈ C[0, b]. Moreover, if a, f ∈ C^q,μ(0, b], q ∈ N , 0 < μ < 1, then y ∈ C^q,ν(0, b], where

ν = max { 1 − α , μ } .

Remark 1

Note that the statement y ∈ C^q,ν(0, b] for the solution y of (1)–(2) in Theorem 1 cannot, in general, be improved. For example, let us consider the following initial value problem:

D Cap 1 2 y ( t ) = π 2 , t ∈ [ 0,1 ] ; y ( 0 ) = 1 .

The exact solution to this problem is y ( t ) = t + 1 , t ∈ [0, 1]. It is clear that y ∈ C^q,1/2(0, 1] for all q ∈ N , that is, the derivatives of y(t) are unbounded at t = 0.

It is easy to see that problem (1)–(2) can be reformulated as an equivalent Volterra integral equation for y:

y ( t ) = y 0 + 1 Γ ( α ) ∫ 0 t ( t − s ) α − 1 [ f ( s ) − a ( s ) y ( s ) ] d s , 0 ≤ t ≤ b ,

y = T y + g ,

(3)

( T y ) ( t ) = − 1 Γ ( α ) ∫ 0 t ( t − s ) α − 1 a ( s ) y ( s ) d s , 0 ≤ t ≤ b ,

(4)

g ( t ) = y 0 + 1 Γ ( α ) ∫ 0 t ( t − s ) α − 1 f ( s ) d s , 0 ≤ t ≤ b .

(5)

Thus, to find a numerical solution to (1)–(2), it is sufficient to construct an approximation method for the integral equation (3).

Collocation method

In order to take into account the potential non-smoothness of the exact solution y = y(t) of (3) at the origin t = 0, we introduce on the interval [0, b] a graded grid Π _N = {t₀, …, t _N }. More exactly, for given N ∈ N , we divide the underlying interval of integration [0, b] into N subintervals [t_j−1, t _j ], j = 1, …, N, with grid points

t j = b j N r , j = 0 , … , N .

(6)

The real number r ≥ 1 characterizes the nonuniformity of the grid. If r = 1, then the grid is uniform; if r > 1, then the grid points are more densely clustered near the left endpoint 0.

Next, for a given integer k ≥ 0, we introduce

S k ( − 1 ) ( Π N ) ≔ v : v | [ t j − 1 , t j ] ∈ π k , j = 1 , … , N .

Here v | [ t j − 1 , t j ] is the restriction of v : [ 0 , b ] → R onto the subinterval [t_j−1, t _j ] ⊂ [0, b] and π _k denotes the set of polynomials of degree not exceeding k. Note that the elements of S k ( − 1 ) ( Π N ) may have jump discontinuities at the interior points t₁, …, t_N−1 of the grid Π _N . We choose m ∈ N points η₁, …, η _m in the interval [0, 1] so that

0 ≤ η 1 < … < η m ≤ 1

(7)

t j p = t j − 1 + η p ( t j − t j − 1 ) , p = 1 , … , m , j = 1 , … , N .

(8)

Approximations y N ∈ S m − 1 ( − 1 ) ( Π N ) to the solution y of (3) can be found by assuming that the following conditions hold:

y N ( t j p ) = ( T y N ) ( t j p ) + g ( t j p ) , p = 1 , … , m , j = 1 , … , N ,

(9)

y N = P N T y N + P N g ,

(10)

( P N v ) ( t j p ) = v ( t j p ) , p = 1 , … , m , j = 1 , … , N , v ∈ C [ 0 , b ] .

If η₁ = 0, then by (P _N v)(t_j1) we denote the right limit lim t → t j − 1 + ( P N v ) ( t ) . If η _m = 1, then by (P _N v)(t _jm ) we denote the left limit lim t → t j − ( P N v ) ( t ) .

The collocation conditions (9) form a system of equations whose exact form is determined by the choice of a basis in the space S m − 1 ( − 1 ) ( Π N ) . We will use the Lagrange fundamental polynomial representation:

y N ( t ) = ∑ l = 1 N ∑ k = 1 m c l k φ l k ( t ) , t ∈ [ 0 , b ] ,

(11)

φ l k ( t ) = 0 f o r t ∉ [ t l − 1 , t l ] , ∏ i = 1 , i ≠ k m t − t l i t l k − t l i f o r t ∈ [ t l − 1 , t l ] , k = 1 , … , m , l = 1 , … , N .

As usual, we take ∏ i = 1 , i ≠ k m t − t l i / t l k − t l i = 1 if m = 1. Note that y _N (t _lk ) = c _lk , k = 1, …, m, l = 1, …, N, and using representation (11) the collocation conditions (9) form the following system of linear algebraic equations with respect to the unknown values {c _jp }:

c j p = ∑ l = 1 N ∑ k = 1 m ( T φ l k ) ( t j p ) c l k + g ( t j p ) , p = 1 , … , m , j = 1 , … , N .

After solving this system of equations, we have found the approximation y _N to y, the solution of (1)–(2). It follows from ¹⁹ that the following convergence result holds (cf.^3,11,14).

Theorem 2

(i) Let α ∈ (0, 1) and a, f ∈ C[0, b]. Let N , m ∈ N and assume that the grid points (6) and the collocation points (8) are used, with collocation parameters η₁, …, η _m satisfying (7).

Then there exists an integer N₀ such that for all N ≥ N₀ equation (10) possesses a unique solution y N ∈ S m − 1 ( − 1 ) ( Π N ) and

sup 0 ≤ t ≤ b | y ( t ) − y N ( t ) | → 0 a s N → ∞ ,

(ii) If, in addition to (i), we assume that a, f ∈ C^m,μ(0, b], where 0 < μ < 1, then for all N ≥ N₀ the error estimate

sup 0 ≤ t ≤ b | y ( t ) − y N ( t ) | ≤ c N − r ( 1 − ν ) f o r 1 ≤ r < m 1 − ν , N − m f o r r ≥ m 1 − ν ,

(12)

Iterated numerical method

Theorem 2 shows that if the functions a and f are smooth enough, and if the grading parameter r is chosen to be sufficiently large, the expected convergence order of the proposed numerical method is of order O(N^−m). This convergence order can be improved by using an iterated method, assuming a little more smoothness of a and f, together with a more precise choice of collocation parameters η₁, …, η _m . More precisely, let y _N be the solution of (10), with N ≥ N₀ (see Theorem 2). Denote

y N i t ≔ T y N + g , N ≥ N 0 ,

(13)

Theorem 3

Let α ∈ (0, 1) and m ∈ N . Assume that a ∈ C^m+1[0, b] and f ∈ C^m+1,μ(0, b], where 0 < μ < 1. Let N ∈ N and assume that the grid points (6) and the collocation parameters η₁, …, η _m satisfying (7) are used. Moreover, assume that the collocation parameters are chosen so that the quadrature approximation

∫ 0 1 F ( x ) d x ≈ ∑ k = 1 m w k F ( η k )

(14)

Then there exists N 0 ∈ N so that for N ≥ N₀ the approximation y N i t defined by (13) is unique and the error estimate

‖ y − y N i t ‖ ∞ ≤ c N − r ( 1 + α − ν ) f o r 1 ≤ r < m + α 1 + α − ν , N − m − α f o r r ≥ m + α 1 + α − ν ,

(15)

Numerical example

Consider the following initial value problem:

D Cap 1 2 y ( t ) + ( 1 − t ) y ( t ) = t 1 2 − t 3 2 − 2 t + 2 + Γ 3 2 , 0 ≤ t ≤ 1 ; y ( 0 ) = 2 .

(16)

We see that (16) is a special problem of (1)–(2) with

α = 1 2 , a ( t ) = 1 − t , f ( t ) = t 1 2 − t 3 2 − 2 t + 2 + Γ 3 2 , y 0 = 2 , b = 1 .

Note that a ∈ C ^m [0, 1] and f ∈ C^m,1/2(0, 1] for all m ∈ N . Therefore by Theorem 1 it follows that y ∈ C^m,ν(0, b], where

ν = max 1 − 1 2 , 1 2 = 1 2 .

(17)

Special case of collocation parameters

We apply the collocation method (see Section 2) and its iterated version (see Section 3) with m = 2, using as collocation parameters the shifted Gauss-Legendre points

η 1 = 3 − 3 6 , η 2 = 1 − η 1 ,

which satisfy the conditions set for collocation parameters in Theorem 3. It follows from Theorem 2 that for sufficiently large N ∈ N we have

sup 0 ≤ t ≤ b | y ( t ) − y N ( t ) | ≤ c 0 N − 0.5 r i f 1 ≤ r < 4 , N − 2 i f r ≥ 4 ,

(18)

where

y ( t ) = t 1 2 + 2 , 0 ≤ t ≤ 1 ,

is the exact solution of problem (16) and c₀ is a positive constant independent of N. Similarly, it follows from Theorem 3 that

‖ y − y N i t ‖ ∞ ≤ c 1 N − r i f 1 ≤ r < 5 2 , N − 2.5 i f r ≥ 5 2 ,

(19)

where c₁ is a positive constant independent of N.

In Table 1 some results of numerical experiments for different values of the parameter r are presented. The actual numerical error ɛ _N for sup_0≤t≤b|y(t) − y _N (t)| and ε N i t for ‖ y − y N i t ‖ ∞ are calculated as follows:

ε N ≔ max j = 1 , … , N max k = 0 , … , 10 | y ( τ j k ) − y N ( τ j k ) | , ε N i t ≔ max j = 1 , … , N max k = 0 , … , 10 | y ( τ j k ) − y N i t ( τ j k ) | ,

where

τ j k ≔ t j − 1 + k ( t j − t j − 1 ) / 10 , k = 0 , … , 10 , j = 1 , … , N ,

with t _j defined by (6). The ratios

ϱ N ≔ ε N / 2 ε N , ϱ N i t ≔ ε N / 2 i t ε N i t ,

characterizing the observed convergence rates, are also presented. Due to (18), the ratios ϱ _N for the non-iterated collocation method for r = 1, r = 2 and r = 4 ought to be approximately 2^0.5 ≈ 1.414, 2¹ = 2 and 2² = 4, respectively. Due to (19), the ratios ϱ N i t for the iterated method for r = 1, r = 2 and r = 4 ought to be approximately 2¹ = 2, 2² = 4 and 2^2.5 = 5.667, respectively. All these ratios are also given in the corresponding rows of Table 1. As we can see from Table 1 the numerical results are in good accord with the theoretical estimates given by Theorem 2 and Theorem 3 (the estimates (12) and (15)).

OPEN IN VIEWER

Table 1.

Numerical results for shifted Gauss points.

r = 1
Collocation method			Iterated method
N	ɛ N	ϱ N	N	ε N i t	ϱ N i t
4	1.44702 × 10−1		4	8.08994 × 10−3
8	1.03464 × 10−1	1.399	8	4.19988 × 10−3	1.926
16	7.38164 × 10−2	1.402	16	2.15907 × 10−3	1.945
32	5.25566 × 10−2	1.405	32	1.10167 × 10−3	1.960
64	3.73555 × 10−2	1.407	64	5.58988 × 10−4	1.971
128	2.65149 × 10−2	1.409	128	2.82463 × 10−4	1.979
256	1.88008 × 10−2	1.410	256	1.42303 × 10−4	1.985
		≈ 1.414			2
r = 2
4	7.38164 × 10−2		4	2.15907 × 10−3
8	3.73555 × 10−2	1.976	8	5.58988 × 10−4	3.863
16	1.88008 × 10−2	1.987	16	1.42303 × 10−4	3.928
32	9.43259 × 10−3	1.993	32	3.59053 × 10−5	3.963
64	4.72453 × 10−3	1.997	64	9.01817 × 10−6	3.981
128	2.36438 × 10−3	1.998	128	2.24132 × 10−6	4.023
256	1.18270 × 10−3	1.999	256	5.60982 × 10−7	3.995
		2			4
r = 4
4	2.56513 × 10−2		4	2.74746 × 10−3
8	6.80136 × 10−3	3.772	8	3.76009 × 10−4	7.307
16	1.86347 × 10−3	3.650	16	7.91378 × 10−5	4.751
32	4.66852 × 10−4	3.992	32	1.39477 × 10−5	5.674
64	1.16773 × 10−4	3.998	64	2.46083 × 10−6	5.668
128	2.91971 × 10−5	4.000	128	4.36171 × 10−7	5.642
256	7.29954 × 10−6	4.000	256	7.72673 × 10−8	5.645
		4			≈ 5.657

General case of collocation parameters

We also briefly highlight the case when the chosen collocation parameters do not satisfy the quadrature condition set in Theorem 3, by applying both methods (for m = 2) with collocation parameters

η 1 = 0.1 , η 2 = 0.9 .

The obtained results are given in Table 2. We see that in this case the iterated method does not attain for r ≥ 5/2 the convergence order O(N^−2.5) predicted by Theorem 3. This shows that the collocation parameter assumption of Theorem 3 cannot be relaxed.

OPEN IN VIEWER

Table 2.

Numerical results for non-Gauss points.

Collocation method			Iterated method
N	ɛ N	ϱ N	N	ε N i t	ϱ N i t
r = 1
4	1.14326 × 10−1		4	1.03667 × 10−2
8	8.15843 × 10−2	1.401	8	5.54147 × 10−3	1.871
16	5.81064 × 10−2	1.404	16	2.86848 × 10−3	1.932
32	4.13139 × 10−2	1.407	32	1.46152 × 10−3	1.963
64	2.93332 × 10−2	1.409	64	7.38548 × 10−4	1.980
128	2.08040 × 10−2	1.410	128	3.71550 × 10−4	1.988
256	1.47427 × 10−2	1.411	256	1.86458 × 10−4	1.993
		≈ 1.414			2
r = 2
4	5.81064 × 10−2		4	3.18598 × 10−3
8	2.93332 × 10−2	1.981	8	9.50767 × 10−4	3.351
16	1.47427 × 10−2	1.990	16	2.54532 × 10−4	3.735
32	7.39120 × 10−3	1.995	32	6.64282 × 10−5	3.832
64	3.70066 × 10−3	1.997	64	1.70902 × 10−5	3.887
128	1.85163 × 10−3	1.999	128	4.35129 × 10−6	3.928
256	9.26126 × 10−4	1.999	256	1.09932 × 10−6	3.959
		2			4
r = 4
4	1.65300 × 10−2		4	4.68951 × 10−3
8	4.35499 × 10−3	3.796	8	1.11609 × 10−3	4.202
16	1.18024 × 10−3	3.690	16	2.84909 × 10−4	3.917
32	3.00346 × 10−4	3.930	32	6.61433 × 10−5	4.307
64	7.65563 × 10−5	3.923	64	1.55383 × 10−5	4.257
128	1.93304 × 10−5	3.960	128	3.76263 × 10−6	4.130
256	4.85373 × 10−6	3.983	256	9.14624 × 10−7	4.114
		4			≈ 5.657

Statements and declarations