Smoothing transformation and spline collocation for nonlinear fractional initial and boundary value problems

Accepted Manuscript Smoothing transformation and spline collocation for nonlinear fractional initial and boundary value problems Arvet Pedas, Enn Tamme, Mikk Vikerpuur PII: DOI: Reference:

S0377-0427(16)30555-6 http://dx.doi.org/10.1016/j.cam.2016.11.022 CAM 10895

To appear in:

Journal of Computational and Applied Mathematics

Received date: 3 June 2016 Revised date: 8 November 2016 Please cite this article as: A. Pedas, E. Tamme, M. Vikerpuur, Smoothing transformation and spline collocation for nonlinear fractional initial and boundary value problems, Journal of Computational and Applied Mathematics (2016), http://dx.doi.org/10.1016/j.cam.2016.11.022 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Smoothing transformation and spline collocation for nonlinear fractional initial and boundary value problemsI Arvet Pedas, Enn Tamme∗, Mikk Vikerpuur Institute of Mathematics and Statistics, University of Tartu Liivi 2, 50409 Tartu, Estonia

Abstract We construct and justify a class of high order methods for the numerical solution of initial and boundary value problems for nonlinear fractional differential equations of the form (D∗α y)(t) = f (t, y(t)) with Caputo type fractional derivatives D∗α y of order α > 0. Using an integral equation reformulation of the underlying problem we first regularize the solution by a suitable smoothing transformation. After that we solve the transformed equation by a piecewise polynomial collocation method on a mildly graded or uniform grid. Optimal global convergence estimates are derived and a superconvergence result for a special choice of collocation parameters is established. To illustrate the reliability of the proposed algorithms two numerical examples are given. Keywords: Nonlinear fractional boundary value problem; Caputo derivative; Weakly singular integral equation; Smoothing transformation; Spline collocation method PACS: 65L60, 65L20, 26A33

1. Introduction In the present paper we study the convergence behavior of a modified spline collocation method for the numerical solution of fractional initial and boundary value problems of the form (D∗α y)(t) = f (t, y(t)), 0 ≤ t ≤ b, (1.1) n0 X

(j)

βij0 y (0) +

j=0

n1 l X X k=1 j=0

βijk y (j) (bk ) = γi ,

i = 0, . . . , n − 1 , n := dαe,

(1.2)

This work was supported by Estonian Science Foundation Grant No. 9104 and by institutional research funding IUT20–57 of the Estonian Ministry of Education and Research. ∗ Corresponding author Email addresses: [email protected] (Arvet Pedas), [email protected] (Enn Tamme), [email protected] (Mikk Vikerpuur) I

Preprint submitted to Elsevier

November 22, 2016

where βij0 , βijk , γi ∈ R := (−∞, ∞), dαe is the smallest integer greater or equal to α ∈ R, n − 1 < α < n , 0 < b1 < · · · < bl ≤ b, l, n ∈ N := {1, 2, . . .}, n0 , n1 ∈ N0 := {0} ∪ N, n0 < n, n1 < n,

(1.3)

f : [0, b] × R → R is a given continuous function, and D∗α y is the Caputo derivative of order α > 0 of an unknown function y. Observe that a special case of problem (1.1)–(1.2) is the following initial value problem: (D∗α y)(t) = f (t, y(t)), 0 ≤ t ≤ b, α > 0, y (i) (0) = γi , i = 0, . . . , n − 1 , n = dαe.

(1.4)

Thus, the results of this paper can be applied also to initial value problems of the form (1.4). The Caputo differentiation operator D∗α of order α > 0 can be defined by the formula (see, e.g, [1]) (D∗α y)(t) := (Dα (y − Qn−1 [y]))(t), t > 0, n = dαe, where Qn−1 [y](s) :=

n−1 (i) X y (0)

i!

i=0

si

and Dα is the Riemann–Liouville fractional differentiation operator of order α > 0 : dn n−α (D y)(t) := n (J y)(t), dt α

t > 0,

n = dαe.

Here J α , the Riemann–Liouville integral operator of order α > 0, is defined by the formula Z t 1 α (J y)(t) := (t − s)α−1 y(s) ds, t > 0, (1.5) Γ(α) 0 where Γ is the Euler gamma function. For α = 0 we set D0 = D∗0 = J 0 := I where I is the identity mapping. If α = n ∈ N then Dn y = D∗n y = y (n) where y (n) is the usual n-th order derivative of y. The integral operator J α , α > 0, is linear, bounded and compact as an operator from L∞ (0, b) into C[0, b] (see, e.g., [2]). Moreover, we have for any y ∈ L∞ (0, b) that (see, e.g., [3]) (J α y)(k) ∈ C[0, b], (J α y)(k) (0) = 0, α > 0, k = 0, . . . , dαe − 1, (1.6) J α J β y = J α+β y,

α > 0,

Dβ J α y = D∗β J α y = J α−β y,

β > 0, 0 < β ≤ α.

(1.7) (1.8)

Differential equations involving differential operators of fractional (non-integer) order have been proved to be a valuable tool in modeling many phenomena in the fields of physics, mechanics, chemistry, engineering, finance and others (see, for example, [4–6]). Mathematical aspects of fractional differential equations and methods of their solution have been discussed 2

by many authors (see [1, 3, 7] and the references cited therein). Various existence and uniqueness results for initial value problems (1.4) can be found in [1] and for some boundary value problems of the form (1.1)–(1.2) in [8]. For most fractional differential equations we cannot provide methods to compute their solutions analytically and it is necessary to use numerical methods. Efficient numerical methods for solving integral and differential equations are spline collocation methods. These methods are also applicable for the numerical solution of fractional differential equations since initial and boundary value problems for fractional differential differential equations can be reformulated as certain integral equations of the second kind with weakly singular kernels (see, e.g., [1]). The convergence of spline collocation methods for weakly singular integral equations is analysed in [9–11] and for fractional differential equations in [12–18]. A review of spectral collocation methods for solving fractional differential equations is given in [19]. From Theorem 2.1 below we see that the derivatives of the solution of problem (1.1)–(1.2) and the corresponding integral equation may be unbounded only near the initial point 0 of the interval of integration [0, b]. Due to the lack of smoothness of the exact solution, piecewise polynomial collocation methods based on uniform grid for solving such equations show slow convergence behavior [20]. In order to construct methods with higher convergence order it is necessary to take into account the possible singular behavior of the exact solution. In particular, we can construct high order collocation methods by using polynomial splines on special graded grids where the grid points are more densely clustered near the singular point t = 0 of the exact solution y(t) of (1.1)–(1.2) [13–16]. However, the use of strongly non-uniform grids may cause an accumulation of rounding errors and unstable behavior of numerical results. To diminish loss of precision we may resort to collocation on the uniform grid using non-polynomial basic functions which reflect the singular behavior of the exact solution. This approach has been applied for solving Volterra integral equations in [21–23] and for solving fractional differential equations in [24–26]. In the present paper we use an alternative approach for diminishing loss of precision. First we introduce an integral equation reformulation of problem (1.1)–(1.2). Then we perform in the integral equation a smoothing change of variables t = b1−ρ τ ρ (ρ ∈ [1, ∞)) so that the singularities of the (usual) derivatives of the exact solution will be milder or disappear. After that we solve the transformed integral equation by a piecewise polynomial collocation method on a mildly graded or uniform grid. The final step of our method is based on a conversion of the obtained spline functions into the (typically non-polynomial) approximate solution for (1.1)–(1.2). Similar ideas for solving Volterra integral equations are used in [27–31] and for solving fractional differential equations in [32–35]. The rest of our paper is arranged as follows. In Section 2 some a priori estimates for higher order derivatives of the exact solution of problem (1.1)–(1.2) are presented (see Theorem 2.1). These estimates will play a key role in the convergence analysis of the proposed algorithms in Sections 5 and 6. In Section 3 some results about piecewise polynomial interpolation on graded grids are presented. In Section 4 a class of smoothing transformations is introduced and based on these transformations and spline collocation techniques a class of numerical methods for solving (1.1)–(1.2) is constructed. In Sections 5 and 6 the attainable order 3

of convergence of proposed algorithms is studied. Optimal global convergence estimates are derived and a superconvergence result for a special choice of collocation parameters is established. Finally, in Section 7 the obtained theoretical results are tested by two numerical examples. The main results of the present paper are formulated in Theorems 2.1, 5.1 and 6.1. 2. Smoothness of the solution Using some ideas from [14] (see also [1, 16]) we find first an integral equation reformulation for the problem (1.1)–(1.2). Let y ∈ C[0, b] be such that D∗α y ∈ C[0, b]. Introduce a new unknown function z := D∗α y. Then (see [1, 3]) y(t) = (J α z)(t) +

