# Cubic Spline Collocation Method for Fractional Differential Equations

Authors: Shui-Ping Yang; Ai-Guo Xiao
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

We discuss the cubic spline collocation method with two parameters for solving the initial value problems (IVPs) of fractional differential equations (FDEs). Some results of the local truncation error, the convergence, and the stability of this method for IVPs of FDEs are obtained. Some numerical examples verify our theoretical results.


---

## Body

## 1. Introduction

In the recent past years, the use of fractional differential equations (FDEs) has gained considerable popularity in some fields (e.g., nonlinear oscillation of earthquake [1], fluid-dynamic traffic model [2], material viscoelastic theory and physics [3–6], etc.). Based on these requirements, the numerical approach has become very important to solve FDEs and analyze the experimental data which are described in a fractional way.The numerical approaches to FDEs have been recently studied by numerous authors [6–15]. However, the state of art is far less advanced for general fractional order differential equations. In the recent years, spline collocation methods including wavelet methods have been successfully applied to some initial value problems (IVPs) and initial-boundary value problems of FDEs [16–23]. Pedas and Tamme [16, 18] discussed spline collocation methods for some classes of IVPs of linear multiterm fractional differential equations and obtained the corresponding convergence results. Li et al. [19] used the higher-order piecewise interpolation for the fractional integral and fractional derivatives, proposed a higher-order algorithm based on the Simpson method for FDEs, and gave the error and stability results. Li [21] applied the cubic B-spline wavelet collocation method to FDEs by the approach that such problems are converted into a system of algebraic equations which is suitable for computer programming.It is worth notice that the spline collocation methods for the IVPs of FDEs are often achieved by their solving of the equivalent IVPs of integral equations with weakly singular kernels. In this paper, the cubic spline collocation method is designed to solve directly the IVPs of general FDEs. And the corresponding theoretical results of the local truncation error, the convergence, and the stability of the cubic spline collocation method for the IVPs of general FDEs are given.This paper is organized as follows. In Section2, we propose the spline cubic collocation method for solving the IVPs of FDEs. In Section 3, the corresponding theoretical results of the convergence and the stability are also given. In Section 4, the theoretical results are identified by some numerical examples.In the following text, we first recall the basic definitions of fractional calculus [6]. Usually, aCDtα denotes the Caputo fractional derivative of order α as
(1)aCDtαy(t)={In-αDny(t)=1Γ(n-α)∫at‍(t-τ)n-α-1dndτny(τ)dτ,t>a,ifn-1<α<n,dndtny(t),ifα=n,
where Dn is the classical differential operator of order n, y(t) is n times continuously differentiable, and Iα denotes the integral operator of order α as
(2)Iαy(t)={1Γ(α)∫at‍(t-τ)α-1y(τ)dτ,t>a,ifα>0,I,ifα=0,
where I is the identity operator.As is well known, there are some different definitions of fractional operator except the Caputo fractional derivative. From a theoretical point of view, the most natural approach is the Riemann-Liouville definition as(3)aRDtαy(t)=DnIn-ay(t)=1Γ(n-α)dndtn∫at‍(t-τ)n-α-1y(τ)dτ,t>a,n-1<α≤n.
The relationship between the Caputo definition and the Riemann-Liouville definition can be given by the following Lemma.Lemma 1 (see [5, 6]).
Letα>0 and n=⌈α⌉. Assume that f(t) is n times continuously differentiable, t0CDtαf and t0RDtαf exist; then
(4)t0RDtαf(t)=t0CDtαf(t)+∑k=0n-1‍fk(t0)Γ(k-α+1)(t-t0)k-α.


## 2. Cubic Spline Collocation Method for FDEs

Consider the initial value problem(5)0CDtαy(t)=f(t,y(t)),t∈I=[0,T],0<α≤1,y(0)=y0,
where f:[0,T]×R→R is a given continuous mapping and satisfies the Lipschitz condition
(6)|f(t,y)-f(t,z)|≤L|y-z|,∀t∈I=[0,T],y,z∈R.
It is well known that the IVP (5) is equivalent to the problem [24]
(7)y(t)=y(0)+1Γ(α)∫0t‍(t-τ)α-1f(τ,y(τ))dτ,y(0)=y0.Letti=ih(i=0,1,…,N,h=T/N,N>0) be the grid points of the uniform partition of [0,T] into the subintervals Ii=[ti-1,ti],i=1,2,…,N. On each subinterval Ii, the cubic C1-spline function S(t,h) can be represented as
(8)S(t,h)=Si(t,h),
where
(9)Si(t,h)=(1+2t-ti-1h)(ti-th)2Si-1(0)+(t-ti-1)(ti-t)2h3Si-1(1)+(1+2ti-th)(t-ti-1h)2Si(0)-(t-ti-1)2(ti-t)h3Si(1),
and S(0,h)=S0(0)=y0, S0(1)=hy′(0).Letξ=(t-ti-1)/h∈[0,1], ξ-=1-ξ=(ti-t)/h. We can write (9) as
(10)Si(t,h)=Si(ti-1+ξh,h)=ξ-2(2ξ+1)Si-1(0)+ξ-2ξSi-1(1)+ξ2(2ξ-+1)Si(0)-ξ2ξ-Si(1),
where Si(0)=S(ti,h)≈y(ti),Si(1)=hS′(ti,h)≈hy′(ti),i=0,1,…,N. Obviously, S(t,h)∈C1[0,T], and S(t,h)≈y(t), for all t∈[0,T].Remark 2.
In this algorithm, the initial valuey′(0) must be provided, but it does not need to be given in problem. That is, if we approximate y′(0), then the approximation of y′(0) is given by using the low-order methods. In this paper, we obtained approximation of y′(0) by using one-order BDF method. That is,
(11)y(h-)≈y~(h-)=y0+h-αΓ(α+1)f(0,y0).
Then
(12)y′(0)≈y~(h-)-y(0)h-,
where h- is a tiny stepsize.By using the collocation conditions in each subintervalIi and (5), we have
(13)0CDtαS(t,h)|t=ti-1+ϕh=f(ti-1+ϕh,S(ti-1+ϕh,h)),ϕ∈{ϕ1,ϕ2},
where ϕ1,ϕ2 are the collocation points, and 0<C0≤ϕ1<ϕ2≤1, C0 is a constant. And according to (10),   we have
(14)(S(ti-1+ϕ1,h)S(ti-1+ϕ2,h))=(ϕ-12(2ϕ1+1)ϕ-12ϕ1ϕ-22(2ϕ2+1)ϕ-22ϕ2)(Si-1(0)Si-1(1))+(ϕ12(2ϕ-1+1)-ϕ12ϕ-1ϕ22(2ϕ-2+1)-ϕ22ϕ-2)(Si(0)Si(1))=C(ϕ1,ϕ2)(Si-1(0)Si-1(1))+D(ϕ1,ϕ2)(Si(0)Si(1)),
where ϕ-j=1-ϕj, i=1,2.Using the definition of the Caputo fractional derivative in the form (1) and using (13),  we have
(15)1Γ(1-α)[∫ti-1ti-1+ϕjh‍S′(τ,h)dτ(ti-1+ϕjh-τ)α+∑n=1i-1‍∫tn-1tn‍S′(τ,h)dτ(ti-1+ϕjh-τ)α]=f(ti-1+ϕjh,S(ti-1+ϕjh,h)),j=1,2.
Through calculation, we can get
(16)-6ϕj2-α(α-3+2ϕ)Si(0)+2ϕj2-α(-3+α+3ϕ)Si(1)+6ϕj2-α(α-3+2ϕ)Si-1(0)+ϕj1-α(-12ϕj+4ϕjα+6ϕj2+6-5α+α2)Si-1(1)+∑n=1i-1‍[-6φ2(n,ϕj)Sn-1(0)+(φ0(n,ϕj)-φ1(n,ϕj)-3φ2(n,ϕj))Sn-1(1)+6φ2(n,ϕj)Si-1(0)+(φ1(n,ϕj)-3φ2(n,ϕj))Sn(1)]=Γ(4-α)hαf(ti-1+ϕjh,S(ti-1+ϕjh)),j=1,2,
where
(17)φ0(n,ϕj)=(i-n+ϕj)1-α(2-α)(3-α)-(i-n+ϕj-1)1-α(3-α),φ1(n,ϕj)=-(i-n+ϕj-1)1-α(2-α)(3-α)-(i-n+ϕj-1)2-α(3-α)+(i-n+ϕj)2-α(3-α),φ2(n,ϕj)=φ1(n,ϕj)+(i-n+ϕj-1)1-α(2-α)(3-α)+2(i-n+ϕj-1)2-α(3-α)+2(i-n+ϕj-1)3-α-2(i-n+ϕj)3-α.
In order to obtain the corresponding iteration formula, we denote(18)A~=(-6ϕ12-α(α-3+2ϕ1)2ϕ12-α(-3+α+3ϕ1)-6ϕ22-α(α-3+2ϕ2)2ϕ22-α(-3+α+3ϕ2)),B~=(6ϕ12-α(α-3+2ϕ1)ϕ11-α(-12ϕ1+4ϕα+6ϕ12+6-5α+α2)6ϕ22-α(α-3+2ϕ2)ϕ21-α(-12ϕ2+4ϕα+6ϕ22+6-5α+α2)),H~(i)=(∑n=0i-1‍[-6φ2(n,ϕ1)Sn-1(0)+(φ0(n,ϕ1)-φ1(n,ϕ1)-3φ2(n,ϕ1))Sn-1(1)+6φ2(n,ϕ1)Sn(0)+(φ1(n,ϕ1)-3φ2(n,ϕ1))Sn(1)]∑n=0i-1‍[-6φ2(n,ϕ2)Sn-1(0)+(φ0(n,ϕ2)-φ1(n,ϕ2)-3φ2(n,ϕ2))Sn-1(1)+6φ2(n,ϕ2)Sn(0)+(φ1(n,ϕ2)-3φ2(n,ϕ2))Sn(1)]).Obviously, it follows from the inequalitiesC0≤ϕ1<ϕ2≤1 that the matrix A~ is invertible, and
(19)S_i=-A~-1B~S_i-1+Γ(4-α)hαA~-1f_i-A~-1H~(i),
where
(20)S_i=(Si(0),Si(1))T,f_i=(fi-1+ϕ1,fi-1+ϕ2)T=(f(ti-1+ϕ1,S(ti-1+ϕ1,h)),f(ti-1+ϕ2,S(ti-1+ϕ2,h)))T.Let(21)A=-A~-1B~,B=A~-1,H(i)=-A~-1H~(i).
Then, it follows that
(22)S_i=AS_i-1+hαΓ(4-α)Bf_i+H(i).
For (22), we can obtain their numerical solutions by using the Newton iterative method.In addition, applying the collocation conditions in each subintervalIi and using (7), we have
(23)S-(t,h)|t=ti-1+ϕh=y0+1Γ(α)∫0ti-1+ϕh(ti-1+ϕh-τ)α-1f(τ,S-(τ,h))dτ,‍ϕ∈{ϕ1,ϕ2},
where ϕ1,ϕ2 are the collocation points, 0<C0≤ϕ1<ϕ2≤1, and
(24)S-(t,h)=(1+2t-ti-1h)(ti-th)2S-i-1(0)+(t-ti-1)(ti-t)2h3S-i-1(1)+(1+2ti-th)(t-ti-1h)2S-i(0)-(t-ti-1)2(ti-t)h3S-i(1),
where S-(0,h)=S-0(0)=y0,S-0(1)=hy′(0),S-i(0)=S-(ti,h)≈y(ti),S-i(1)=hS-′(ti,h)≈hy′(ti), i=1,2,…,N.In this paper, for convenience, the cubic spline collocation method based on the direct discretization with (13) is called direct spline collocation method (DSCM for short); and the cubic spline collocation method based on the fractional order integral equations (7) and (27) is called indirect spline collocation method (ISCM for short).


## 3. Theoretical Results

Lemma 3 (see [25]).
Suppose thatzn≥0,n=0,1,…,N, and satisfy
(25)zn≤hαc0∑i=0n-1‍(n-j)α-1zj+c1,n=0,1,…,N,
where 0<α≤1,c1≥0,c0>0 is independent of h>0. Then
(26)zn≤c1Eα(c0Γ(α)(nh)α),n=0,1,…,N,
where Eα is a Mittag-Leffler function.Lemma 4.
Assume that the functiong(t)∈C1[0,T]. Then, for all ϕ∈[0,1],t∈[0,T-h], one has
(27)∫tt+h(s-τ)α-1g(τ)dτ=g(t+ϕh)∫tt+h‍(s-τ)α-1dτ+o(h),
where 0<α≤1, 0<s≤T is a constant and satisfies t+h=s or s-(t+h)≥C0h, C0∈(0,1], o(·) is infinitesimal of higher order.Proof.
Using the Taylor formula, we have(28)g(τ)=g(t+ϕh)+g′(t+ϕh)(τ-t-ϕh)+o(h),∀τ∈[t,t+h].
Applying the integral mean value theorem, we can obtain
(29)∫tt+h‍(s-τ)α-1g(τ)dτ=∫tt+h‍(s-τ)α-1×(g(t+ϕh)+g′(t+ϕh)×(τ-t-ϕh)+o(h)g′)dτ=g(t+ϕh)∫tt+h‍(s-τ)α-1dτ+g′(t+ϕh)∫tt+h‍(s-τ)α-1(τ-t-ϕh)dτ+o(h)=g(t+ϕh)∫tt+h‍(s-τ)α-1dτ+g′(t+ϕh)(ξ-t-ϕh)∫tt+h‍(s-τ)α-1dτ+o(h)=g(t+ϕh)∫tt+h‍(s-τ)α-1dτ+g′(t+ϕh)(ξ-t-ϕh)×[1α((s-t)α-(s-t-h)α)]+o(h),
where ξ∈[t,t+h]. Since g(t)∈C1[0,T], g′(t) is bounded on [0,T]. We have
(30)∫tt+h‍(s-τ)α-1g(τ)dτ=g(t+ϕh)∫tt+h‍(s-τ)α-1dτ+O(h)1α[(s-t)α-(s-t-h)α]+o(h).
Consider the following.(i)
Ifs=t+h, then
(31)∫tt+h‍(s-τ)α-1g(τ)dτ=g(t+ϕh)∫tt+h‍(s-τ)α-1dτ+O(h)1α[(t+h-t)α]+o(h)=g(t+ϕh)∫tt+h‍(s-τ)α-1dτ+O(h)1αhα+o(h)=g(t+ϕh)∫tt+h‍(s-τ)α-1dτ+O(h1+α)+o(h)=g(t+ϕh)∫tt+h‍(s-τ)α-1dτ+o(h).(ii)
IfC0h≤s-(t+h), then
(32)(s-τ)α-1≤(s-(t+h))α-1≤C0α-1hα-1,∀τ∈[t,t+h],1α[(s-t)α-(s-t-h)α]=∫tt+h‍(s-τ)α-1dτ≤∫tt+h‍C0α-1hα-1dτ=C0α-1hα.
Thus,(33)1α[(s-t)α-(s-t-h)α]=o(1).
By means of (30), Lemma 4 is established obviously.Theorem 5 (the local truncation error).
If the analytical solution of the problem (5) y(t)∈C4[0,T], then the local truncation error of DSCM is O(h4).Proof.
Applying the Taylor expansion and the definition ofS(t,h), we have
(34)y(t)=(1+2t-ti-1h)(ti-th)2y(ti-1)y(t)=+(t-ti-1)(ti-t)2h3hy′(ti-1)y(t)=+(1+2ti-th)(t-ti-1h)2y(ti)y(t)=-(t-ti-1)2(ti-t)h3hy′(ti)+O((t-ti-1)4)y(t)=y-(t)+O((t-ti-1)4),y′(t)=y-′(t)+O((t-ti-1)3),∀t∈[ti-1,ti].
From the definition of the Caputo fractional derivative, we obtain
(35)1Γ(1-α)[∫ti-1ti-1+ϕjh‍(y-′(τ)+O((τ-ti-1)3))dτ(ti-1+ϕjh-τ)α1Γ(1-α)+∑n=1i-1‍∫tn-1tn‍y′(τ)dτ(ti-1+ϕjh-τ)α]=f(ti-1+ϕjh,y-(ti-1+ϕjh)+O(ϕjh)4),j=1,2.
That is,
(36)1Γ(1-α)∫ti-1ti-1+ϕjh‍y-′(τ)dτ(ti-1+ϕjh-τ)α=f(ti-1+ϕjh,y-(ti-1+ϕjh)+O(ϕjh)4)-1Γ(1-α)[∑n=1i-1‍∫tn-1tn‍y′(τ)dτ(ti-1+ϕjh-τ)α-1Γ(1-α)+∫ti-1ti-1+ϕjh‍O((τ-ti-1)3)dτ(ti-1+ϕjh-τ)α],
where j=1,2. Moreover,
(37)(y(ti)hy′(ti))=A(y(ti-1)hy′(ti-1))+hαΓ(4-α)B×(f(ti-1+ϕ1h,y-(ti-1+ϕ1h)+O(ϕ1h)4)f(ti-1+ϕ2h,y-(ti-1+ϕ2h)+O(ϕ2h)4))-hαA~-1{∫ti-1ti-1+ϕ1h‍O((τ-ti-1)3)dτ(ti-1+ϕ1h-τ)α(1Γ(1-α)∫ti-1ti-1+ϕ1h‍O((τ-ti-1)3)dτ(ti-1+ϕ1h-τ)α1Γ(1-α)∫ti-1ti-1+ϕ2h‍O((τ-ti-1)3)dτ(ti-1+ϕ2h-τ)α)(∑n=1i-1‍∫tn-1tn‍y′(τ)dτ(ti-1+ϕ1h-τ)α∑n=1i-1‍∫tn-1tn‍y′(τ)dτ(ti-1+ϕ2h-τ)α)+(1Γ(1-α)∫ti-1ti-1+ϕ1h‍O((τ-ti-1)3)dτ(ti-1+ϕ1h-τ)α1Γ(1-α)∫ti-1ti-1+ϕ2h‍O((τ-ti-1)3)dτ(ti-1+ϕ2h-τ)α)}.
SupposingS(t,h)=y(t), for all t∈[0,ti-1], and substituting y(t) into (13), we have
(38)(y(ti)hy′(ti))=A(y(ti-1)hy′(ti-1))+hαΓ(4-α)B×(f(ti-1+ϕ1h,S(ti-1+ϕ1h,h))f(ti-1+ϕ2h,S(ti-1+ϕ2h,h)))+H-(i)+Ti,(39)H-(i)=-hαA~-1(∑n=1i-1‍∫tn-1tn‍S′(τ,h)dτ(ti-1+ϕ1h-τ)α∑n=1i-1‍∫tn-1tn‍S′(τ,h)dτ(ti-1+ϕ2h-τ)α)=-hαA~-1(∑n=1i-1‍∫tn-1tn‍y′(τ)dτ(ti-1+ϕ1h-τ)α∑n=1i-1‍∫tn-1tn‍y′(τ)dτ(ti-1+ϕ2h-τ)α),Si(t,h)=(1+2t-ti-1h)(ti-th)2y(ti-1)+(t-ti-1)(ti-t)2h3hy′(ti-1)+(1+2ti-th)(t-ti-1h)2y(ti)-(t-ti-1)2(ti-t)h3hy′(ti),
where Ti is the local truncation error.
Combining (38), (37) with the Lipschitz condition yields
(40)|Ti|≤hαLΓ(4-α)∥B∥(|O((ϕ1h)4)||O((ϕ2h)4)|)+hα∥A~-1∥(1Γ(2-α)(ϕ1h)1-α|O(h3)|1Γ(2-α)(ϕ2h)1-α|O(h3)|),∥Ti∥≤C^h4,
where C^ is an appropriate positive constant. And the norm ∥·∥ is the 1-norm of the matrix (·) in this paper.To prove the convergence of DSCM for the problem (5), we first give the following lemmas.Lemma 6.
If the analytical solution of the problem (5) y(t)∈C4[0,T] and the matrixes C(ϕ1,ϕ2) and D(ϕ1,ϕ2) in (14) satisfy that D-1(ϕ1,ϕ2) exists, and
(41)∥D-1(ϕ1,ϕ2)C(ϕ1,ϕ2)∥<1,
then the numerical solutions of ISCM for problem (5) S-(t,h) satisfy that limh→0S-(t,h)=y(t).Proof.
For (23), we have
(42)S-(t,h)|t=ti+ϕjh=y0+1Γ(α)∫0ti+ϕjh‍(ti+ϕjh-τ)α-1f(τ,S-(τ,h))dτ,=y0+1Γ(α)∫0ti+ϕjh‍(ti+ϕjh-τ)α-1222221j=1,2,S-(0,h)=y0,S-(ti+ϕjh,h)=y0+1Γ(α)[∑n=1i‍∫titi+ϕjh(ti+ϕjh-τ)α-1f(τ,S-(τ,h))dτ‍+∑n=1i‍∫tn-1tn‍(ti+ϕjh-τ)α-1f(τ,S-(τ,h))dτ]j=1,2.
It follows from (7) that the analytical solution y(t) of the problem (5) satisfies
(43)y(ti+ϕjh)=y(0)+1Γ(α)[∑n=1i‍∫titi+ϕjh‍(ti+ϕjh-τ)α-1f(τ,y(τ))dτ+∑n=1i‍∫tn-1tn‍(ti+ϕjh-τ)α-1f(τ,y(τ))dτ],j=1,2.
Moreover, (44) and (45) yield
(44)S-(ti+ϕjh,h)-y(ti+ϕjh)=1Γ(α){∫titi+ϕjh‍(ti+ϕjh-τ)α-1×[f(τ,S-(τ,h))-f(τ,y(τ))]dτ+∑n=1i‍∫tn-1tn‍(ti+ϕjh-τ)α-1×[f(τ,S-(τ,h))-f(τ,y(τ))]dτ∫titi+ϕjh‍}=hαΓ(α)×{∫0ϕj‍(ϕj-ξ)α-1[f(ti+ξh,S-(ti+ξh,h))-f(ti+ξh,y(ti+ξh))S-]dξ+∑n=0i-1‍∫01‍(i+ϕj-n-ξ)α-1×[f(tn+ξh,S-(tn+ξh,h))-f(tn+ξh,y(tn+ξh))(tn+ξh,S-(tn+ξh,h))]dξ∫0ϕj‍},j=1,2.
By using the Lipschitz condition, we have(45)|S-(ti+ϕjh,h)-y(ti+ϕjh)|=|1Γ(α){∫titi+ϕjh‍(ti+ϕjh-τ)α-1×[f(τ,S-(τ,h))-f(τ,y(τ))]dτ+∑n=1i‍∫tn-1tn‍(ti+ϕjh-τ)α-1×[f(τ,S-(τ,h))-f(τ,y(τ))]dτ∑n=1i‍}|≤1Γ(α){∑n=1i‍∫titi+ϕjh‍(ti+ϕjh-τ)α-1×|f(τ,S-(τ,h))-f(τ,y(τ))|dτ+∑n=1i‍∫tn-1tn‍(ti+ϕjh-τ)α-1×|f(τ,S-(τ,h))-f(τ,y(τ))|dτ∑n=1i‍}≤LΓ(α){∫titi+ϕjh‍(ti+ϕjh-τ)α-1|S-(τ,h)-y(τ)|dτ+∑n=1i‍∫tn-1tn‍(ti+ϕjh-τ)α-1|S-(τ,h)-y(τ)|dτ}=hαLΓ(α){∫0ϕj‍(ϕj-τ)α-1|S-(ti+τh,h)-y(ti+τh)|dτ+∑n=0i-1‍∫01‍(i-n+ϕj-τ)α-1×|S-(tn+τh,h)-y(tn+τh)|dτ∑n=0i-1‍},j=1,2.
It follows from the definition ofS-(t,h) that
(46)(y(ti+hu)y(ti+hv))=(u-2(2u+1)u-2uv-2(2v+1)v-2v)(y(ti)hy′(ti))+(u2(2u-+1)-u2u-v2(2v-+1)-v2v-)(y(ti+1)hy′(ti+1))=C(u,v)(yihyi′)+D(u,v)(yi+1hyi+1′)+O(h4),∀u,v∈[0,1],
where u-=1-u, v-=1-v. Obviously, for (14) and (46), the matrixes C(u,v) and D(u,v) are continuous on [0,1]×[0,1]; thus, we have
(47)(S-(ti+hu,h)-y(ti+hu)S-(ti+hv,h)-y(ti+hv))=(u-2(2u+1)u-2vv-2(2v+1)v-2v)(S-i(0)-y(ti)S-i(1)-hy′(ti))+(u2(2u-+1)-u2u-v2(2v-+1)-v2v-)(S-i+1(0)-y(ti+1)S-i+1(1)-hy′(ti+1))=C(u,v)(S-i(0)-y(ti)S-i(1)-hy′(ti))+D(u,v)(S-i+1(0)-y(ti+1)S-i+1(1)-hy′(ti+1))+O(h4),∀u,v∈[0,1].
Therefore, there exist λ1>0, λ2>0 such that
(48)λ1=sup0≤u<v≤1∥C(u,v)∥,λ2=sup0≤u<v≤1∥D(u,v)∥,λ=λ1+λ2.
Obviously,(49)∫0ϕj‍(ϕj-ξ)α-1dξ=-1α(ϕj-ξ)α|0ϕj=1α(ϕj)α≤1α,∀ϕj∈[C0,1].
Applying the integral mean theorem, we have
(50)∫01‍(i-n+ϕj-ξ)α-1dξ=(i-n+ϕj-θ)α-1(0≤θ≤1)≤(i-n)α-1(1+ϕj-θi-n)α-1≤(i-n)α-1(1+ϕj-1i-n)α-1.
It follows from the inequalities 0<C0≤ϕ1<ϕ2≤1 and 0≤n≤i-1 that
(51)∫01‍(i-n+ϕj-ξ)α-1dξ≤(i-n)α-1(1+ϕj-1i-n)α-1≤(i-n)α-1(1+ϕj-1)α-1≤(i-n)α-1C0α-1,j=1,2.
Combining (45), (48), and (49) with (51) and defining
(52)ϵi=(|S-i(0)-y(ti)||S-i(1)-hy′(ti)|),ϵi′=(|S-(ti+ϕ1,h)-y(ti+ϕ1)||S-(ti+ϕ2,h)-y(ti+ϕ2)|),Xi=max0≤n≤i+1{∥ϵn∥},
where ti+ϕj=ti+ϕjh, j=1,2, we can obtain
(53)∥ϵi′∥≤hαL[1Γ(α)(λ1∥ϵi∥+λ2∥ϵi+1∥)α+C^1h4Γ(α+1)((i+1)h)α+1Γ(α)∑n=0i-1‍C0α-1(i-n)α-1×(λ1∥ϵn∥+λ2∥ϵn+1∥)1Γ(α)]≤hαL[1Γ(α+1)λXi+1Γ(α)∑n=0i-1‍C0α-1×(i-n)α-1λXn∑n=0i-1‍]+C^2h4,
where C^1,C^2 are some positive constants. If the matrix D(ϕ1,ϕ2) is invertible, then it follows from (47) that we have
(54)(S-i+1(0)-y(ti+1)S-i+1(1)-hy′(ti+1))=-D-1(ϕ1,ϕ2)C(ϕ1,ϕ2)(S-i(0)-y(ti)S-i(1)-hy′(ti))+D-1(ϕ1,ϕ2)(S-(ti+hϕ1,h)-y(ti+hϕ1)S-(ti+hϕ2,h)-y(ti+hϕ2))+O(h4).
Hence,
(55)∥ϵi+1∥≤∥D-1(ϕ1,ϕ2)C(ϕ1,ϕ2)∥∥ϵi∥+∥D-1(ϕ1,ϕ2)∥∥ϵi′∥+C^3h4≤∥D-1(ϕ1,ϕ2)C(ϕ1,ϕ2)∥Xi+∥D-1(ϕ1,ϕ2)∥∥ϵi′∥+C^3h4,
where C^3 is an appropriate positive constant. Take ∥ϵℓ+1∥=max0≤n≤i+1∥ϵn∥=Xi, 0≤ℓ≤i. Then Xℓ=max0≤n≤ℓ+1∥ϵn∥=∥ϵℓ+1∥, Xi=Xℓ, and
(56)Xℓ=∥ϵℓ+1∥≤∥D-1(ϕ1,ϕ2)C(ϕ1,ϕ2)∥∥ϵℓ∥+∥D-1(ϕ1,ϕ2)∥∥ϵℓ′∥+C^3h4≤∥D-1(ϕ1,ϕ2)C(ϕ1,ϕ2)∥Xℓ+∥D-1(ϕ1,ϕ2)∥∥ϵℓ′∥+C^3h4.
Since∥D-1(ϕ1,ϕ2)C(ϕ1,ϕ2)∥<1, we have
(57)(1-∥D-1(ϕ1,ϕ2)C(ϕ1,ϕ2)∥)Xℓ≤∥D-1(ϕ1,ϕ2)∥∥ϵℓ′∥+C^3h4,Xℓ≤C1∥ϵℓ′∥+C^4h4,
where C1=∥D-1(ϕ1,ϕ2)∥/(1-∥D-1(ϕ1,ϕ2)C(ϕ1,ϕ2)∥), C^4 is an appropriate positive constant. According to (53), we obtain
(58)Xℓ≤C1hαL[∑n=0ℓ-1‍1Γ(α+1)λXℓ+1Γ(α)×∑n=0ℓ-1‍C0α-1(ℓ-n)α-1λXn1Γ(α)∑n=0ℓ-1‍]+C^5h4,(1-C1hαLλΓ(α+1))Xℓ≤C1hαL(1Γ(α)∑n=0ℓ-1‍C0α-1(ℓ-n)α-1λXn)+C^5h4,
where C^5>0 is an appropriate constant. Obviously, there exist C~∈(0,1) and h1>0 such that C1hαLλ/Γ(α+1)≤C~<1 for h<h1. Thus,
(59)Xℓ≤C1hαLλ(1-C~)(1Γ(α)∑n=0ℓ-1‍C0α-1(ℓ-n)α-1Xn)+C^6h4.
Let C=C1C0α-1λ/(1-C~)Γ(α); then
(60)∥ϵi+1∥≤Xℓ≤Chα(∑n=0ℓ-1‍(ℓ-n)α-1Xn)+C^6h4,
where C^6>0 is an appropriate constant. Applying Lemma 3, we get
(61)∥ϵi+1∥≤Xℓ≤C^6h4Eα(CΓ(α)(ℓh)α)≤C^6h4Eα(CΓ(α)Tα).
By using the convergence of the Mittag-Leffler function Eα(z) [6], we have
(62)Xℓ⟶0,∥ϵi+1∥⟶0(h⟶0).
By means of (24), (47), we obtain
(63)sup0≤t≤T∥S-(t,h)-y(t)∥≤max0≤i≤N-1{sup0≤u<v≤1∥C(u,v)∥∥ϵi∥+sup0≤u<v≤1∥D(u,v)∥∥ϵi+1∥}+C^7h4≤(λ1+λ2)Eα(CΓ(α)Tα)C^6h4+C^7h4≤[C^7+(λ1+λ2)C^6Eα(CΓ(α)Tα)]h4,
where C^7 is a positive constant. Consequently, we have
(64)limh→0S-(t,h)=y(t).Remark 7.
(
1
) If ϕ1≠ϕ2, then
(65)D(ϕ1,ϕ2)=(ϕ12(2ϕ-1+1)-ϕ12ϕ-1ϕ22(2ϕ-2+1)-ϕ22ϕ-2)
is invertible, and(66)D-1(ϕ1,ϕ2)=(--1+ϕ2ϕ12(-ϕ2+ϕ1)-1+ϕ1ϕ22(-ϕ2+ϕ1)--3+2ϕ2ϕ12(-ϕ2+ϕ1)-3+2ϕ1ϕ22(-ϕ2+ϕ1)).
(
2
)
Consider (67)D-1(ϕ1,ϕ2)C(ϕ1,ϕ2)=(-(ϕ1ϕ2+ϕ1+ϕ2)(-1+ϕ2)(-1+ϕ1)ϕ12ϕ22-(-1+ϕ2)(-1+ϕ1)ϕ1ϕ22ϕ12-3ϕ1+2ϕ1ϕ2-3ϕ2+2ϕ22ϕ12ϕ22-ϕ1ϕ2-2ϕ1+3-2ϕ2ϕ1ϕ2).When∥D-1(ϕ1,ϕ2)C(ϕ1,ϕ2)∥<1, the range of values ϕ1,ϕ2 can be obtained by using Algorithm 8, which can be shown in Figures 1 and 2, respectively.The range of valuesϕ1,ϕ2 when ∥D-1C∥1<1 (a) or ∥D-1C∥2<1 (b).
(a)
(b)Figure 2
The range of valuesϕ1,ϕ2 when ∥D-1C∥∞<1.Algorithm 8.
For0≤ϕ1≤1
Forϕ1<ϕ2≤1
if  ∥D-1(ϕ1,ϕ2)C(ϕ1,ϕ2)∥<1
recordϕ1,ϕ2;
end
end
end
plot(ϕ1,ϕ2).Lemma 9.
If the analytical solution of the problem (5) y(t)∈C4[0,T], the spline functions S(t,h) and S-(t,h) defined as (9) and (24) on the same grids are the numerical solutions of the problem (5) obtained by DSCM and ISCM, respectively; then we have the following conclusions.(i)
Iflimh→0S-(t,h)=y(t), then limh→0S(t,h)=y(t), where y(t) is the analytical solution of the problem (5).(ii)
Assume thatyS(t,h),zS(t,h) and yS-(t,h),zS-(t,h) are the numerical solutions obtained by DSCM and ISCM for the problem
(68)0CDtαy(t)=f(t,y(t)),t∈I=[0,T],0<α≤1,y(0)=y0,and the perturbed problem
(69)0CDtαz(t)=f(t,z(t)),t∈I=[0,T],0<α≤1,z(0)=z0,
where f satisfies the Lipschitz condition (6), respectively. Then there exists the constant h0>0 such that, for all h<h0, when |yS-(tn,h)-zS-(tn,h)|≤C1|y0-z0|, there exists the constant C>0 such that |yS(tn,h)-zS(tn,h)|≤C|y0-z0|; here C1>0 is a constant, n=1,2,…,N.Proof.
Firstly, we give the proof of the conclusion (i). Applying DSCM, we have(70)0CDtαS(t,h)|t=ti-1+ϕh=f(ti-1+ϕh,S(ti-1+ϕh,h)),i=1,2…,N.
Since S(t,h)∈C1[0,T], we obtain
(71)S(ti-1+ϕh,h)=y0+1Γ(α)∫0ti-1+ϕh‍(ti-1+ϕh-τ)α-10CDταS(τ,h)dτ=y0+1Γ(α)[∑n=1i-1‍∫tn-1tn‍(ti-1+ϕh-τ)α-10CDταS(τ,h)dτ+∫ti-1ti-1+ϕh‍(ti-1+ϕh-τ)α-10CDταS(τ,h)dτ∑n=1i-1‍].
It follows from Lemma 4 that
(72)S(ti-1+ϕh,h)=y0+1Γ(α)[∑n=1i-1‍(0CDταS(τ,h)|τ=tn-1+ϕ×∫tn-1tn‍(ti-1+ϕh-τ)α-1dτ+o(h))i+0CDταS(τ,h)|τ=ti-1+ϕhi×∫ti-1ti-1+ϕh‍(ti-1+ϕh-τ)α-1dτ+o(h)∑n=1i-1‍∑n=1i-1‍]=y0+1Γ(α)[∑n=1i-1‍0CDταS(τ,h)|τ=tn-1+ϕi×∫tn-1tn(ti-1+ϕh-τ)α-1dτi+0CDταS(τ,h)|τ=ti-1+ϕhi×∫ti-1ti-1+ϕh(ti-1+ϕh-τ)α-1dτ∑n=1i-1‍]+o(1),
where o(1)→0  as  h→0.
By means of (70) and (72), we obtain
(73)S(ti-1+ϕh,h)=y0+1Γ(α)[∑n=1i-1f(tn-1+ϕ,S(tn-1+ϕ,h))i×∫tn-1tn(ti-1+ϕh-τ)α-1dτi+f(ti-1+ϕh,S(ti-1+ϕh,h))i×∫ti-1ti-1+ϕh(ti-1+ϕh-τ)α-1dτ∑]+o(1).
Applying Lemma 4, we get
(74)S(ti-1+ϕh,h)=y0+1Γ(α)[∑n=1i-1(∫tn-1tn(ti-1+ϕh-τ)α-1×f(τ,S(τ,h))dτ+o(h)∫)+∫ti-1ti-1+ϕh‍(ti-1+ϕh-τ)α-1×f(τ,S(τ,h))dτ+o(h)∑n=1i-1(∫tn-1tn(ti-1+ϕh-τ)α-1]+o(1)=y0+1Γ(α)[∑n=1i-1∫tn-1tn(ti-1+ϕh-τ)α-1f(τ,S(τ,h))dτ+∫ti-1ti-1+ϕh‍(ti-1+ϕh-τ)α-1×f(τ,S(τ,h))dτ∑]+o(1),
where there exists a constant C^8>0, γ>0 such that |o(1)|=C^8hγ. Based on the definitions of S(t,h),S-(t,h), it follows that
(75)|S(ti-1+ϕjh,h)-S-(ti-1+ϕjh,h)|≤|1Γ(α){∫ti-1ti-1+ϕjh(ti-1+ϕjh-τ)α-1×[f(τ,S(τ,h))-f(τ,S-(τ,h))]dτ+∑n=1i-1‍∫tn-1tn‍(ti-1+ϕjh-τ)α-1×[f(τ,S(τ,h))-f(τ,S-(τ,h))]dτ∫ti-1ti-1+ϕjh(ti-1+ϕjh-τ)α-1∫}|+C^8hγ≤1Γ(α){∫ti-1ti-1+ϕjh‍(ti-1+ϕjh-τ)α-1×|f(τ,S(τ,h))-f(τ,S-(τ,h))|dτ+∑n=1i-1∫tn-1tn(ti-1+ϕjh-τ)α-1×|f(τ,S(τ,h))-f(τ,S-(τ,h))|dτ∫ti-1ti-1+ϕjh‍}+C^8hγ≤LΓ(α){∫ti-1ti-1+ϕjh‍(ti-1+ϕjh-τ)α-1×|S(τ,h)-S-(τ,h)|dτ+∑n=1i-1‍∫tn-1tn‍(ti-1+ϕjh-τ)α-1×|S(τ,h)-S-(τ,h)|dτ∫ti-1ti-1+ϕjh‍}+C^8hγ,j=1,2.
Based on the similar proof to that of Lemma6, it follows from Lemma 3 and the inequality ∥D-1(ϕ1,ϕ2)C(ϕ1,ϕ2)∥<1 that
(76)∥ϵi∥≤C^9hγEα(CΓ(α)Tα),i=1,2,…,N,
where
(77)ϵi=(|Si(0)-S-i(0)||Si(1)-S-i(1)|),
and C is defined in Lemma 6, C^9>0 is a constant. By means of the convergence of the Mittag-Leffler function Eα(z) [6], we get
(78)∥ϵi∥⟶0,h⟶0;limh→0S(t,h)=limh→0S-(t,h).
From Lemma 6, we have
(79)limh→0S-(t,h)=y(t),limh→0S(t,h)=limh→0S-(t,h)=y(t).
Now, we give the proof of the conclusion (ii) of Lemma6. According to the conclusion (i), we get
(80)|yS(tn,h)-zS(tn,h)|=|yS(tn,h)-yS-(tn,h)+yS-(tn,h)-zS-(tn,h)+zS-(tn,h)-zS(tn,h)|≤|yS(tn,h)-yS-(tn,h)|+|yS-(tn,h)-zS-(tn,h)|+|zS-(tn,h)-zS(tn,h)|≤max0≤i≤N|yE(ti,h)|+max0≤i≤N|zE(ti,h)|+C1|y0-z0|,n=1,2,…,N,
where yE(t,h)=yS(t,h)-yS-(t,h), zE(t,h)=zS(t,h)-zS-(t,h).
According to the proof of the conclusion (i), there existsh0>0 such that, for all h<h0,
(81)max0≤i≤N|yE(ti,h)|+max0≤i≤N|zE(ti,h)|≤|y0-z0|.
Then,
(82)|yS(tn,h)-zS(tn,h)|≤(1+C1)|y0-z0|,n=1,2,…,N.
Obviously, take C=1+C1; hence, the conclusion (ii) is right.Theorem 10 (convergence of DSCM).
If the analytical solution of the problem (5) y(t)∈C4[0,T] and the matrixes C(ϕ1,ϕ2) and D(ϕ1,ϕ2) of (14) satisfy that D-1(ϕ1,ϕ2) exists, and
(83)∥D-1(ϕ1,ϕ2)C(ϕ1,ϕ2)∥<1,
then DSCM is convergent, and limh→0S(t,h)=y(t).Proof.
This result follows directly from Lemmas6 and 9.Now, we consider the stability of DSCM.Theorem 11 (stability of DSCM).
If the analytical solution of the problem (5) y(t)∈C4[0,T] and the matrixes C(ϕ1,ϕ2) and D(ϕ1,ϕ2) in (14) satisfy that D-1(ϕ1,ϕ2) exists and
(84)∥D-1(ϕ1,ϕ2)C(ϕ1,ϕ2)∥<1,
then ISCM is stable; that is if ISCM is, applied to solve the problem (68) and the perturbed problem (69), respectively, then, for all 0<h<h1 (h1>0 is a constant), one has(i)
(85)∥ϵi∥≤C|y0-z0|,i=1,2,…,N,
where
(86)ϵi=(|yS-i(0)-zS-i(0)||yS-i(1)-zS-i(1)|),C>0 is a constant, {yS-i(0),yS-i(1)} are the numerical solutions of the problem (68), and {zS-i(0),zS-i(1)} are the numerical solutions of the perturbed problem (69), and(ii)
(87)sup0≤t≤T∥yS-(t,h)-zS-(t,h)∥≤C′|y0-z0|,
where C′>0 is a constant, yS-(t,h) is the numerical solution of the problem (68), and zS-(t,h) is the numerical solution of the problem (69).
Moreover, DSCM is stable.Proof.
According to Lemma9, if ISCM is stable, then DSCM is stable. In the following text, we give the proof of stability of ISCM.
By using ISCM to solve the problem (68) and the problem (69), we have
(88)yS-(t,h)|t=ti+ϕjh=y0+1Γ(α)∫0ti+ϕjh‍(ti+ϕjh-τ)α-1×f(τ,yS-(τ,h))dτ,j=1,2,yS-(0,h)=y0,zS-(t,h)|t=ti+ϕjh=z0+1Γ(α)∫0ti+ϕjh‍(ti+ϕjh-τ)α-1×f(τ,zS-(τ,h))dτ,j=1,2,zS-(0,h)=z0.
The previous equations can be written as
(89)yS-(ti+ϕjh,h)=y0+1Γ(α)[∫titi+ϕjh‍(ti+ϕjh-τ)α-1f(τ,yS-(τ,h))dτ+∑n=1i‍∫tn-1tn‍(ti+ϕjh-τ)α-1f(τ,yS-(τ,h))dτ],j=1,2,(90)zS-(ti+ϕjh,h)=z0+1Γ(α)[∑n=1i‍∫titi+ϕjh‍(ti+ϕjh-τ)α-1f(τ,zS-(τ,h))dτ+∑n=1i‍∫tn-1tn‍(ti+ϕjh-τ)α-1f(τ,zS-(τ,h))dτ],j=1,2.
Thus,
(91)yS-(ti+ϕjh,h)-zS-(ti+ϕjh,h)=(y0-z0)+1Γ(α)×{∫titi+ϕjh‍(ti+ϕjh-τ)α-1×[f(τ,yS-(τ,h))-f(τ,zS-(τ,h))]dτ+∑n=1i‍∫tn-1tn‍(ti+ϕjh-τ)α-1×[f(τ,yS-(τ,h))-f(τ,zS-(τ,h))]dτ∫titi+ϕjh‍}=(y0-z0)+hαΓ(α)×{∫0ϕj‍(ϕj-ξ)α-1×[f(ti+ξh,yS-(ti+ξh,h))-f(ti+ξh,zS-(ti+ξh,h))]dξ+∑n=0i-1‍∫01‍(i+ϕj-n-ξ)α-1×[f(tn+ξh,yS-(tn+ξh,h))-f(tn+ξh,zS-(tn+ξh,h))]dξ∑n=0i-1‍},j=1,2.
Applying the Lipschitz condition, we get
(92)|yS-(ti+ϕjh,h)-zS-(ti+ϕjh,h)|≤|y0-z0|+|1Γ(α){∫titi+ϕjh‍(ti+ϕjh-τ)α-1×[f(τ,yS-(τ,h))-f(τ,zS-(τ,h))]dτ+∑n=1i‍∫tn-1tn‍(ti+ϕjh-τ)α-1×[f(τ,yS-(τ,h))-f(τ,zS-(τ,h))]dτ∫titi+ϕjh‍}|≤|y0-z0|+1Γ(α){∫titi+ϕjh‍(ti+ϕjh-τ)α-1×|f(τ,yS-(τ,h))-f(τ,zS-(τ,h))|dτ+∑n=1i‍∫tn-1tn‍(ti+ϕjh-τ)α-1×|f(τ,yS-(τ,h))-f(τ,zS-(τ,h))|dτ∫titi+ϕjh‍}≤|y0-z0|+LΓ(α){∫titi+ϕjh‍(ti+ϕjh-τ)α-1×|yS-(τ,h)-zS-(τ,h)|dτ+∑n=1i‍∫tn-1tn‍(ti+ϕjh-τ)α-1×|yS-(τ,h)-zS-(τ,h)|dτ∫titi+ϕjh‍}=|y0-z0|+hαLΓ(α){∫0ϕj‍(ϕj-τ)α-1×|yS-(ti+τh,h)-zS-(ti+τh,h)|dτ+∑n=0i-1‍∫01‍(i-n+ϕj-τ)α-1×|yS-(tn+τh,h)-zS-(tn+τh,h)|dτ∫0ϕj‍},
where j=1,2. Note that the definitions of yS-(t,h),zS-(t,h). Let
(93)λ1=sup0≤u<v≤1∥C(u,v)∥,λ2=max0≤u<v≤1∥D(u,v)∥,λ=λ1+λ2,ϵi=(|yS-i(0)-zS-i(0)||yS-i(1)-zS-i(1)|),Xi=max0≤n≤i+1{∥ϵn∥},ϵi′=(|yS-(ti+ϕ1,h)-zS-(ti+ϕ1,h)||yS-(ti+ϕ2,h)-zS-(ti+ϕ2,h)|),
where ti+ϕj=ti+ϕjh,j=1,2.
It follows from (92), (93), (49), and (51) that
(94)∥ϵi′∥≤|y0-z0|+hαL[1Γ(α)(λ1∥ϵi∥+λ2∥ϵi+1∥)α+1Γ(α)∑n=0i-1‍C0α-1(i-n)α-1×(λ1∥ϵn∥+λ2∥ϵn+1∥)1Γ(α)]≤|y0-z0|+hαL[∑n=0i-1‍1Γ(α+1)(λXi)+1Γ(α)×∑n=0i-1‍C0α-1(i-n)α-1(λXn)∑n=0i-1‍1Γ(α+1)]≤|y0-z0|+hαL[∑n=0i-1‍1Γ(α+1)λXi+1Γ(α)×∑n=0i-1‍C0α-1(i-n)α-1λXn1Γ(α+1)].
Take∥ϵℓ+1∥=max0≤n≤i+1{∥ϵn∥}=Xi; hence, Xi=Xℓ,1≤ℓ≤i. According to the fact that D(ϕ1,ϕ2) is invertible and ∥D-1(ϕ1,ϕ2)C(ϕ1,ϕ2)∥<1, we have
(95)Xℓ=∥ϵℓ+1∥≤∥D-1(ϕ1,ϕ2)C(ϕ1,ϕ2)∥∥ϵℓ∥+∥D-1(ϕ1,ϕ2)∥∥ϵℓ′∥≤∥D-1(ϕ1,ϕ2)C(ϕ1,ϕ2)∥Xℓ+∥D-1(ϕ1,ϕ2)∥∥ϵℓ′∥,(1-∥D-1(ϕ1,ϕ2)C(ϕ1,ϕ2)∥)Xℓ≤∥D-1(ϕ1,ϕ2)∥∥ϵℓ′∥.
Moreover,(96)Xℓ≤C1∥ϵℓ′∥,
where C1=∥D-1(ϕ1,ϕ2)∥/(1-∥D-1(ϕ1,ϕ2)C(ϕ1,ϕ2)∥). By means of (94), we have
(97)Xℓ≤C1|y0-z0|+C1hαL[∑n=0ℓ-1‍1Γ(α+1)λXℓ+1Γ(α)∑n=0ℓ-1‍C0α-1(ℓ-n)α-1λXn],(98)(1-C1hαLλΓ(α+1))Xℓ≤C1|y0-z0|+C1Lhα(1Γ(α)∑n=0ℓ-1‍C0α-1(ℓ-n)α-1λXn).
Obviously, there existC~∈(0,1) and h1>0 such that C1hαLλ/Γ(α+1)≤C~<1 as h<h1. Then
(99)Xℓ≤C1(1-C~)(|y0-z0|+C1hαLΓ(α)∑n=0ℓ-1‍C0α-1(ℓ-n)α-1Xn).
Let C-=C1C0α-1Lλ/(1-C~)Γ(α). We obtain
(100)Xℓ≤C1(1-C~)|y0-z0|+C-hα∑n=0ℓ-1‍(ℓ-n)α-1Xn.
Applying Lemma3 yields
(101)Xℓ≤C1(1-C~)|y0-z0|Eα(C-Γ(α)(ℓh)α)≤C1(1-C~)|y0-z0|Eα(C-Γ(α)Tα).
Using the convergence of Mittag-Leffler function Eα(z) and denoting
(102)C=C1(1-C~)Eα(C-Γ(α)Tα),
we obtain
(103)∥ϵi+1∥≤Xi=Xℓ≤C|y0-z0|.
Moreover, from the definitions ofyS-(t,h) and zS-(t,h), we have
(104)(yS-(ti+hu,h)-zS-(ti+hu,h)yS-(ti+hv,h)-zS-(ti+hv,h))=(u-2(2u+1)u-2vv-2(2v+1)v-2v)(yS-i(0)-zS-i(0)yS-i(1)-zS-i(1))+(u2(2u-+1)-u2u-v2(2v-+1)-v2v-)(yS-i+1(0)-zS-i+1(0)yS-i+1(1)-zS-i+1(1))=C(u,v)(yS-i(0)-zS-i(0)yS-i(1)-zS-i(1))+D(u,v)(yS-i+1(0)-zS-i+1(0)yS-i+1(1)-zS-i+1(1)),∀u,v∈[0,1].
By means of (103), we get
(105)sup0≤t≤T∥yS-(t,h)-zS-(t,h)∥≤max0≤i≤N-1{sup0≤u<v≤1∥C(u,v)∥∥ϵi∥+sup0≤u<v≤1∥D(u,v)∥∥ϵi+1∥}≤(λ1+λ2)C|y0-z0|≤C′|y0-z0|,
where C′=(λ1+λ2)C; that is, ISCM is stable. According to Lemma 9, DSCM is stable.


## 4. Illustrative Examples

In order to demonstrate our theoretical results, we apply DSCM to the problem (5) and present some numerical examples in this section. Let
(106)∥error∥L1=h∑i=1N‍|S(ti,h)-y(ti)|,∥error∥L2=(h∑i=1N‍|S(ti,h)-y(ti)|2)1/2,orderL1=log1/2∥error∥L1(h/2)∥error∥L1(h),orderL2=log1/2∥error∥L2(h/2)∥error∥L2(h),
where ∥error∥(h) is the error function which depends on the stepsize h.Example 1.
Consider the initial value problem [10, 21]
(107)0CDtαy(t)=-y(t)+t2+2t2-αΓ(3-α),t∈I=[0,1],0<α<1,y(0)=0.
The exact solution is y(t)=t2.For different values ofα∈(0,1), the numerical solutions for problem (107) are obtained by using DSCM, fast wavelet collocation method (FWCM for short) reported in [21], and the method reported in [10]. When α=0.5, the time stepsize h=1/160; the absolute errors of DSCM, FWCM, and the method reported in [10] are shown in Table 1. When α=0.75, h=1/160; the absolute errors of these methods are shown in Table 2. The numerical solutions and the exact solution are shown in Figure 3. From the numerical results, the results obtained by DSCM are better than by FWCM and the method reported in [10] in terms of accuracy if the exact solution is sufficiently smooth. Therefore, DSCM is a valid method in solving fractional differential equation.Table 1
The absolute errors of DSCM (ϕ1=0.85,ϕ2=0.99) and the methods reported in [10, 21]; α=0.5.
t
Error by DSCM
Error by [21]
Error by [10]
0.031250
2.168404
×
1
0
-
18
0.059585760872
×
1
0
-
12
0.043283783288
×
1
0
-
4
0.093750
1.734723
×
1
0
-
18
0.043017672757
×
1
0
-
12
0.115137021309
×
1
0
-
4
0.187500
6.938894
×
1
0
-
18
0.053908266740
×
1
0
-
12
0.134008681836
×
1
0
-
4
0.281250
0.000000
×
1
0
-
19
0.024466539905
×
1
0
-
12
0.145485346092
×
1
0
-
4
0.375000
2.775558
×
1
0
-
17
0.003941291737
×
1
0
-
12
0.153757162245
×
1
0
-
4
0.468750
0.000000
×
1
0
-
19
0.064837024638
×
1
0
-
12
0.160223410662
×
1
0
-
4
0.562500
0.000000
×
1
0
-
19
0.054012350148
×
1
0
-
12
0.165528139359
×
1
0
-
4
0.656250
5.551115
×
1
0
-
17
0.009825473768
×
1
0
-
12
0.170022170524
×
1
0
-
4
0.750000
1.110223
×
1
0
-
16
0.017985612999
×
1
0
-
12
0.173918021056
×
1
0
-
4
0.843750
0.000000
×
1
0
-
19
0.031863400807
×
1
0
-
12
0.177354280767
×
1
0
-
4
0.937500
1.110223
×
1
0
-
16
0.201283434365
×
1
0
-
12
0.180426383865
×
1
0
-
4Table 2
The absolute errors of DSCM (ϕ1=0.9,ϕ2=0.98) and the methods reported in [10, 21]; α=0.75.
t
Error by DSCM
Error by [21]
Error by [10]
0.03125
7.155734
×
1
0
-
18
0.330128888907
×
1
0
-
13
0.043283783288
×
1
0
-
4
0.093750
1.040834
×
1
0
-
17
0.347308987125
×
1
0
-
13
0.115137021309
×
1
0
-
4
0.187500
4.163336
×
1
0
-
17
0.022967738822
×
1
0
-
13
0.134008681836
×
1
0
-
4
0.281250
9.714451
×
1
0
-
16
0.198174809896
×
1
0
-
13
0.145485346092
×
1
0
-
4
0.375000
2.220446
×
1
0
-
16
0.293931545769
×
1
0
-
13
0.153757162245
×
1
0
-
4
0.468750
3.885781
×
1
0
-
16
0.221489493413
×
1
0
-
13
0.160223410662
×
1
0
-
4
0.562500
7.216450
×
1
0
-
16
0.194289029309
×
1
0
-
13
0.165528139359
×
1
0
-
4
0.656250
9.436896
×
1
0
-
16
0.170974345793
×
1
0
-
13
0.170022170524
×
1
0
-
4
0.750000
1.221245
×
1
0
-
15
0.091038288019
×
1
0
-
13
0.173918021056
×
1
0
-
4
0.843750
1.665335
×
1
0
-
15
0.119904086659
×
1
0
-
13
0.177354280767
×
1
0
-
4
0.937500
2.220446
×
1
0
-
15
0.275335310107
×
1
0
-
13
0.180426383865
×
1
0
-
4Figure 3
The error between the numerical solutions and the exact solution of Example1; α=0.75.Example 2.
Consider the initial value problem [12]
(108)0CDtαy(t)=40320Γ(9-α)t8-α-3Γ(5+α/2)Γ(5-α/2)t4-α/2+94Γ(1+α)+(32tα/2-t4)3-y(t)3/2,t∈I=[0,1],y(0)=0,0<α≤1.
The exact solution is y(t)=t8-3t4+α/2+(9/4)tα. Obviously, y(t)∉C4[0,1]. In order to verify our theoretical results, we modify problem (108) as follows:
(109)0CDtαy(t)=40320Γ(9-α)t8-α-3Γ(7+α/2)Γ(7-α/2)t6-α/2+332Γ(5+α)t4+(32t2+α/2-t4)3-y(t)3/2,t∈I=[0,1],y(0)=0,0<α≤1.
The exact solution is y(t)=t8-3t6+α/2+(9/4)t4+α.Tables3 and 4 list the errors and the error orders of DSCM with α=0.5 and α=0.75 and show that the errors between the numerical solutions and the exact solution are very small, respectively. Figures 4 and 5 illustrate high accuracy of DSCM with ϕ1=0.9, ϕ2=0.98 and show that the error is also very small. All of the numerical results show that DSCM for solving nonlinear FDEs is convergent and the method is robust.Table 3
The error between the numerical solutions and the exact solution of the problem (109); α=0.5, ϕ1=0.88,ϕ2=0.98.
h
∥
error
∥
L
1
Orde
r
L
1
∥
error
∥
L
2
Orde
r
L
2
1
/
10
1.935883
×
1
0
-
5
2.379899
×
1
0
-
5
1
/
20
1.762835
×
1
0
-
6
3.457022411
2.223867
×
1
0
-
6
3.419757934
1
/
40
1.615949
×
1
0
-
7
3.447443874
2.045199
×
1
0
-
7
3.442757379
1
/
80
1.477809
×
1
0
-
8
3.450810893
1.868714
×
1
0
-
8
3.452123534
1
/
160
1.340725
×
1
0
-
9
3.462413611
1.697493
×
1
0
-
9
3.460568256Table 4
The error between the numerical solutions and the exact solution of the problem (109); α=0.75, ϕ1=0.95,ϕ2=0.99.
h
∥
error
∥
L
1
Orde
r
L
1
∥
error
∥
L
2
Orde
r
L
2
1
/
10
5.391748
×
1
0
-
5
6.474353
×
1
0
-
5
1
/
20
5.313670
×
1
0
-
6
3.342972530
6.727347
×
1
0
-
6
3.264483517
1
/
40
5.541285
×
1
0
-
7
3.261416159
7.058210
×
1
0
-
7
3.252663415
1
/
80
5.812611
×
1
0
-
8
3.252962303
7.434173
×
1
0
-
8
3.247058188
1
/
160
6.111707
×
1
0
-
9
3.249539075
7.834071
×
1
0
-
9
3.246338153The error between the numerical solutions and the exact solution of the problem (109); h=0.05, α=0.000001.
(a)
(b)The error between the numerical solutions and the exact solution of the problem (109); α=0.9, h=0.05.
(a)
(b)Example 3.
Consider the initial value problem(110)0CDtαy(t)+ty2(t)=f(t),t∈I=[0,T],0<α≤1,y(0)=θ1
and the perturbed problem
(111)0CDtαz(t)+tz2(t)=f(t),t∈I=[0,T],0<α≤1,z(0)=θ2.We use the cubic spline collocation method to solve the problems (110) and (111), respectively. Selecting f(t)=tsin(t), T=20, θ1=0.1, θ2=0.2, we obtain the numerical results given in Figure 6. When f(t)=tsin(t), T=20, θ1=2.0, θ2=-2.0, the numerical results are given in Figure 7. When f(t)=tet, T=20, θ1=1.0, θ2=2.0, the numerical results are shown in Figure 8.The numerical results of the problems (110) and (111); h=0.5, T=20.0, α=0.8, θ1=0.1, θ2=0.2, ϕ1=0.9, ϕ2=0.98.
(a)
(b)The numerical results of the problems (110) and (111); h=0.5, α=0.8, θ1=2.0, θ2=-2.0, ϕ1=0.9, ϕ2=0.98.
(a)
(b)The numerical results of the problems (110) and (111); h=0.5, α=0.8, θ1=1.0, θ2=2.0, ϕ1=0.9, ϕ2=0.98.
(a)
(b)From these figures, we can see that the absolute errors of the numerical solutions of the problems (110) and (111) decrease and finally tend to 0 as t increases. Thus, we can draw the conclusion that the cubic spline collocation method for nonlinear FDEs is stable. The numerical results verify our theoretical results.


## 5. Conclusion

In this paper, the cubic spline collocation method with two parameters is successfully applied to the IVPs of general FDEs. The result of the local truncation error of this method is given. And the convergence and stability results of the cubic spline collocation method for the fractional order integral equations which is equivalent to the IVPs of general FDEs are obtained. By using the relationship between the numerical solutions from the cubic spline method for the IVPs of general FDEs and the numerical solutions obtained from the cubic spline method for the corresponding equivalent IVPs of fractional order integral equations, we also obtain some results of the convergence and the stability of the method for the IVPs of general FDEs. Some numerical examples successfully verify our theoretical results and show that the given method is efficient.



---

*Source: 864025-2013-08-07.xml*