n−1 X

cλ t λ ,

t ∈ [0, b], ,

λ=0

n = dαe ∈ N,

(2.1)

where cλ ∈ R (λ = 0, . . . , n − 1) are arbitrary constants. The function y of the form (2.1) satisfies the boundary value conditions (1.2) if and only if (see (1.8) and (1.6)) n0 X j=0

βij0 j! cj +

n1 l X X



βijk (J

k=1 j=0

α−j

z)(bk ) +

n−1 X λ=j

 λ! λ−j b cλ = γi , (λ − j)! k

i = 0, . . . , n − 1.

We rewrite these conditions in the form  j n1 n−1  l X l X X X X j! j−λ βiλk bk cj = γi − j! βij0 + βijk (J α−j z)(bk ), (j − λ)! j=0 k=1 λ=0 k=1 j=0

i = 0, . . . , n − 1,

(2.2)

setting βij0 = 0 for j > n0 and βijk = 0 for j > n1 k = 1, . . . , l. Clearly, (2.2) is a linear system of n equations with respect to c0 , . . . , cn−1 . Let M :=

j! βij0 +

j l X X k=1

j! bj−λ βiλk k (j − λ)! λ=0

!n−1

i,j=0

be the matrix of the system (2.2). In the sequel we assume that the matrix M is regular. Observe that M is regular if and only if from all polynomials y of degree n − 1 only y = 0 satisfies the homogeneous boundary conditions n0 X j=0

βij0 y (j) (0) +

n1 l X X

βijk y (j) (bk ) = 0 ,

k=1 j=0

i = 0, . . . , n − 1 ,

(2.3)

corresponding to the conditions (1.2) with γi = 0, i = 0, . . . , n − 1. Indeed, substituting (2.1) with z = 0 into (2.3) we obtain a homogeneous system with respect to c0 , c1 , . . . , cn−1 . 4

This system coincides with the system (2.2) with γi = 0 (i = 0, . . . , n − 1) and z = 0. Therefore, the corresponding to (2.2) homogeneous system has only the trivial solution c0 = c1 = · · · = cn−1 = 0 (and thus M is regular) if and only if from all polynomials y of degree n − 1 only y = 0 satisfies (2.3). −1 Let M −1 = (pλ,µ )n−1 , the solution of the system (2.2) λ,µ=0 be the inverse of M . Using M can be written in the form c λ = δλ −

n1 l X X

κλjk (J α−j z)(bk ) ,

λ = 0, . . . , n − 1,

k=1 j=0

where δλ :=

n−1 X

pλµ γµ ,

κλjk :=

µ=0

n−1 X

pλµ βµjk .

(2.4)

µ=0

Therefore a function y in the form (2.1) satisfies the conditions (1.2) if and only if it can be expressed by the formula y = Gz + Q (2.5) where α

(Gz)(t) := (J z)(t) −

n−1 X

λ

t

κλjk (J α−j z)(bk ),

k=1 j=0

λ=0

Q(t) :=

n1 l X X

n−1 X

δλ t λ ,

λ=0

t ∈ [0, b] .

t ∈ [0, b] ,

(2.6)

(2.7)

Remark 2.1. In the case of an initial value problem (1.4) we have G = J α,

Q(t) =

n−1 X γλ λ=0

λ!

tλ .

Suppose now that y ∗ ∈ C[0, b] is a solution of the boundary value problem (1.1)–(1.2) such that z ∗ := D∗α y ∗ ∈ C[0, b]. Then it follows from the observations above that y ∗ has the form (see (2.5)) y ∗ = Gz ∗ + Q where G and Q are defined by the formulas (2.6) and (2.7), respectively. Inserting (2.5) into (1.1), we see that z ∗ = D∗α y ∗ satisfies the equation z(t) = f (t, (Gz)(t) + Q(t)),

t ∈ [0, b] .

(2.8)

Conversely, it is easy to show that if z ∗ ∈ C[0, b] is a solution of Eq. (2.8) then y ∗ = Gz ∗ +Q is a solution to (1.1)–(1.2). In this sense Eq. (2.8) is equivalent to the boundary value problem (1.1)–(1.2) and we can use it in constructing of high order methods for the numerical solution of problem (1.1)–(1.2). Observe that (2.8) is a nonlinear integral equation with respect to z. 5

An important question that arises by studying the attainable order of the convergence of a numerical method are the smoothness properties of the exact solution of a fractional differential equation. Some information about the smoothness properties of the solution can be obtained by using asymptotic expansions of the solution in fractional powers with respect to the independent variable (see, e.g., [1]). In the present paper we use another approach: in order to characterize the singularities of the derivatives of a solution of problem (1.1)–(1.2) at t = 0, we introduce a weighted space of smooth functions C q,ν (0, b] (cf., e.g., [2]). For given q ∈ N and ν ∈ R, ν < 1, by C q,ν (0, b] we denote the set of continuous functions y : [0, b] → R which are q times continuously differentiable in (0, b] and such that for all t ∈ (0, b] and i = 1, . . . , q the following estimate holds:   if i < 1 − ν   1 (i) y (t) ≤ c 1 + | log t| if i = 1 − ν .  1−ν−i  t if i > 1 − ν Here c = c(y) is a positive constant. Equipped with a suitable norm the set C q,ν (0, b] becomes a Banach space (see, e.g., [2, 35]). Clearly, C q [0, b] ⊂ C q,ν (0, b] ⊂ C m,µ (0, b] ⊂ C[0, b],

q ≥ m ≥ 1,

ν ≤ µ < 1.

(2.9)

Note that a function of the form y(t) = g1 (t) tµ +g2 (t) is included in C q,ν (0, b] if µ ≥ 1−ν > 0 and gj ∈ C q [0, b] , j = 1, 2. For the proof of Theorem 2.1 below we need following result from [15]. Lemma 2.1. Assume that the following conditions for Eq. (1.1) are fulfilled: ¯ → R is a continuous function which is q times (q ∈ N) (i) α > 0, α 6∈ N and f : Ω continuously differentiable in Ω where ¯ := {(t, y) : t ∈ [0, b], y ∈ R}; Ω

Ω := {(t, y) : t ∈ (0, b], y ∈ R},

(2.10)

(ii) there exist a monotonically increasing function ψ : [0, ∞) → R and a real number ν ∈ [1 − α, 1) such that for all nonnegative integers i and j with i + j ≤ q and for all (t, y) ∈ Ω   i+j 1 if i < 1 − ν   ∂ ≤ ψ(|y|) 1 + | log t| if i = 1 − ν f (t, y) ; (2.11) ∂ti ∂y j  1−ν−i  t if i > 1 − ν

if α ∈ (0, 1), then we assume in addition to (2.11) that for all nonnegative integers i and j with i + j ≤ q and for all (t, y1 ), (t, y2 ) ∈ Ω it holds i+j   ∂ 1 if i = 0 ∂ti ∂y j [f (t, y1 ) − f (t, y2 )] ≤ ψ(max{|y1 |, |y2 |})|y1 − y2 | t1−ν−i if i > 0 ; (2.12) (iii) for arbitrary given constants δi ∈ R (i = 0, . . . , n − 1) Eq. (1.1) possesses a solution y ∗∗ ∈ C[0, b] such that D∗α y ∗∗ ∈ C[0, b] and (y ∗∗ )(i) (0) = δi ,

i = 0, . . . , n − 1,

Then y ∗∗ ∈ C q,ν (0, b] and D∗α y ∗∗ ∈ C q,ν (0, b]. 6

n = dαe.

The regularity of a solution to (1.1)–(1.2) can be characterized by following theorem. Theorem 2.1. Let the conditions (i) and (ii) of Lemma 2.1 be fulfilled. Moreover, assume that (1.3) is true and from all polynomials y of degree dαe − 1 only y = 0 satisfies the conditions (2.3). Finally, suppose that the boundary value problem (1.1)–(1.2) possesses a solution y ∗ ∈ C[0, b] such that D∗α y ∗ ∈ C[0, b]. Then y ∗ ∈ C q,ν (0, b] and D∗α y ∗ ∈ C q,ν (0, b]. Proof. Let y ∗ ∈ C[0, b] be a solution to (1.1)–(1.2) such that D∗α y ∗ ∈ C[0, b]. Denote z ∗ := D∗α y ∗ and c∗λ

:= δλ −

l X n−1 X

κλjk (J α−j z ∗ )(bk ),

k=1 j=0

λ = 0, . . . , n − 1,

n = dαe,

where δλ and κλjk are defined by (2.4). Then (see (2.1)), ∗

α ∗

y (t) = (J z )(t) +

n−1 X λ=0

c∗λ tλ ,

t ∈ [0, b].

Consequently, y ∗ is a solution to Eq. (1.1) which satisfies the following initial conditions (see (1.8) and (1.6)): (y ∗ )(i) (0) = i! c∗i , i = 0, . . . , n − 1. This together with Lemma 2.1 yields the assertions of Theorem 2.1. 3. Piecewise polynomial interpolation For N ∈ N and 1 ≤ r < ∞ , let ΠN := {t0 , . . . , tN } be a partition (a graded grid) of the interval [0, b] with the grid points  r j tj = b , j = 0, 1, . . . , N , (3.1) N where r ∈ [1, ∞) is the so called grading exponent. If r = 1, then the grid points (3.1) are distributed uniformly; for r > 1 the grid points (3.1) are more densely clustered near the left endpoint of the interval [0,b]. (−1) For a given integer m ≥ 0 we denote by Sm (ΠN ) the standard space of piecewise polynomial functions on [0, b] associated with the grid ΠN :  (−1) Sm (ΠN ) = v : v (tj−1 ,tj ) ∈ πm , j = 1, . . . , N . Here v (tj−1 ,tj ) is the restriction of v : [0, b] → R onto the subinterval (tj−1 , tj ) and πm denotes (−1)

the set of polynomials of degree not exceeding m. Note that the elements of Sm (ΠN ) may have jump discontinuities at the interior points t1 , . . . , tN −1 of the grid ΠN . 7

In every interval [tj−1 , tj ], j = 1, . . . , N , we define m ∈ N collocation points tj1 , . . . , tjm by formula tjk = tj−1 + ηk (tj − tj−1 ) , k = 1, . . . , m, j = 1, . . . , N, (3.2)

where η1 . . . , ηm are some fixed (collocation) parameters which do not depend on j and N and satisfy 0 ≤ η1 < η2 < . . . < ηm ≤ 1 . (3.3) (−1)

For given N, m ∈ N let PN : C[0, b] → Sm−1 (ΠN ) be an interpolation operator such that (−1)

PN v ∈ Sm−1 (ΠN ), (PN v)(tjk ) = v(tjk ), k = 1, . . . , m, j = 1, . . . , N,

(3.4)

for any continuous function v ∈ C[0, b]. If η1 = 0, then by (PN v)(tj1 ) we denote the right limit limt→tj−1 ,t>tj−1 (PN v)(t). If ηm = 1, then by (PN v)(tjm ) we denote the left limit limt→tj ,t
N ∈ N,

with a constant c which is independent of N. Moreover, we have kPN z − zk∞ → 0

for every

z ∈ C[0, b]

as

N → ∞,

(3.5)

where v ∈ L∞ (0, b).

kvk∞ := sup |v(t)|, 0
(0, b], m ∈ N, −∞ < ν < 1, then  −r(1−ν) f or  N −m kPN z − zk∞ ≤ c N (1 + log N ) f or  −m N f or

If z ∈ C

m,ν

m 1 ≤ r < 1−ν m r = 1−ν =1 m r = 1−ν > 1 or r >

m 1−ν

  

.

(3.6)

Here r ∈ [1, ∞) is introduced in (3.1) and c is a constant not depending on N . 4. Smoothing transformation and numerical solutions Let us consider a change of variables t = b1−ρ τ ρ ,

τ ∈ [0, b],

(4.1)

depending on a parameter ρ ∈ [1, ∞). Then τ = b(ρ−1)/ρ t1/ρ ∈ [0, b] for t ∈ [0, b]. In the case ρ = 1 it follows from (4.1) that t = τ . We are interested in transformations (4.1) with ρ > 1 since such transformations possess a smoothing property for z ∈ C[0, b] ∩ C q (0, b] (q ∈ N) with singularities of derivatives of z(t) at t = 0. Actually, we obtain from [29, 30] the following result. 8

Lemma 4.1. Let z ∈ C q,ν (0, b] (q ∈ N, −∞ < ν < 1) and let zρ (τ ) := z(b1−ρ τ ρ ), τ ∈ [0, b], where ρ ∈ [1, ∞) if ν ∈ (0, 1) and ρ ∈ N if ν ≤ 0. Then zρ ∈ C q,νρ (0, b] where νρ := 1 − ρ(1 − ν). For the numerical solution of the boundary value problem (1.1)–(1.2) we use the following method. First we construct an equation of the form (see (2.8)) z = Tz

(4.2)

where (T z)(t) := f (t, (Gz)(t) + Q(t)), t ∈ [0, b], with G and Q given by (2.6) and (2.7), respectively. Introducing in (1.5) the change of variables t = b1−ρ τ ρ , s = b1−ρ σ ρ , τ, σ ∈ [0, b], ρ ∈ [1, ∞), we obtain that (J α z)(t) = (Jρα zρ )(τ ),

t = b1−ρ τ ρ ,

τ ∈ [0, b],

α > 0,

(4.3)

where zρ (τ ) := z(b1−ρ τ ρ ) and (Jρα zρ )(τ )

ρb(1−ρ)α := Γ(α)

Z

0

τ

(τ ρ − σ ρ )α−1 σ ρ−1 zρ (σ) dσ, τ ∈ [0, b],

α > 0.

(4.4)

Using in (4.2) the change of variables (4.1), we get for zρ (τ ) = z(b1−ρ τ ρ ) an integral equation of the form zρ = Tρ zρ (4.5) where (Tρ zρ )(τ ) := f (b1−ρ τ ρ , (Gρ zρ )(τ ) + Qρ (τ )) , n1 n−1 l X X X (1−ρ)λ ρλ α (Gρ zρ )(τ ) := (Jρ zρ )(τ ) − b τ κλjk (Jρα−j zρ )(bkρ ), Qρ (τ ) := Q(b

τ ),

(4.7)

k=1 j=0

λ=0

1−ρ ρ

(4.6)

bkρ :=

1/ρ b(ρ−1)/ρ bk

∈ (0, b],

τ ∈ [0, b],

with κλjk , Q and Jρα given by (2.4), (2.7) and (4.4), respectively. (−1) Approximations zρ,N ∈ Sm−1 (ΠN ) (m, N ∈ N) to the solution zρ∗ of Eq. (4.5) we find by collocation method from the conditions zρ,N (tjk ) = (Tρ zρ,N )(tjk ),

k = 1, . . . , m, j = 1, . . . , N,

(4.8)

with {tjk }, defined by (3.2). If η1 = 0, then by zρ,N (tj1 ) we denote the right limit limt→tj−1 ,t>tj−1 zρ,N (t). If ηm = 1, then zρ,N (tjm ) denotes the left limit limt→tj ,t
(4.9)

where Gρ is defined by (4.7) and Qρ (τ ) = Q(b1−ρ τ ρ ), τ ∈ [0, b]. Further, approximations yρ,N to yρ∗ we find by the formula yρ,N = Gρ zρ,N + Qρ

(4.10)

(−1)

where zρ,N ∈ Sm−1 (ΠN ) is determined by (4.8). Finally, approximations yN (t) to the solution y ∗ (t) of problem (1.1)–(1.2) we find by setting yN (t) := yρ,N (b(ρ−1)/ρ t1/ρ ),

t ∈ [0, b]. (−1)

Note that the conditions (4.8) for finding zρ,N ∈ Sm−1 (ΠN ) have an operator equation representation zρ,N = PN Tρ zρ,N (4.11) where PN and Tρ are defined by (3.4) and (4.6), respectively. The collocation conditions (4.8) form a system of equations whose exact form is determined (−1) by the choice of a basis in the space Sm−1 (ΠN ). If η1 > 0 or ηm < 1, then we can use the Lagrange fundamental polynomial representation: zρ,N (τ ) =

m N X X

cβγ Lβγ (τ ) ,

β=1 γ=1

τ ∈ [0, b] ,

(4.12)

where Lβγ (τ ) = 0 for τ 6∈ [tβ−1 , tβ ] and m Y

Lβγ (τ ) =

i=1,i6=γ

τ − tβi tβγ − tβi

for τ ∈ [tβ−1 , tβ ], γ = 1, . . . , m, β = 1, . . . , N.

(−1)

Then zρ,N ∈ Sm−1 (ΠN ) and zρ,N (tjk ) = cjk , k = 1, . . . , m, j = 1, . . . , N . Searching the solution to (4.8) in the form (4.12), we obtain a system of nonlinear algebraic equations with respect to the coefficients {cjk }: cjk = f



b1−ρ tρjk ,

N X m X

(Gρ Lβγ )(tjk )cβγ

β=1 γ=1

 + Qρ (tjk ) ,

k = 1, . . . , m, j = 1, . . . , N. (4.13)

Remark 4.1. In the case of initial value problem (1.4) Gρ = Jρα (see Remark 2.1). Since (Jρα Lβγ )(tjk ) = 0 if β > j, in this case the coefficients cj1 , . . . , cjm can be found for every j = 1, . . . , N from the system cjk = f



b1−ρ tρjk ,

j m X X

(Jρα Lβγ )(tjk )cβγ

β=1 γ=1

with m unknowns cjk , k = 1, . . . , m. 10

 + Qρ (tjk ) ,

k = 1, . . . , m,

(4.14)

Having {cjk } in hand, we get with the help of (4.10) and (4.12) that yρ,N (τ ) =

N X m X

cβγ (Gρ Lβγ )(τ ) + Qρ (τ ),

β=1 γ=1

τ ∈ [0, b].

(4.15)

Note that this algorithm can also be used in the case if in (3.3) η1 = 0 and ηm = 1. In this case we have tjm = tj+1,1 = tj and cjm = cj+1,1 (j = 1, . . . , N − 1), and hence in the system (4.12) there are (m − 1)N + 1 equations and unknowns. 5. Convergence analysis First of all we introduce some notations and definitions. Let X be a Banach space with a norm kxk, x ∈ X. A sequence {AN }N ∈N of operators AN ∈ L(X, X) is called compactly converging to A ∈ L(X; X) (we write AN → A compactly) if AN x → Ax as N → ∞ for every x ∈ X and for any bounded sequence {xN }N ∈N , xN ∈ X, it follows that the sequence {AN xN }N ∈N is relatively compact in X (i.e. every subsequence {AN xN }N ∈N0 ⊂N contains a subsequence {AN xN }N ∈N00 ⊂N0 converging in X). Let us consider the nonlinear equations x = Sx

(5.1)

and x = SN x,

N ∈ N,

(5.2)

where S : B → X and SN : B → X are nonlinear operators defined on an open set B ⊂ X. We recall that S : B → X is called Fr´echet differentiable at x0 ∈ B if there exists a linear operator S 0 (x0 ) ∈ L(X, X) such that kSx − Sx0 − S 0 (x0 )(x − x0 )k/kx − x0 k → 0 as kx − x0 k → 0; in this case S 0 (x0 ) is the (unique) Fr´echet derivative of S at x0 . To study the convergence behavior of our method we use the following result from the approximation theory by Vainikko (see Theorem 2 in [36]). Lemma 5.1. Let the following conditions be fulfilled: 10 equation (5.1) has a solution x∗ ∈ B, and the operator S is Fr´ echet differentiable at x∗ ; 20 there is a positive number δ such that the operators SN (N ∈ N) are Fr´ echet differentiable ∗ in the ball kx − x k ≤ δ which is assumed to be contained in B, and for any ε > 0 there is a δε ∈ (0, δ] such that for every N ∈ N 0 0 kSN (x) − SN (x∗ )kL(X,X) ≤ ε whenever

kx − x∗ k ≤ δε ;

30 kSN x∗ − Sx∗ k → 0 as N → ∞; 0 0 40 SN (x∗ ) → S 0 (x∗ ) compactly, whereby SN (x∗ ) ∈ L(X, X) (N ∈ N) are compact and the 0 ∗ homogeneous equation x = S (x )x has in X only the trivial solution x = 0. 11

Then there exist N0 ∈ N and δ0 ∈ (0, δ] such that equation (5.2) has for N ≥ N0 a unique solution xN in the ball kx − x∗ k ≤ δ0 . Thereby xN → x∗ as N → ∞ and the following error estimate holds: kxN − x∗ k ≤ c kSN x∗ − Sx∗ k, N ≥ N0 . (5.3) Here c is a positive constant not depending on N . The following lemma presents some properties of Jρα which follow from (4.3) and from the corresponding properties of the Riemann-Liouville integral operator J α defined by (1.5). Lemma 5.2. Let α > 0 and ρ ≥ 1 be some given real numbers. Then Jρα defined by (4.4) is linear, bounded and compact as an operator from L∞ (0, b) into C[0, b]. Moreover, we have for any x ∈ L∞ (0, b) that Jρα Jρβ x = Jρα+β x,

α > 0,

β > 0,

ρ ≥ 1.

Now we are ready to prove the convergence of method (4.10)–(4.11) and study its attainable order of convergence for arbitrary collocation parameters η1 , . . . , ηm satisfying (3.3). Theorem 5.1. Suppose (1.3) and let problem (1.1)–(1.2) have a solution y ∗ ∈ C[0, b] such ¯ → R be a continuous function such that ∂ f (t, y) is that z ∗ := D∗α y ∗ ∈ C[0, b]. Let f : Ω ∂y ∂2 ¯ continuous for (t, y) ∈ Ω, 2 f (t, y) is continuous for (t, y) ∈ Ω and ∂y

j ∂ ≤ ψ(|y|), f (t, y) ∂y j

(t, y) ∈ Ω,

j = 0, 1, 2.

¯ are defined by (2.10) and ψ : [0, ∞) → R is a monotonically increasing Here Ω and Ω function. Moreover, assume that from all polynomials y of degree n − 1 (n = dαe) only y = 0 satisfies the homogeneous boundary conditions (2.3) and from all solutions y ∈ C[0, b] of the linear homogeneous (fractional) differential equation (D∗α y)(t) =


∂ f t, y ∗ (t) y(t), ∂y

t ∈ [0, b],

(5.4)

only y = 0 satisfies the conditions (2.3). Finally, let m ∈ N and assume that the collocation points (3.2) with grid points (3.1) and arbitrary parameters η1 , . . . , ηm satisfying (3.3) are used. Then there exist N0 ∈ N and δ0 > 0 such that, for all N ≥ N0 , Eq. (4.11) possesses a (−1) unique solution zρ,N ∈ Sm−1 (ΠN ) in the ball kx − zρ∗ k∞ ≤ δ0 where zρ∗ (τ ) = z ∗ (b1−ρ τ ), τ ∈ [0, b], z ∗ = D∗α y ∗ and ρ ≥ 1. Moreover, kyN − y ∗ k∞ → 0

as N → ∞ ,

(5.5)

where yN (t) = yρ,N (b(ρ−1)/ρ t1/ρ ), t ∈ [0, b], and yρ,N is defined with the help of zρ,N by (4.10). 12

If, in addition, the assumptions of Theorem 2.1 are fulfilled with q = m and some ν ∈ [1 − α, 1), then for all N ≥ N0 and r ≥ 1 the following error estimate holds:  −ρr(1−ν)  m f or 1 ≤ ρr < 1−ν  N  m =1 kyN − y ∗ k∞ ≤ c N −m (1 + log N ) f or ρr = 1−ν . (5.6)  −m m m  N f or ρr = 1−ν > 1 or ρr > 1−ν

Here c is a constant not depending on N , ρ ∈ [1, ∞) if ν ∈ (0, 1) and ρ ∈ N if ν ≤ 0.

Proof. To apply Lemma 5.1 we consider (4.5) and (4.11) as operator equations (5.1) and (5.2) in the space X := L∞ (0, b) where S := Tρ and SN := PN Tρ with Tρ and PN , defined by (4.6) and (3.4), respectively. 0 First we find the Fr´echet derivative Tρ (x0 ) for Tρ at x0 ∈ L∞ (0, b). If x and x0 belong to L∞ (0, b), then, due to Taylor formula,

∂ f t, (Gρ x0 )(τ ) + Qρ (τ ) Gρ (x − x0 ) (τ ) ∂y 2

2 1 ∂ + f t, (Gρ x0 )(τ ) + ξ(τ )(Gρ (x − x0 ))(τ ) + Qρ (τ ) Gρ (x − x0 ) (τ ) , 2 2 ∂y (Tρ x)(τ ) − (Tρ x0 )(τ ) =

where ξ(τ ) ∈ [0, 1] and t = b1−ρ τ ∈ (0, b]. From this it follows that # "
∂ f (b1−ρ τ, y) Tρ0 (x0 )x (τ ) = (Gρ x)(τ ) , ∂y y=(Gρ x0 )(τ )+Qρ (τ )

(5.7)

(5.8)

where τ ∈ [0, b] and x0 , x ∈ L∞ (0, b). Observe that Tρ0 (x0 ) ∈ L(L∞ (0, b), C[0, b]) if x0 ∈ L∞ (0, b). From (5.7) we see that
0 SN (x0 )x = (PN Tρ )0 (x0 )z = PN Tρ0 (x0 )x , x0 , x ∈ L∞ (0, b),

with Tρ0 (x0 )x, defined by (5.8). It follows from Lemma 5.2 that Gρ (defined by (4.7)) is linear, bounded and compact as an operator from L∞ (0, b) into C[0, b]. Since from all solutions y of Eq. (5.4) only y = 0 satisfies (2.3), Eq. x = Tρ0 (zρ∗ )x has in L∞ (0, b) only the trivial solution x = 0. Using these observations and Lemma 5.2 we can check that the operators S = Tρ and SN = PN Tρ satisfy the conditions 10 − 40 of Lemma 5.1. Hence, in agreement with Lemma 5.1, there exist N0 ∈ N and δ0 > 0 such that, for all N ≥ N0 , Eq. (4.11) possesses a unique solution zρ,N in the ball kx − zρ∗ k∞ ≤ δ0 , and kzρ,N − zρ∗ k∞ ≤ c kPN zρ∗ − zρ∗ k∞ , with a positive constant c which is independent of N . Since yρ,N = Gρ zρ,N + Qρ and yρ∗ = Gρ zρ∗ + Qρ , we obtain that kyN − y ∗ k∞ ≤ c kzρ,N − zρ∗ k∞ ≤ c1 kPN zρ∗ − zρ∗ k∞ , 13

(5.9)

with some positive constants c and c1 which are independent of N . This together with zρ∗ ∈ C[0, b] and (3.5) yields the convergence (5.5). If the assumptions of Theorem 2.1 are fulfilled with q = m and ν ∈ [1 − α, 1), then z ∗ ∈ C m,ν (0, b], and we get by Lemma 4.1 that zρ∗ ∈ C m,νρ (0, b], νρ = 1 − ρ(1 − ν). This together with (5.9) and (3.6) yields the estimate (5.6). 6. Superconvergence results It follows from Theorem 5.1 that in the case of sufficiently smooth f (t, y), using sufficiently large values of ρr, for method (4.10)–(4.11) by every choice of collocation parameters 0 ≤ η1 < · · · < ηm ≤ 1 a convergence of order O(N −m ) can be expected. From Theorem 6.1 below we see that by a careful choice of parameters η1 , . . . , ηm it is possible to establish a faster convergence of our method. In order to study the superconvergence properties of our method by a special choice of collocation parameters we need the following lemma from [35]. Lemma 6.1. Let z ∈ C m+1,ν (0, b] with some m ∈ N and ν ∈ (−∞, 1). Let N ∈ N, α ∈ (0, 1], r ∈ [1, ∞), ρ ∈ [1, ∞) if ν ∈ (0, 1) and ρ ∈ N if ν ≤ 0, zρ (τ ) = z(b1−ρ τ ρ ), τ ∈ [0, b]. Let {tj }, PN and Jρα be defined by (3.1), (3.4) and (4.4), respectively. Moreover, assume that a quadrature approximation Z 1 m X wk g(ηk ) (6.1) g(x) dx ≈ 0

k=1

with the knots {ηk } satisfying (3.3) and appropriate weights {wk } is exact for all polynomials g of degree m. Then the following estimate holds:   EN (m, α, ν, ρr) if 0<α<1 α kJρ (PN zρ − zρ )k∞ ≤ c . (6.2) ∗ EN (m, ν, ρr) if α=1 Here c is a constant not depending on N , and  −ρr(α+1−ν) N if    N −m−α (1 + log N ) if EN (m, α, ν, ρr) := N −m−α if    or  −ρr(2−ν) N    N −m−1 (1 + log N )2 ∗ (m, ν, ρr) := EN N −m−1 (1 + log N )    −m−1 N

if if if if

m+α 1 ≤ ρr < α+1−ν m+α ρr = α+1−ν =1 m+α ρr = α+1−ν > 1 m+α ρr > α+1−ν

1 ≤ ρr < m+1 2−ν ρr = m+1 =1 2−ν ρr = m+1 >1 2−ν m+1 ρr > 2−ν

      

      

.

,

(6.3)

(6.4)

Theorem 6.1. Assume that the following conditions are fulfilled: (i) problem (1.1)–(1.2) has a solution y ∗ ∈ C[0, b] and (1.3) holds; (ii) m ∈ N and problem (1.1)–(1.2) satisfies the conditions of Theorem 2.1 with q = m + 1 and ν ∈ [1 − α, 1); 14

(iii) from all solutions y ∈ C[0, b] to Eq. (5.4) only y = 0 satisfies the boundary conditions (2.3); (iv) ρ is a smoothing parameter (see (4.1)) such that ρ ∈ [1, ∞) if ν ∈ (0, 1) and ρ ∈ N if ν ≤ 0; (v) PN is defined by (3.4) where the collocation points (3.2) are generated by grid points (3.1) and by parameters 0 ≤ η1 < η2 < · · · < ηm ≤ 1 so that a quadrature approximation (6.1) is exact for all polynomials of degree m. Then there exist N0 ∈ N and δ0 > 0 such that for all N ≥ N0 Eq. (4.11) possesses a unique solution zρ,N in the ball kx−zρ∗ k ≤ δ0 where zρ∗ (τ ) = z ∗ (b1−ρ τ ), τ ∈ [0, b], z ∗ = D∗α y ∗ . Moreover, for N ≥ N0 the following error estimates hold:  −ρr(α−n1 +1−ν) m+α−n1  N if 1 ≤ ρr < α−n      N −m−α+n1 (1 + log N ) if ρr = m+α−n11 +1−ν  = 1 ∗ α−n1 +1−ν kyN − y k∞ ≤ c (6.5) m+α−n1 N −m−α+n1 if ρr = α−n >1   +1−ν   1   m+α−n1 or ρr > α−n 1 +1−ν for α − n1 < 1 and

kyN − y ∗ k∞

 −ρr(2−ν) N    N −m−1 (1 + log N )2 ≤c N −m−1 (1 + log N )    −m−1 N

if if if if

1 ≤ ρr < m+1 2−ν ρr = m+1 = 1 2−ν m+1 ρr = 2−ν > 1 ρr > m+1 2−ν

      

(6.6)

for α − n1 ≥ 1. In (6.5) and (6.6) yN (t) = yρ,N (b(ρ−1)/ρ t1/ρ ), t ∈ [0, b], with yρ,N given by (4.10), r ∈ [1, ∞) is given by (3.1) and c is a positive constant not depending on N . ˆ0 ∈ N and δˆ0 > 0 Proof. Since q = m + 1 ≥ 2, it follows from Theorem 5.1 that there exist N ˆ such that for N ≥ N0 Eq. (4.11) has a unique solution zρ,N in the ball kx − zρ∗ k∞ ≤ δˆ0 . Denote ˆ0 , zˆρ,N := Tρ zρ,N , N ≥ N (6.7) with Tρ defined by the formula (4.6). Taking into account (4.11), we obtain that PN zˆρ,N = zρ,N . Substituting zρ,N = Pρ,N zˆρ,N into (6.7) we see that zˆρ,N is a solution of the equation zˆρ,N = Tρ PN zˆρ,N ,

ˆ0 . N ≥N

(6.8)

Consider the equations (4.5) and (6.8) as operator equations (5.1) and (5.2) in the space X := C[0, b] where S := Tρ and SN := Tρ PN . In a similar way as we obtained the formula (5.8), we get for the Fr´echet derivative of SN = Tρ PN at x0 ∈ C[0, b] the following formula: # "
∂ 0 (x0 )x (τ ) = f (b1−ρ τ, y) (Gρ PN x)(τ ), τ ∈ [0, b], x ∈ C[0, b]. SN ∂y y=(Gρ PN x0 )(τ )+Qρ (τ )

Moreover, we can check that the operators S = Tρ and SN = Tρ PN satisfy the conditions ˆ0 and δ0 > 0 10 − 40 of Lemma 5.1. From this Lemma it follows that there exist N0 ≥ N 15

such that, for N ≥ N0 , Eq. (6.8) possesses a unique solution zˆρ,N in the ball kx − zρ∗ k∞ ≤ δ0 whereby (6.9) kˆ zρ,N − zρ∗ k∞ ≤ c kTρ PN zρ∗ − Tρ zρ∗ k∞ , N ≥ N0 . Here and below by c and also c1 we denote positive constants which are independent of N . As
(Tρ PN zρ∗ )(τ ) − (Tρ zρ∗ )(τ ) = f b1−ρ τ, (Gρ PN zρ∗ )(τ ) + Qρ (τ )
− f b1−ρ τ, (Gρ zρ∗ )(τ ) + Qρ (τ ) , τ ∈ [0, b],

then with the help of Lagrange formula we obtain that

kTρ PN zρ∗ − Tρ zρ∗ k∞ ≤ c kGρ (PN zρ∗ − zρ∗ )k∞

(6.10)

Using (4.7) and Lemma 5.2 we get kGρ (PN zρ∗ − zρ∗ )k∞ ≤ kJρα (PN zρ∗ − zρ∗ )k∞

n1 l X X α−j
Jρ (PN zρ∗ − zρ∗ ) (bkρ ) ≤ c1 kJρα−n1 (PN zρ∗ − zρ∗ )k∞ +c

(6.11)

k=1 j=0

1/ρ

where bkρ = b(ρ−1)/ρ bk

∈ [0, b]. From (6.9)–(6.11) it follows that

kˆ zρ,N − zρ∗ k∞ ≤ c kJρα−n1 (PN zρ∗ − zρ∗ )k∞ ,

N ≥ N0 .

(6.12)

Further, as yρ,N = Gρ zρ,N + Qρ , yρ∗ = Gρ zρ∗ + Qρ and zρ,N = PN zˆρ,N , then kyN − y ∗ k∞ = kyρ,N − yρ∗ k∞ = kGρ (zρ,N − zρ∗ )k∞ ≤ kGρ PN (ˆ zρ,N − zρ∗ )k∞ + kGρ (PN zρ∗ − zρ∗ )k∞ ,

N ≥ N0 .

This together with (6.11) and (6.12) yields kyN − y ∗ k∞ ≤ c kJρα−n1 (PN zρ∗ − zρ∗ )k∞ ,

N ≥ N0 .

(6.13)

Because of Theorem 2.1 we have z ∗ ∈ C m+1,ν (0, b] and, due to Lemma 5.2, kJρα−n1 (PN zρ∗ − zρ∗ )k∞ ≤ c kJρ1 (PN zρ∗ − zρ∗ )k∞ ,

α − n1 ≥ 1,

N ≥ N0 .

This together with (6.13) and Lemma 6.1 yields the estimates (6.5) and (6.6). 7. Numerical experiments In this section, we present some numerical experiments to demonstrate the accuracy of the modified spline collocation method (4.10)–(4.11) and compare the actual convergence rate with the theoretical error estimates of Theorem 6.1. For obtaining numerical results Fortran-98 language has been used. 16

Example 1. Consider the following initial value problem:
D∗1/3 y (t) = f (t, y(t)), y(0) = 1, 0 ≤ t ≤ 1,

where

f (t, y) :=

(7.1)

1 Γ(5/3) 1/3 1 3 ty − t(t2/3 + 1)3 + t . 10 10 Γ(4/3)

The function f (t, y) satisfies the conditions (2.11) and (2.12) with ν = 2/3 and arbitrary ∗ 2/3 q ∈ N. An exact solution of this  problem is y (t) = t + 1 (see Figure 1(a)). Its Caputo  1/3 derivative is z ∗ (t) := D∗ y ∗ (t) = Γ(5/3) t1/3 . Γ(4/3) 1/3

To solve (7.1) with method (4.10)–(4.11) we set z := D∗ y. Then a solution of (7.1) can be presented in the form y ∗ (t) = (J 1/3 z ∗ )(t) + 1 where z ∗ is a solution of the integral equation
z(t) = f t, (J 1/3 z)(t) + 1 , t ∈ [0, 1]. Using a change of variables t = τ ρ , τ ∈ [0, 1], ρ ≥ 1, and notations yρ (τ ) := y(τ ρ ), zρ (τ ) := z(τ ρ ),

τ ∈ [0, 1], ρ ≥ 1,

(7.2)

we get that yρ∗ (τ ) = (Jρ1/3 zρ∗ )(τ ) + 1,

τ ∈ [0, 1],

where zρ∗ = z ∗ (τ ρ ) is a solution of the integral equation
zρ (τ ) = f τ ρ , (Jρ1/3 zρ )(τ ) + 1 ,

Here

(Jρ1/3 zρ )(τ )

ρ = Γ(1/3)

Z

0

τ

τ ∈ [0, 1].

(τ ρ − σ ρ )−2/3 σ ρ−1 zρ (σ) dσ,

(7.3)

τ ∈ [0, 1].

Since τ = t1/ρ , we have y ∗ (t) = yρ∗ (t1/ρ ), t ∈ [0, 1]. For the numerical solution of Eq. (7.3) we use spline collocation method (4.8). We determine approximations zρ,N to the solution zρ∗ of (7.3) in the form (4.12), finding the coefficients {cjk } from collocation conditions (4.14): cjk

  j m X X ρ 1/3 (Jρ Lβγ )(tjk )cβγ + 1 , = f tjk ,

k = 1, . . . , m, j = 1, . . . , N.

(7.4)

β=1 γ=1

For calculations we use grid points tj = (j/N )r , j = 0, . . . , N, r ≥ 1, and collocation points tjk =√tj−1 + ηk (tj − tj−1 ), k = 1, . . . , m, j = 1, . . . , N, where in the case m = 2, η1 = √ (3 − 3)/6, η2 = 1 − η1 , and in the case m = 3, η1 = (5 − 15)/10, η2 = 1/2, η3 = 1 − η1 . These parameters {ηk } are the the node points of the Gaussian quadrature rule (6.1) which in the case m = 2 is exact for all polynomials of degree 3 and in the case m = 3 is exact for all polynomials of degree 5. 17

Weakly singular integrals in (7.4) are found by the package QUADPACK and the coefficients {cjk } in (4.12) are determined from the nonlinear system (7.4) by Newton method. Using these coefficients we find approximations yρ,N to the solution y ∗ of problem (7.1) by the formula (4.15): yρ,N (τ ) =

N X m X

cβγ (Jρ1/3 Lβγ )(τ ) + 1,

β=1 γ=1

τ ∈ [0, 1].

In Tables 7.1–7.3 some results of numerical experiments for different values of the parameters N , ρ and r are presented. The results in Tables 7.1 and 7.2 are calculated using uniform grids (r = 1) and the results in Table 7.3 using nonuniform grids (r > 1). The quantities εN in these tables are the approximate values of the norm kyN − y ∗ k∞ , calculated as follows: εN := max

max |yN (sjk ) − y ∗ (sjk )|.

(7.5)

j=1,...,N k=0,...,10

Here sjk := tj−1 + k(tj − tj−1 )/10, k = 0, . . . , 10, j = 1, . . . , N , and yN (t) = yρ,N (t1/ρ ), 0 ≤ t ≤ 1. The ratios %N := εN/2 /εN , characterizing the observed convergence rate, are also presented. Moreover, in the Figure 1(b) the graph of the error yN (t) − y ∗ (t) for the problem (7.1) in the case m = 3, r = 2, ρ = 2, N = 64 is displayed. We see that the approximate solution yN (t) oscillates around the exact solution y ∗ (t) for t ∈ [0, 1]. Table 7.1: Numerical results for problem (7.1) in the case m = 2.

N 4 8 16 32 64 128 256

r = 1, ρ = 1 εN %N −2 2.27 · 10 1.42 · 10−2 1.592 8.95 · 10−3 1.589 5.64 · 10−3 1.588 3.55 · 10−3 1.588 2.24 · 10−3 1.587 1.41 · 10−3 1.587 1.587

r = 1, ρ = 2 εN %N −3 2.80 · 10 1.12 · 10−3 2.513 4.43 · 10−4 2.519 1.76 · 10−4 2.520 6.98 · 10−5 2.520 2.77 · 10−5 2.520 1.10 · 10−5 2.520 2.520

r = 1, ρ = 3.5 εN %N −4 6.10 · 10 1.21 · 10−4 5.039 2.40 · 10−5 5.040 4.77 · 10−6 5.040 9.46 · 10−7 5.040 1.88 · 10−7 5.040 3.73 · 10−8 5.040 5.040

r = 1, ρ = 4 εN %N −4 9.17 · 10 1.93 · 10−4 4.754 3.98 · 10−5 4.844 7.95 · 10−6 5.011 1.57 · 10−6 5.072 3.08 · 10−7 5.093 6.03 · 10−8 5.103 5.040

Since for problem (7.1) α = 1/3, n1 = 0 and ν = 2/3, it follows from (6.5) that, for ρ ≥ 1, r ≥ 1 and sufficiently large N , in case m = 2,  −2ρr/3  N if 1 ≤ ρr ≤ 3.5 εN ≤ c , N −7/3 if ρr ≥ 3.5 18

Table 7.2: Numerical results for problem (7.1) in the case m = 3.

N 4 8 16 32 64 128 256

r = 1, ρ = 1 εN %N 7.27 · 10−3 4.57 · 10−3 1.591 2.88 · 10−3 1.589 1.81 · 10−3 1.588 1.14 · 10−3 1.588 7.19 · 10−4 1.587 4.53 · 10−4 1.587 1.587

r = 1, ρ = 2 εN %N 7.09 · 10−4 2.81 · 10−4 2.521 1.12 · 10−4 2.520 4.43 · 10−5 2.520 1.76 · 10−5 2.520 6.98 · 10−6 2.520 2.77 · 10−6 2.520 2.520

r = 1, ρ = 4 εN %N 8.43 · 10−5 1.33 · 10−5 6.349 2.09 · 10−6 6.350 3.29 · 10−7 6.350 5.19 · 10−8 6.350 8.17 · 10−9 6.350 1.29 · 10−9 6.350 6.350

r = 1, ρ = 5 εN %N 3.68 · 10−5 3.65 · 10−6 10.08 3.62 · 10−7 10.08 3.60 · 10−8 10.08 3.57 · 10−9 10.08 4.50 · 10−10 7.93 −10 7.90 · 10 0.57 10.08

Table 7.3: Numerical results for problem (7.1) in the case m = 3.

N 4 8 16 32 64 128 256

r = 2, ρ = 1 εN %N 2.88 · 10−3 1.14 · 10−3 2.521 4.53 · 10−4 2.520 1.80 · 10−4 2.520 7.13 · 10−5 2.520 2.83 · 10−5 2.520 1.12 · 10−5 2.520 2.520

r = 2, ρ = 2 r = 2, ρ = 2.5 r = 3, ρ = 2 εN %N εN %N εN %N 1.48 · 10−4 5.42 · 10−5 2.62 · 10−4 2.34 · 10−5 6.329 5.69 · 10−6 9.53 2.86 · 10−5 9.16 −6 −7 −6 3.69 · 10 6.349 5.74 · 10 9.91 3.06 · 10 9.34 5.81 · 10−7 6.350 5.73 · 10−8 10.02 3.16 · 10−7 9.68 −8 −9 −8 9.16 · 10 6.350 5.70 · 10 10.05 3.21 · 10 9.84 1.44 · 10−8 6.350 5.92 · 10−10 9.62 3.26 · 10−9 9.83 −9 −10 −10 2.27 · 10 6.350 6.07 · 10 0.97 5.47 · 10 5.97 6.350 10.08 10.08

(a)

(b) 2

Figure 1: (a) The exact solution y ∗ (t) = t 3 + 1 for (7.1). (b) The error yN (t) − y ∗ (t) for (7.1) in the case m = 3, r = 2, ρ = 2, N = 64.

19

and in case m = 3, εN ≤ c



N −2ρr/3 if N −10/3 if

1 ≤ ρr ≤ 5 ρr ≥ 5



.

Due to these error estimates the ratios %N in the case m = 2 for ρr = 1, ρr = 2 and ρr ≥ 3.5 ought to be approximately 22/3 ≈ 1.587, 24/3 ≈ 2.520 and 27/3 ≈ 5.040, respectively. In the case m = 3 for ρr = 1, ρr = 2, ρ = 4 and ρr ≥ 5 the ratios %N ought to be approximatively 22/3 ≈ 1.587, 24/3 ≈ 2.520, 28/3 ≈ 6.350 and 210/3 ≈ 10.079, respectively. These values are given in the last row of Tables 7.1–7.3. We see from Tables 7.1–7.3 that the theoretical estimate (6.5) of Theorem 6.1 is in good accordance with the actual convergence rate of the proposed method. The irregular behavior of the ratios %N in the Tables 7.2–7.3 for εN < 10−9 can be explained by round-off errors in the computations. Example 2. Consider the following boundary value problem:
D∗1.5 y (t) = f (t, y(t)), y(0) = −1, y 0 (1) = 1.9, 0 ≤ t ≤ 1, (7.6) where

f (t, y) :=


2 Γ(2.9) 0.4 1 2 1 1.9 y − t −1 + t . 2 2 Γ(1.4)

The exact solution of this problem is y ∗ (t) = t1.9 − 1 (see Figure 2(a)). Its derivative z ∗ (t) := (D∗1.5 y ∗ ) (t) = Γ(2.9) t0.4 . Γ(1.4) Denote z := (D∗1.5 y). Then a solution y ∗ of (7.6) can be presented in the form y ∗ = Gz ∗ +Q where

(Gz)(t) := J 1,5 z (t) − t J 0.5 z (1), Q(t) := −1 + 1.9 t, t ∈ [0, 1], and z ∗ is a solution of the integral equation

z(t) = f (t, (Gz) (t) + Q(t)) , yρ∗

t ∈ [0, 1].

Using the
change of variables t = τ ρ , τ ∈ [0, 1], ρ ≥ 1, and notations (7.2) we get that = Gρ zρ∗ + Qρ where

(Gρ zρ )(τ ) := Jρ1,5 zρ (τ ) − τ ρ Jρ0.5 zρ (1),

Qρ (τ ) := −1 + 1.9 τ ρ ,

τ ∈ [0, 1],

and zρ∗ (τ ) = z ∗ (τ ρ ) is a solution of the integral equation

zρ (τ ) = f τ ρ , Gρ zρ∗ (τ ) + Qρ (τ ) , τ ∈ [0, 1].

(7.7)

Approximations zρ,N to the solution zρ∗ of (7.7) we find in the form (4.12) where the coefficients {cjk } are found from the nonlinear system (4.13) by Newton method. At that the same grid points {tj } and collocation points {tjk } as in Example 1 are used. After that the approximations yρ,N to the solution y ∗ of problem (7.6) by the formula (4.15) are found. In Tables 7.4 and 7.5 the errors εN calculated by (7.5) and their ratios %N := εN/2 /εN , are presented. In Figure 2(b) the graph of the error yN (t) − y ∗ (t) for (7.6) in the case m = 3, r = 1, ρ = 4, N = 64 is displayed. 20

Since for problem (7.6) α = 1.5, n1 ρ ≥ 1, r ≥ 1 and sufficiently large N , in  −0.9ρr N εN ≤ c N −2.5

= 1 and ν = 0.6, it follows from (6.5) that, for case m = 2,  ≈ 2.78 if 1 ≤ ρr ≤ 25 9 , if ρr ≥ 25 9

and in case m = 3, εN ≤ c



N −0.9ρr if N −3.5 if

1 ≤ ρr ≤ ρr ≥ 35 9

35 9

≈ 3.89



.

Table 7.4: Numerical results for problem (7.6) in the case m = 2.

N 4 8 16 32 64 128 256

r = 1, ρ = 1 εN %N −3 1.44 · 10 5.39 · 10−4 2.68 2.10 · 10−4 2.56 8.07 · 10−5 2.60 3.08 · 10−5 2.62 1.17 · 10−5 2.63 4.43 · 10−6 2.64 1.87

r = 1, ρ = 2 εN %N −4 3.84 · 10 6.07 · 10−5 6.33 1.01 · 10−5 5.99 1.73 · 10−6 5.84 3.00 · 10−7 5.77 5.23 · 10−8 5.74 9.15 · 10−9 5.72 3.48

r = 1, ρ = 3 r = 1, ρ = 4 εN %N εN %N −4 −3 6.47 · 10 3.05 · 10 1.03 · 10−4 6.28 4.79 · 10−4 6.38 1.74 · 10−5 5.91 8.06 · 10−5 5.94 3.02 · 10−6 5.77 1.39 · 10−5 5.78 5.28 · 10−7 5.72 2.44 · 10−6 5.72 9.27 · 10−8 5.69 4.28 · 10−7 5.70 1.63 · 10−8 5.69 7.52 · 10−8 5.69 5.66 5.66

Table 7.5: Numerical results for problem (7.6) in the case m = 3.

N 4 8 16 32 64 128 256

r = 1, ρ = 1 εN %N 5.78 · 10−4 2.14 · 10−4 2.70 8.06 · 10−5 2.66 3.05 · 10−5 2.65 1.15 · 10−5 2.64 4.36 · 10−6 2.64 1.65 · 10−6 2.64 1.87

r = 1, ρ = 2 εN %N 4.64 · 10−6 2.92 · 10−7 15.90 3.74 · 10−8 7.80 6.40 · 10−9 5.85 1.04 · 10−9 6.18 1.60 · 10−10 6.49 2.09 · 10−11 7.64 3.48

r = 1, ρ = 3 εN %N 9.16 · 10−6 7.29 · 10−7 12.56 6.16 · 10−8 11.84 5.32 · 10−9 11.57 4.57 · 10−10 11.66 3.45 · 10−11 13.25 3.41 · 10−12 10.09 6.50

r = 1, ρ = 4 εN %N 2.02 · 10−5 1.63 · 10−6 12.35 1.40 · 10−7 11.69 1.22 · 10−8 11.48 1.06 · 10−9 11.43 8.61 · 10−11 12.36 2.34 · 10−12 36.74 11.31

Due to these error estimates the ratios %N in the case m = 2 for ρr = 1, ρr = 2 and ρr ≥ 2.78 could be approximately 20.9 ≈ 1.87, 21.8 ≈ 3.48 and 22.5 ≈ 5.66, respectively. 21

(a)

(b)

Figure 2: (a) The exact solution y ∗ (t) = t1.9 − 1 for (7.6). (b) The error yN (t) − y ∗ (t) for (7.6) with m = 3, r = 1, ρ = 4, N = 64.

In the case m = 3 for ρr = 1, ρr = 2, ρ = 3 and ρr ≥ 3.89 the ratios %N could be approximatively 20.9 ≈ 1.87, 21.8 ≈ 3.48, 22.7 ≈ 6.50 and 23.5 ≈ 11.31, respectively. These values are given in the last row of Tables 7.4 and 7.5. As we can see from these tables the actual convergence rate accords with the estimate (6.5). Even for ρ = 1, 2 and in the case m = 3 also for ρ = 3 the convergence is substantially faster than it is predicted by the estimate (6.5). From the presented examples we see that the proposed algorithms for the numerical solution of nonlinear fractional initial and boundary value problems of the form (1.1)–(1.2) are effective and with high order accuracy. The obtained experimental convergence rates in Example 1 support the theoretical ones. In particular, it follows from Tables 7.1–7.3 that, in general, the attainable order of global convergence of the proposed algorithms on conditions of Theorem 6.1 cannot be improved. On the other hand, it follows from Tables 7.4 and 7.5 that for some special problems the actual convergence rate of approximate solutions may be better than the theoretical convergence rate given by Theorem 6.1. References [1] K. Diethelm. The Analysis of Fractional Differential Equations, volume 2004 of Lecture Notes in Mathematics. Springer, Berlin, 2010. [2] H. Brunner, A. Pedas, and G. Vainikko. Piecewise polynomial collocation methods for linear volterra integro-differential equations with weakly singular kernels. SIAM J. Numer. Anal., 39:957–982, 2001. [3] A. A. Kilbas, H. M. Srivastava, and J. J. Trujillo. Theory and Applications of Fractional Differential Equations,, volume 204 of North-Holland Mathematics Studies. Elsevier, Amsterdam, 2006. [4] D. Baleanu, K. Diethelm, E. Scalas, and J. J. Trujillo. Fractional Calculus: Models and Numerical Methods. World Scientific Publishing Co. Pte. Ltd, Singapore, 2012. [5] R. Herrmann. Fractional Calculus: An Introduction for Physicists. World Scientific Publishing Co. Pte. Ltd, Hackensack, second edition edition, 2014. [6] F. Mainardi. Fractional Calculus and Waves in Linear Viscoelasticity. Imperial College Press, London, 2010. [7] I. Podlubny. Fractional Differential Equations. Academic Press, San Diego, 1999. [8] R. P. Agarwal, M. Benchohra, and S. Hamani. A survey on existence results for boundary value

22

[9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32]

problems of nonlinear fractional differential equations and inclusions. Acta Appl. Math, 109:973–1033, 2010. G. Vainikko. Multidimensional Weakly Singular Integral Equations, volume 1549 of Lecture Notes in Mathematics. Springer, Berlin, 1993. H. Brunner. Collocation Methods for Volterra Integral and Related Functional Equations, volume 15 of Cambridge Monographs on Applied and Computational Mathematics. Cambridge University Press, Cambridge, 2004. H. Brunner, A. Pedas, and G. Vainikko. The piecewise polynomial collocation method for nonlinear weakly singular volterra equations. Math. Comp., pages 1079–1095, 1999. A. Pedas and E. Tamme. On the convergence of spline collocation methods for solving fractional differential equations. J. Comput. Appl. Math., 235:3502–3514, 2011. A. Pedas and E. Tamme. Spline collocation methods for linear multi-term fractional differential equations. J. Comput. Appl. Math., 236:167–176, 2011. A. Pedas and E. Tamme. Piecewise polynomial collocation for linear boundary value problems of fractional differential equations. J. Comput. Appl. Math., 236:3349–3359, 2012. A. Pedas and E. Tamme. Numerical solution of nonlinear fractional differential equations by spline collocation methods. J. Comput. Appl. Math., 255:216–230, 2014. A. Pedas and E. Tamme. Spline collocation for nonlinear fractional boundary value problems. Appl. Math. Comput., 244:502–513, 2014. A. Pedas, E. Tamme, and M. Vikerpuur. Piecewise polynomial collocation for a class of fractional integro-diffrential equations. In C. Constanda and A. Kirsch, editors, Integral Methods in Science and Engineering, pages 471–482. Springer International Publishing, Switzerland, 2015. N. Kopteva and M. Stynes. An efficient collocation method for a caputo two-point boundary value problem. BIT Numer. Math., 55:1105–1123, 2015. A. H. Bhrawy, T. M. Taha, and J. A. T. Machado. A review of operational matrices and spectral techniques for fractional calculus. Nonlinear Dynam., 81:1023–1052, 2015. H. Brunner and P. J. van der Houven. The Numerical Solution of Volterra Equations. North-Holland, Amsterdam, 1986. H. Brunner. Nonpolynomial spline collocation for volterra equations with weakly singular kernels. SIAM J. Numer. Anal., 20:1106–1119, 1983. Y. Cao, T. Herdman, and Y. Xu. A hybrid collocation method for volterra integral equations with weakly singular kernels. SIAM J. Numer. Anal., 41:364–381, 2003. M. Huang and Y. Xu. Superconvergence of the iterated hybrid collocation method for weakly singular volterra integral equations. J. Integral Equations Appl., 18:83–116, 2006. X. Ma and C. Huang. Numerical solution of fractional integro-differential equations by a hybrid collocation method. Appl. Math Comput., 219:6750–6760, 2013. N. J. Ford, M. L. Morgado, and M. Rebelo. Nonpolynomial collocation approximation of solutions to fractional differential equations,. Fract. Calc. Appl. Anal., 16:874–891, 2013. N. J. Ford, M. L. Morgado, and M. Rebelo. A nonpolynomial collocation method for fractional terminal value problems. J. Comput. Appl. Math., 275:392–402, 2015. T. Diogo, S. McKee, and T. Tang. Collocation methods for second-kind volterra integral equations with weakly singular kernels,. Proc. Roy. Soc. Edinburgh, 124A:199–210, 1994. P. Baratella and A. P. Orsi. A new approach to the numerical solution of weakly singular volterra integral equations. J. Comput. Appl. Math., 163:401–418, 2004. A. Pedas and G. Vainikko. Smoothing transformation and piecewise polynomial collocation for weakly singular volterra integral equations. Computing, 73:271–293, 2004. M. Kolk and A. Pedas. Numerical solution of volterra integral equations with weakly singular kernels which may have a boundary singularity. Math. Model. Anal., 14:79–89, 2009. M. Kolk, A. Pedas, and G. Vainikko. High-order methods for volterra integral equations with general weak singularities,. Numer. Funct. Anal. Optim., 30:1002–1024, 2009. J. Zhao, J. Xiao, and N. J. Ford. Collocation methods for fractional integro-differential equations with

23

weakly singular kernels,. Numer. Algorithms, 65:723–743, 2014. [33] F. Ghoreishi and P. Mokhtary. Spectral collocation method for multi-order fractional differential equations. Int. J. Comput. Methods, 11(05):1350072, 2014. [34] M. Kolk, A. Pedas, and E. Tamme. Modified spline collocation for linear fractional differential equations. J. Comput. Appl. Math., 283:28–40, 2015. [35] M. Kolk, A. Pedas, and E. Tamme. Smoothing transformation and spline collocation for linear fractional boundary value problems. Appl. Math. Comput., 283:234–250, 2016. [36] G. Vainikko. Approximative methods for nonlinear equations (two approaches to the convergence problem). Nonlin. Anal., 2:647–687, 1978.

24