Cubic Spline Collocation Method for Fractional Differential Equations

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.

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 Section 2, 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, ^𝐶_𝑎𝐷^𝛼_𝑡 denotes the Caputo fractional derivative of order 𝛼 as ^𝐶_𝑎𝐷^𝛼_𝑡𝑦(𝑡)=𝐼^𝑛−𝛼𝐷^𝑛𝑦(𝑡)=1Γ(𝑛−𝛼)^𝑡_𝑎(𝑡−𝜏)^{𝑛−𝛼−1}𝑑^𝑛𝑑𝜏^𝑛𝑦(𝜏)𝑑𝜏,𝑡>𝑎,if𝑑𝑛−1<𝛼<𝑛,_𝑛𝑑𝑡^𝑛𝑦(𝑡),if𝛼=𝑛, (1) where 𝐷^𝑛 is the classical differential operator of order 𝑛, 𝑦(𝑡) is 𝑛 times continuously differentiable, and 𝐼^𝛼 denotes the integral operator of order 𝛼 as 𝐼^𝛼𝑦(𝑡)=1Γ(𝛼)^𝑡_𝑎(𝑡−𝜏)^𝛼−1𝑦(𝜏)𝑑𝜏,𝑡>𝑎,if𝛼>0,𝐼,if𝛼=0, (2) where 𝐼 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 ^𝑅_𝑎𝐷^𝛼_𝑡𝑦(𝑡)=𝐷^𝑛𝐼^𝑛−𝑎𝑦(𝑡)=1Γ()𝑑𝑛−𝛼𝑑𝑡^𝑛𝑛𝑡_𝑎(𝑡−𝜏)^{𝑛−𝛼−1}𝑦(𝜏)𝑑𝜏,𝑡>𝑎,𝑛−1<𝛼≤𝑛. (3) 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 𝑛=⌈𝛼⌉. Assume that 𝑓(𝑡) is 𝑛 times continuously differentiable, ^𝐶_𝑡0𝐷^𝛼_𝑡𝑓 and ^𝑅_𝑡0𝐷^𝛼_𝑡𝑓 exist; then ^𝑅_𝑡0𝐷^𝛼_𝑡𝑓(𝑡)=^𝐶_𝑡0𝐷^𝛼_𝑡𝑓(𝑡)+^𝑛−1_𝑘=0𝑓^𝑘𝑡₀Γ(𝑘−𝛼+1)𝑡−𝑡₀^𝑘−𝛼. (4)

2. Cubic Spline Collocation Method for FDEs

Consider the initial value problem ^𝐶₀𝐷^𝛼_𝑡𝑦(𝑡)(𝑡)[]𝑦(0)=𝑓𝑡,𝑦,𝑡∈𝐼=0,𝑇,0<𝛼≤1,=𝑦₀, (5) where 𝑓∶[0,𝑇]×𝑅→𝑅 is a given continuous mapping and satisfies the Lipschitz condition 𝑓𝑡,𝑦−𝑓(𝑡,𝑧)≤𝐿𝑦−𝑧,∀𝑡∈𝐼=[0,𝑇],𝑦,𝑧∈𝑅. (6) It is well known that the IVP (5) is equivalent to the problem [24] 𝑦(𝑡)(0)+1=𝑦Γ(𝛼)^𝑡₀(𝑡−𝜏)^𝛼−1𝑓(𝜏)𝑦(0)𝜏,𝑦𝑑𝜏,=𝑦₀. (7)

Let 𝑡_𝑖=𝑖ℎ(𝑖=0,1,…,𝑁,ℎ=𝑇/𝑁,𝑁>0) be the grid points of the uniform partition of [0,𝑇] into the subintervals 𝐼_𝑖=[𝑡_𝑖−1,𝑡_𝑖],𝑖=1,2,…,𝑁. On each subinterval 𝐼_𝑖, the cubic 𝐶¹-spline function 𝑆(𝑡,ℎ) can be represented as 𝑆(𝑡,ℎ)=𝑆_𝑖(𝑡,ℎ), (8) where 𝑆_𝑖()=𝑡,ℎ1+2𝑡−𝑡_𝑖−1ℎ𝑡_𝑖−𝑡ℎ²𝑆⁽⁰⁾_𝑖−1+𝑡−𝑡_𝑖−1𝑡_𝑖−𝑡²ℎ³𝑆⁽¹⁾_𝑖−1+𝑡1+2^𝑖−𝑡ℎ𝑡−𝑡_𝑖−1ℎ²𝑆_𝑖⁽⁰⁾−𝑡−𝑡_𝑖−1²𝑡_𝑖−𝑡ℎ³𝑆_𝑖⁽¹⁾, (9) and 𝑆(0,ℎ)=𝑆₀⁽⁰⁾=𝑦₀, 𝑆₀⁽¹⁾=ℎ𝑦^′(0).

Let 𝜉=(𝑡−𝑡_𝑖−1)/ℎ∈[0,1], 𝜉=1−𝜉=(𝑡_𝑖−𝑡)/ℎ. We can write (9) as 𝑆_𝑖(𝑡,ℎ)=𝑆_𝑖𝑡_𝑖−1=+𝜉ℎ,ℎ𝜉²(2𝜉+1)𝑆⁽⁰⁾_𝑖−1+𝜉²𝜉𝑆⁽¹⁾_𝑖−1+𝜉²2𝜉+1𝑆_𝑖⁽⁰⁾−𝜉²𝜉𝑆_𝑖⁽¹⁾, (10) where 𝑆_𝑖⁽⁰⁾=𝑆(𝑡_𝑖,ℎ)≈𝑦(𝑡_𝑖),𝑆_𝑖⁽¹⁾=ℎ𝑆^′(𝑡_𝑖,ℎ)≈ℎ𝑦^′(𝑡_𝑖),𝑖=0,1,…,𝑁. Obviously, 𝑆(𝑡,ℎ)∈𝐶¹[0,𝑇], and 𝑆(𝑡,ℎ)≈𝑦(𝑡), for all 𝑡∈[0,𝑇].

Remark 2. In this algorithm, the initial value 𝑦^′(0) must be provided, but it does not need to be given in problem. That is, if we approximate 𝑦^′(0), then the approximation of 𝑦^′(0) is given by using the low-order methods. In this paper, we obtained approximation of 𝑦^′(0) by using one-order BDF method. That is, 𝑦ℎ≈𝑦ℎ=𝑦₀+ℎ^𝛼Γ()𝑓𝛼+10,𝑦₀. (11) Then 𝑦^′(0)≈𝑦ℎ−𝑦(0)ℎ, (12) where ℎ is a tiny stepsize.

By using the collocation conditions in each subinterval 𝐼_𝑖 and (5), we have ^𝐶₀𝐷^𝛼_𝑡𝑆(𝑡,ℎ)_{𝑡=𝑡𝑖−1+𝜙ℎ}=𝑓𝑡_𝑖−1+𝜙ℎ,𝑆𝑡_𝑖−1+𝜙ℎ,ℎ,𝜙∈𝜙₁,𝜙₂, (13) where 𝜙₁,𝜙₂ are the collocation points, and 0<𝐶₀≤𝜙₁<𝜙₂≤1, 𝐶₀ is a constant. And according to (10),   we have 𝑆𝑡_{𝑖−1+𝜙1}𝑆𝑡,ℎ_{𝑖−1+𝜙2}=,ℎ𝜙²₁2𝜙₁+1𝜙²₁𝜙₁𝜙²₂2𝜙₂+1𝜙²₂𝜙₂𝑆⁽⁰⁾_𝑖−1𝑆⁽¹⁾_𝑖−1+𝜙²₁2𝜙₁+1−𝜙²₁𝜙₁𝜙²₂2𝜙₂+1−𝜙²₂𝜙₂𝑆_𝑖⁽⁰⁾𝑆_𝑖⁽¹⁾=𝐶𝜙₁,𝜙₂𝑆⁽⁰⁾_𝑖−1𝑆⁽¹⁾_𝑖−1+𝐷𝜙₁,𝜙₂𝑆_𝑖⁽⁰⁾𝑆_𝑖⁽¹⁾, (14) where 𝜙_𝑗=1−𝜙_𝑗, 𝑖=1,2.

Using the definition of the Caputo fractional derivative in the form (1) and using (13),  we have 1Γ(1−𝛼)^{𝑡𝑖−1+𝜙𝑗ℎ}_𝑡𝑖−1𝑆^′(𝜏,ℎ)𝑑𝜏𝑡_𝑖−1+𝜙_𝑗ℎ−𝜏^𝛼+^𝑖−1_𝑛=1^𝑡𝑛_𝑡𝑛−1𝑆^′(𝜏,ℎ)𝑑𝜏𝑡_𝑖−1+𝜙_𝑗ℎ−𝜏^𝛼=𝑓𝑡_𝑖−1+𝜙_𝑗ℎ,𝑆𝑡_𝑖−1+𝜙_𝑗ℎ,ℎ,𝑗=1,2. (15) Through calculation, we can get −6𝜙_𝑗^2−𝛼𝛼−3+2𝜙𝑆_𝑖⁽⁰⁾+2𝜙_𝑗^2−𝛼−3+𝛼+3𝜙𝑆_𝑖⁽¹⁾+6𝜙_𝑗^2−𝛼𝛼−3+2𝜙𝑆⁽⁰⁾_𝑖−1+𝜙_𝑗^1−𝛼−12𝜙_𝑗+4𝜙_𝑗𝛼+6𝜙²_𝑗+6−5𝛼+𝛼²𝑆⁽¹⁾_𝑖−1+^𝑖−1_𝑛=1−6𝜑₂𝑛,𝜙_𝑗𝑆⁽⁰⁾_𝑛−1+𝜑₀𝑛,𝜙_𝑗−𝜑₁𝑛,𝜙_𝑗−3𝜑₂𝑛,𝜙_𝑗𝑆⁽¹⁾_𝑛−1+6𝜑₂𝑛,𝜙_𝑗𝑆⁽⁰⁾_𝑖−1+𝜑₁𝑛,𝜙_𝑗−3𝜑₂𝑛,𝜙_𝑗𝑆_𝑛⁽¹⁾=Γ(4−𝛼)ℎ^𝛼𝑓𝑡_𝑖−1+𝜙_𝑗ℎ,𝑆𝑡_𝑖−1+𝜙_𝑗ℎ,𝑗=1,2, (16) where 𝜑₀𝑛,𝜙_𝑗=𝑖−𝑛+𝜙_𝑗^1−𝛼()−2−𝛼)(3−𝛼𝑖−𝑛+𝜙_𝑗−1^1−𝛼(),𝜑3−𝛼₁𝑛,𝜙_𝑗=−𝑖−𝑛+𝜙_𝑗−1^1−𝛼()−2−𝛼)(3−𝛼𝑖−𝑛+𝜙_𝑗−1^2−𝛼()+3−𝛼𝑖−𝑛+𝜙_𝑗^2−𝛼(),𝜑3−𝛼₂𝑛,𝜙_𝑗=𝜑₁𝑛,𝜙_𝑗+𝑖−𝑛+𝜙_𝑗−1^1−𝛼(2−𝛼)(3−𝛼)+2𝑖−𝑛+𝜙_𝑗−1^2−𝛼(3−𝛼)+2𝑖−𝑛+𝜙_𝑗−1^3−𝛼−2𝑖−𝑛+𝜙_𝑗^3−𝛼. (17) In order to obtain the corresponding iteration formula, we denote𝐴=−6𝜙₁^2−𝛼𝛼−3+2𝜙₁2𝜙₁^2−𝛼−3+𝛼+3𝜙₁−6𝜙₂^2−𝛼𝛼−3+2𝜙₂2𝜙₂^2−𝛼−3+𝛼+3𝜙₂,𝐵=6𝜙₁^2−𝛼𝛼−3+2𝜙₁𝜙₁^1−𝛼−12𝜙₁+4𝜙𝛼+6𝜙²₁+6−5𝛼+𝛼²6𝜙₂^2−𝛼𝛼−3+2𝜙₂𝜙₂^1−𝛼−12𝜙₂+4𝜙𝛼+6𝜙²₂+6−5𝛼+𝛼²,𝐻(𝑖)=^{𝑖−1𝑛=0}−6𝜑₂𝑛,𝜙₁𝑆⁽⁰⁾_𝑛−1+𝜑₀𝑛,𝜙₁−𝜑₁𝑛,𝜙₁−3𝜑₂𝑛,𝜙₁𝑆⁽¹⁾_𝑛−1+6𝜑₂𝑛,𝜙₁𝑆_𝑛⁽⁰⁾+𝜑₁𝑛,𝜙₁−3𝜑₂𝑛,𝜙₁𝑆_𝑛⁽¹⁾_{𝑖−1𝑛=0}−6𝜑₂𝑛,𝜙₂𝑆⁽⁰⁾_𝑛−1+𝜑₀𝑛,𝜙₂−𝜑₁𝑛,𝜙₂−3𝜑₂𝑛,𝜙₂𝑆⁽¹⁾_𝑛−1+6𝜑₂𝑛,𝜙₂𝑆_𝑛⁽⁰⁾+𝜑₁𝑛,𝜙₂−3𝜑₂𝑛,𝜙₂𝑆_𝑛⁽¹⁾. (18)

Obviously, it follows from the inequalities 𝐶₀≤𝜙₁<𝜙₂≤1 that the matrix 𝐴 is invertible, and 𝑆_𝑖=−𝐴⁻¹𝐵𝑆_𝑖−1+Γ(4−𝛼)ℎ^𝛼𝐴⁻¹𝑓_𝑖−𝐴⁻¹𝐻(𝑖), (19) where 𝑆_𝑖=𝑆_𝑖⁽⁰⁾,𝑆_𝑖^(1)𝑇,𝑓_𝑖=𝑓_{𝑖−1+𝜙1},𝑓_{𝑖−1+𝜙2}^𝑇=𝑓𝑡_{𝑖−1+𝜙1},𝑆𝑡_{𝑖−1+𝜙1},ℎ,𝑓𝑡_{𝑖−1+𝜙2},𝑆𝑡_{𝑖−1+𝜙2},ℎ^𝑇. (20)

Let 𝐴=−𝐴⁻¹𝐵,𝐵=𝐴⁻¹,𝐻(𝑖)=−𝐴⁻¹𝐻(𝑖). (21) Then, it follows that 𝑆_𝑖=𝐴𝑆_𝑖−1+ℎ^𝛼Γ(4−𝛼)𝐵𝑓_𝑖+𝐻(𝑖). (22) For (22), we can obtain their numerical solutions by using the Newton iterative method.

In addition, applying the collocation conditions in each subinterval 𝐼_𝑖 and using (7), we have 𝑆(𝑡,ℎ)_{𝑡=𝑡𝑖−1+𝜙ℎ}=𝑦₀+1Γ(𝛼)^𝑡𝑖−1₀^+𝜙ℎ𝑡_𝑖−1+𝜙ℎ−𝜏^𝛼−1𝑓𝜏,𝑆()𝜙𝜏,ℎ𝑑𝜏,𝜙∈₁,𝜙₂, (23) where 𝜙₁,𝜙₂ are the collocation points, 0<𝐶₀≤𝜙₁<𝜙₂≤1, and 𝑆()=𝑡,ℎ1+2𝑡−𝑡_𝑖−1ℎ𝑡_𝑖−𝑡ℎ²𝑆⁽⁰⁾_𝑖−1+𝑡−𝑡_𝑖−1𝑡_𝑖−𝑡²ℎ³𝑆⁽¹⁾_𝑖−1+𝑡1+2^𝑖−𝑡ℎ𝑡−𝑡_𝑖−1ℎ²𝑆_𝑖⁽⁰⁾−𝑡−𝑡_𝑖−1²𝑡_𝑖−𝑡ℎ³𝑆_𝑖⁽¹⁾, (24) where 𝑆(0,ℎ)=𝑆₀⁽⁰⁾=𝑦₀,𝑆₀⁽¹⁾=ℎ𝑦^′(0),𝑆_𝑖⁽⁰⁾=𝑆(𝑡_𝑖,ℎ)≈𝑦(𝑡_𝑖),𝑆_𝑖⁽¹⁾=ℎ𝑆^′(𝑡_𝑖,ℎ)≈ℎ𝑦^′(𝑡_𝑖), 𝑖=1,2,…,𝑁.

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 that 𝑧_𝑛≥0,𝑛=0,1,…,𝑁, and satisfy 𝑧_𝑛≤ℎ^𝛼𝑐₀^𝑛−1_𝑖=0𝑛−𝑗^𝛼−1𝑧_𝑗+𝑐₁,𝑛=0,1,…,𝑁, (25) where 0<𝛼≤1,𝑐₁≥0,𝑐₀>0 is independent of ℎ>0. Then 𝑧_𝑛≤𝑐₁𝐸_𝛼𝑐₀Γ(𝛼)(𝑛ℎ)^𝛼,𝑛=0,1,…,𝑁, (26) where 𝐸_𝛼 is a Mittag-Leffler function.

Lemma 4. Assume that the function 𝑔(𝑡)∈𝐶¹[0,𝑇]. Then, for all 𝜙∈[0,1],𝑡∈[0,𝑇−ℎ], one has _𝑡^𝑡+ℎ(𝑠−𝜏)^𝛼−1𝑔(𝜏)𝑑𝜏=𝑔𝑡+𝜙ℎ_𝑡^𝑡+ℎ(𝑠−𝜏)^𝛼−1𝑑𝜏+𝑜(ℎ), (27) where 0<𝛼≤1, 0<𝑠≤𝑇 is a constant and satisfies 𝑡+ℎ=𝑠 or 𝑠−(𝑡+ℎ)≥𝐶₀ℎ, 𝐶₀∈(0,1], 𝑜(⋅) is infinitesimal of higher order.

Proof. Using the Taylor formula, we have 𝑔(𝜏)=𝑔𝑡+𝜙ℎ+𝑔^′(ℎ),[].𝑡+𝜙ℎ𝜏−𝑡−𝜙ℎ+𝑜∀𝜏∈𝑡,𝑡+ℎ (28) Applying the integral mean value theorem, we can obtain _𝑡^𝑡+ℎ(𝑠−𝜏)^𝛼−1𝑔(𝜏)=𝑑𝜏_𝑡𝑡+ℎ(𝑠−𝜏)^𝛼−1×𝑔𝑡+𝜙ℎ+𝑔^′×(ℎ)𝑡+𝜙ℎ𝜏−𝑡−𝜙ℎ+𝑜𝑑𝜏=𝑔𝑡+𝜙ℎ_𝑡^𝑡+ℎ(𝑠−𝜏)^𝛼−1𝑑𝜏+𝑔^′𝑡+𝜙ℎ_𝑡^𝑡+ℎ(𝑠−𝜏)^𝛼−1𝜏−𝑡−𝜙ℎ𝑑𝜏+𝑜(ℎ)=𝑔𝑡+𝜙ℎ_𝑡^𝑡+ℎ(𝑠−𝜏)^𝛼−1𝑑𝜏+𝑔^′𝑡+𝜙ℎ𝜉−𝑡−𝜙ℎ_𝑡^𝑡+ℎ(𝑠−𝜏)^𝛼−1𝑑𝜏+𝑜(ℎ)=𝑔𝑡+𝜙ℎ_𝑡^𝑡+ℎ(𝑠−𝜏)^𝛼−1𝑑𝜏+𝑔^′×1𝑡+𝜙ℎ𝜉−𝑡−𝜙ℎ𝛼(𝑠−𝑡)^𝛼−(𝑠−𝑡−ℎ)^𝛼+𝑜(ℎ), (29) where 𝜉∈[𝑡,𝑡+ℎ]. Since 𝑔(𝑡)∈𝐶¹[0,𝑇], 𝑔^′(𝑡) is bounded on [0,𝑇]. We have _𝑡^𝑡+ℎ(𝑠−𝜏)^𝛼−1𝑔(𝜏)𝑑𝜏=𝑔𝑡+𝜙ℎ_𝑡^𝑡+ℎ(𝑠−𝜏)^𝛼−1(ℎ)1𝑑𝜏+𝑂𝛼(𝑠−𝑡)^𝛼−(𝑠−𝑡−ℎ)^𝛼+𝑜(ℎ). (30)
Consider the following. (i) If 𝑠=𝑡+ℎ, then _𝑡^𝑡+ℎ(𝑠−𝜏)^𝛼−1𝑔(𝜏)𝑑𝜏=𝑔𝑡+𝜙ℎ_𝑡^𝑡+ℎ(𝑠−𝜏)^𝛼−1(ℎ)1𝑑𝜏+𝑂𝛼(𝑡+ℎ−𝑡)^𝛼+𝑜(ℎ)=𝑔𝑡+𝜙ℎ_𝑡^𝑡+ℎ(𝑠−𝜏)^𝛼−1(ℎ)1𝑑𝜏+𝑂𝛼ℎ^𝛼+𝑜(ℎ)=𝑔𝑡+𝜙ℎ_𝑡^𝑡+ℎ(𝑠−𝜏)^𝛼−1ℎ𝑑𝜏+𝑂^1+𝛼+𝑜(ℎ)=𝑔𝑡+𝜙ℎ_𝑡^𝑡+ℎ(𝑠−𝜏)^𝛼−1𝑑𝜏+𝑜(ℎ). (31)(ii) If 𝐶₀ℎ≤𝑠−(𝑡+ℎ), then (𝑠−𝜏)^𝛼−1≤(𝑠−(𝑡+ℎ))^𝛼−1≤𝐶₀^𝛼−1ℎ^𝛼−1[],1,∀𝜏∈𝑡,𝑡+ℎ𝛼(𝑠−𝑡)^𝛼−(𝑠−𝑡−ℎ)^𝛼=_𝑡^𝑡+ℎ(𝑠−𝜏)^𝛼−1𝑑𝜏≤_𝑡^𝑡+ℎ𝐶₀^𝛼−1ℎ^𝛼−1𝑑𝜏=𝐶₀^𝛼−1ℎ^𝛼. (32)
Thus, 1𝛼(𝑠−𝑡)^𝛼−(𝑠−𝑡−ℎ)^𝛼=𝑜(1). (33) By means of (30), Lemma 4 is established obviously.

Theorem 5 (the local truncation error). If the analytical solution of the problem (5) 𝑦(𝑡)∈𝐶⁴[0,𝑇], then the local truncation error of DSCM is 𝑂(ℎ⁴).

Proof. Applying the Taylor expansion and the definition of 𝑆(𝑡,ℎ), we have 𝑦(𝑡)=1+2𝑡−𝑡_𝑖−1ℎ𝑡_𝑖−𝑡ℎ²𝑦𝑡_𝑖−1+𝑡−𝑡_𝑖−1𝑡_𝑖−𝑡²ℎ³ℎ𝑦^′𝑡_𝑖−1+𝑡1+2^𝑖−𝑡ℎ𝑡−𝑡_𝑖−1ℎ²𝑦𝑡_𝑖−𝑡−𝑡_𝑖−1²𝑡_𝑖−𝑡ℎ³ℎ𝑦^′𝑡_𝑖+𝑂𝑡−𝑡_𝑖−1⁴=𝑦(𝑡)+𝑂𝑡−𝑡_𝑖−1⁴,𝑦^′(𝑡)=𝑦^′(𝑡)+𝑂𝑡−𝑡_𝑖−1³,∀𝑡∈𝑡_𝑖−1,𝑡_𝑖. (34) From the definition of the Caputo fractional derivative, we obtain 1Γ(1−𝛼)^{𝑡𝑖−1+𝜙𝑗ℎ}_𝑡𝑖−1𝑦^′(𝜏)+𝑂𝜏−𝑡_𝑖−1³𝑑𝜏𝑡_𝑖−1+𝜙_𝑗ℎ−𝜏^𝛼+^𝑖−1_𝑛=1^𝑡𝑛_𝑡𝑛−1𝑦^′(𝜏)𝑑𝜏𝑡_𝑖−1+𝜙_𝑗ℎ−𝜏^𝛼=𝑓𝑡_𝑖−1+𝜙_𝑗ℎ,𝑦𝑡_𝑖−1+𝜙_𝑗ℎ+𝑂𝜙_𝑗ℎ⁴,𝑗=1,2. (35) That is, 1Γ(1−𝛼)^{𝑡𝑖−1+𝜙𝑗ℎ}_𝑡𝑖−1𝑦^′(𝜏)𝑑𝜏𝑡_𝑖−1+𝜙_𝑗ℎ−𝜏^𝛼=𝑓𝑡_𝑖−1+𝜙_𝑗ℎ,𝑦𝑡_𝑖−1+𝜙_𝑗ℎ+𝑂𝜙_𝑗ℎ⁴−1Γ(1−𝛼)^𝑖−1_𝑛=1^𝑡𝑛_𝑡𝑛−1𝑦^′(𝜏)𝑑𝜏𝑡_𝑖−1+𝜙_𝑗ℎ−𝜏^𝛼+^{𝑡𝑖−1+𝜙𝑗ℎ}_𝑡𝑖−1𝑂𝜏−𝑡_𝑖−1³𝑑𝜏𝑡_𝑖−1+𝜙_𝑗ℎ−𝜏^𝛼, (36) where 𝑗=1,2. Moreover, 𝑦𝑡_𝑖ℎ𝑦^′𝑡_𝑖𝑦𝑡=𝐴^𝑖−1ℎ𝑦^′𝑡_𝑖−1+ℎ^𝛼Γ()𝐵×𝑓𝑡4−𝛼_𝑖−1+𝜙₁ℎ,𝑦𝑡_𝑖−1+𝜙₁ℎ+𝑂𝜙₁ℎ⁴𝑓𝑡_𝑖−1+𝜙₂ℎ,𝑦𝑡_𝑖−1+𝜙₂ℎ+𝑂𝜙₂ℎ⁴−ℎ^𝛼𝐴^{−1𝑖−1𝑛=1𝑡𝑛𝑡𝑛−1}𝑦^′(𝜏)𝑑𝜏𝑡_𝑖−1+𝜙₁ℎ−𝜏^𝛼_{𝑖−1𝑛=1𝑡𝑛𝑡𝑛−1}𝑦^′(𝜏)𝑑𝜏𝑡_𝑖−1+𝜙₂ℎ−𝜏^𝛼+1Γ(1−𝛼)^{𝑡𝑖−1+𝜙1ℎ}_𝑡𝑖−1𝑂𝜏−𝑡_𝑖−1³𝑑𝜏𝑡_𝑖−1+𝜙₁ℎ−𝜏^𝛼1Γ(1−𝛼)^{𝑡𝑖−1+𝜙2ℎ}_𝑡𝑖−1𝑂𝜏−𝑡_𝑖−1³𝑑𝜏𝑡_𝑖−1+𝜙₂ℎ−𝜏^𝛼. (37)
Supposing 𝑆(𝑡,ℎ)=𝑦(𝑡), for all 𝑡∈[0,𝑡_𝑖−1], and substituting 𝑦(𝑡) into (13), we have 𝑦𝑡_𝑖ℎ𝑦^′𝑡_𝑖𝑦𝑡=𝐴^𝑖−1ℎ𝑦^′𝑡_𝑖−1+ℎ^𝛼Γ()𝐵×𝑓𝑡4−𝛼_𝑖−1+𝜙₁ℎ,𝑆𝑡_𝑖−1+𝜙₁𝑓𝑡ℎ,ℎ_𝑖−1+𝜙₂ℎ,𝑆𝑡_𝑖−1+𝜙₂+ℎ,ℎ𝐻(𝑖)+𝑇_𝑖,𝐻(𝑖)=−ℎ^𝛼𝐴^{−1𝑖−1𝑛=1𝑡𝑛𝑡𝑛−1}𝑆^′(𝜏,ℎ)𝑑𝜏𝑡_𝑖−1+𝜙₁ℎ−𝜏^𝛼_{𝑖−1𝑛=1𝑡𝑛𝑡𝑛−1}𝑆^′(𝜏,ℎ)𝑑𝜏𝑡_𝑖−1+𝜙₂ℎ−𝜏^𝛼=−ℎ^𝛼𝐴^{−1𝑖−1𝑛=1𝑡𝑛𝑡𝑛−1}𝑦^′(𝜏)𝑑𝜏𝑡_𝑖−1+𝜙₁ℎ−𝜏^𝛼_{𝑖−1𝑛=1𝑡𝑛𝑡𝑛−1}𝑦^′(𝜏)𝑑𝜏𝑡_𝑖−1+𝜙₂ℎ−𝜏^𝛼,𝑆_𝑖(𝑡,ℎ)=1+2𝑡−𝑡_𝑖−1ℎ𝑡_𝑖−𝑡ℎ²𝑦𝑡_𝑖−1+𝑡−𝑡_𝑖−1𝑡_𝑖−𝑡²ℎ³ℎ𝑦^′𝑡_𝑖−1+𝑡1+2^𝑖−𝑡ℎ𝑡−𝑡_𝑖−1ℎ²𝑦𝑡_𝑖−𝑡−𝑡_𝑖−1²𝑡_𝑖−𝑡ℎ³ℎ𝑦^′𝑡_𝑖, (38)(39) where 𝑇_𝑖 is the local truncation error.
Combining (38), (37) with the Lipschitz condition yields 𝑇_𝑖≤ℎ^𝛼(𝐵‖𝑂𝜙𝐿Γ4−𝛼)‖¹ℎ⁴𝑂𝜙₂ℎ⁴+ℎ^𝛼𝐴⁻¹1Γ()𝜙2−𝛼¹ℎ^1−𝛼𝑂ℎ³1Γ()𝜙2−𝛼²ℎ^1−𝛼𝑂ℎ³,𝑇_𝑖≤𝐶ℎ⁴, (40) where 𝐶 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) 𝑦(𝑡)∈𝐶⁴[0,𝑇] and the matrixes 𝐶(𝜙₁,𝜙₂) and 𝐷(𝜙₁,𝜙₂) in (14) satisfy that 𝐷⁻¹(𝜙₁,𝜙₂) exists, and 𝐷⁻¹𝜙₁,𝜙₂𝐶𝜙₁,𝜙₂<1, (41) then the numerical solutions of ISCM for problem (5) 𝑆(𝑡,ℎ) satisfy that lim_ℎ→0𝑆(𝑡,ℎ)=𝑦(𝑡).

Proof. For (23), we have 𝑆(𝑡,ℎ)_{𝑡=𝑡𝑖+𝜙𝑗ℎ}=𝑦₀+1Γ(𝛼)^{𝑡𝑖+𝜙𝑗ℎ}₀𝑡_𝑖+𝜙_𝑗ℎ−𝜏^𝛼−1𝑓𝜏,𝑆(𝜏,ℎ)𝑑𝜏,𝑗=1,2,𝑆(0,ℎ)=𝑦₀,𝑆𝑡_𝑖+𝜙_𝑗ℎ,ℎ=𝑦₀+1Γ(𝛼)^{𝑡𝑖+𝜙𝑗ℎ}_𝑡𝑖𝑡_𝑖+𝜙_𝑗ℎ−𝜏^𝛼−1𝑓𝜏,𝑆()+𝜏,ℎ𝑑𝜏_{𝑖𝑛=1𝑡𝑛𝑡𝑛−1}𝑡_𝑖+𝜙_𝑗ℎ−𝜏^𝛼−1𝑓𝜏,𝑆(𝜏,ℎ)𝑑𝜏𝑗=1,2. (42) It follows from (7) that the analytical solution 𝑦(𝑡) of the problem (5) satisfies 𝑦𝑡_𝑖+𝜙_𝑗ℎ(0)+1=𝑦Γ(𝛼)^{𝑡𝑖+𝜙𝑗ℎ}_𝑡𝑖𝑡_𝑖+𝜙_𝑗ℎ−𝜏^𝛼−1𝑓(𝜏)+𝜏,𝑦𝑑𝜏_{𝑖𝑛=1𝑡𝑛𝑡𝑛−1}𝑡_𝑖+𝜙_𝑗ℎ−𝜏^𝛼−1𝑓𝜏,𝑦(𝜏)𝑑𝜏,𝑗=1,2. (43)
Moreover, (44) and (45) yield 𝑆𝑡_𝑖+𝜙_𝑗ℎ,ℎ−𝑦𝑡_𝑖+𝜙_𝑗ℎ=1Γ(𝛼)^{𝑡𝑖+𝜙𝑗ℎ}_𝑡𝑖𝑡_𝑖+𝜙_𝑗ℎ−𝜏^𝛼−1×𝑓𝜏,𝑆()(𝜏)+𝜏,ℎ−𝑓𝜏,𝑦𝑑𝜏_{𝑖𝑛=1𝑡𝑛𝑡𝑛−1}𝑡_𝑖+𝜙_𝑗ℎ−𝜏^𝛼−1×𝑓𝜏,𝑆((𝜏=ℎ𝜏,ℎ)−𝑓𝜏,𝑦)𝑑𝜏_𝛼Γ(𝛼)×^𝜙𝑗₀𝜙_𝑗−𝜉^𝛼−1𝑓𝑡_𝑖+𝜉ℎ,𝑆𝑡_𝑖𝑡+𝜉ℎ,ℎ−𝑓_𝑖+𝜉ℎ,𝑦𝑡_𝑖++𝜉ℎ𝑑𝜉_{𝑖−1𝑛=010}𝑖+𝜙_𝑗−𝑛−𝜉^𝛼−1×𝑓𝑡_𝑛+𝜉ℎ,𝑆𝑡_𝑛𝑡+𝜉ℎ,ℎ−𝑓_𝑛+𝜉ℎ,𝑦𝑡_𝑛+𝜉ℎ𝑑𝜉,𝑗=1,2. (44)
By using the Lipschitz condition, we have 𝑆𝑡_𝑖+𝜙_𝑗ℎ,ℎ−𝑦𝑡_𝑖+𝜙_𝑗ℎ=1Γ(𝛼)^{𝑡𝑖+𝜙𝑗ℎ}_𝑡𝑖𝑡_𝑖+𝜙_𝑗ℎ−𝜏^𝛼−1×𝑓𝜏,𝑆()(𝜏)+𝜏,ℎ−𝑓𝜏,𝑦𝑑𝜏_{𝑖𝑛=1𝑡𝑛𝑡𝑛−1}𝑡_𝑖+𝜙_𝑗ℎ−𝜏^𝛼−1×𝑓𝜏,𝑆()(𝜏)≤1𝜏,ℎ−𝑓𝜏,𝑦𝑑𝜏Γ(𝛼)^{𝑡𝑖+𝜙𝑗ℎ}_𝑡𝑖𝑡_𝑖+𝜙_𝑗ℎ−𝜏^𝛼−1×𝑓𝜏,𝑆()(𝜏)+𝜏,ℎ−𝑓𝜏,𝑦𝑑𝜏_{𝑖𝑛=1𝑡𝑛𝑡𝑛−1}𝑡_𝑖+𝜙_𝑗ℎ−𝜏^𝛼−1×𝑓𝜏,𝑆()(𝜏)≤𝐿𝜏,ℎ−𝑓𝜏,𝑦𝑑𝜏Γ(𝛼)^{𝑡𝑖+𝜙𝑗ℎ}_𝑡𝑖𝑡_𝑖+𝜙_𝑗ℎ−𝜏^𝛼−1𝑆()(𝜏+𝜏,ℎ−𝑦)𝑑𝜏_{𝑖𝑛=1𝑡𝑛𝑡𝑛−1}𝑡_𝑖+𝜙_𝑗ℎ−𝜏^𝛼−1𝑆()(𝜏=ℎ𝜏,ℎ−𝑦)𝑑𝜏_𝛼𝐿Γ(𝛼)^𝜙𝑗₀𝜙_𝑗−𝜏^𝛼−1𝑆𝑡_𝑖+𝜏ℎ,ℎ−𝑦𝑡_𝑖++𝜏ℎ𝑑𝜏_{𝑖−1𝑛=010}𝑖−𝑛+𝜙_𝑗−𝜏^𝛼−1×𝑆𝑡_𝑛+𝜏ℎ,ℎ−𝑦𝑡_𝑛+𝜏ℎ𝑑𝜏,𝑗=1,2. (45)
It follows from the definition of 𝑆(𝑡,ℎ) that 𝑦𝑡_𝑖𝑦𝑡+ℎ𝑢_𝑖=+ℎ𝑣𝑢²(2𝑢+1)𝑢²𝑢𝑣²(2𝑣+1)𝑣²𝑣𝑦𝑡_𝑖ℎ𝑦^′𝑡_𝑖+𝑢²(2𝑢+1)−𝑢²𝑢𝑣²(2𝑣+1)−𝑣²𝑣𝑦𝑡_𝑖+1ℎ𝑦^′𝑡_𝑖+1()𝑦=𝐶𝑢,𝑣^𝑖ℎ𝑦^′_𝑖()𝑦+𝐷𝑢,𝑣^𝑖+1ℎ𝑦^′_𝑖+1+𝑂ℎ⁴,∀𝑢,𝑣∈[0,1], (46) where 𝑢=1−𝑢, 𝑣=1−𝑣. Obviously, for (14) and (46), the matrixes 𝐶(𝑢,𝑣) and 𝐷(𝑢,𝑣) are continuous on [0,1]×[0,1]; thus, we have 𝑆𝑡_𝑖+ℎ𝑢,ℎ−𝑦𝑡_𝑖+ℎ𝑢𝑆𝑡_𝑖+ℎ𝑣,ℎ−𝑦𝑡_𝑖=+ℎ𝑣𝑢²(2𝑢+1)𝑢²𝑣𝑣²(2𝑣+1)𝑣²𝑣𝑆_𝑖⁽⁰⁾−𝑦𝑡_𝑖𝑆_𝑖⁽¹⁾−ℎ𝑦^′𝑡_𝑖+𝑢²(2𝑢+1)−𝑢²𝑢𝑣²(2𝑣+1)−𝑣²𝑣𝑆⁽⁰⁾_𝑖+1−𝑦𝑡_𝑖+1𝑆⁽¹⁾_𝑖+1−ℎ𝑦^′𝑡_𝑖+1=𝐶(𝑢,𝑣)𝑆_𝑖⁽⁰⁾−𝑦𝑡_𝑖𝑆_𝑖⁽¹⁾−ℎ𝑦^′𝑡_𝑖+𝐷(𝑢,𝑣)𝑆⁽⁰⁾_𝑖+1−𝑦𝑡_𝑖+1𝑆⁽¹⁾_𝑖+1−ℎ𝑦^′𝑡_𝑖+1+𝑂ℎ⁴,∀𝑢,𝑣∈[0,1]. (47) Therefore, there exist 𝜆₁>0, 𝜆₂>0 such that 𝜆₁=sup_{0≤𝑢<𝑣≤1}‖𝐶(𝑢,𝑣)‖,𝜆₂=sup_{0≤𝑢<𝑣≤1}‖𝐷(𝑢,𝑣)‖,𝜆=𝜆₁+𝜆₂. (48)
Obviously, ^𝜙𝑗₀𝜙_𝑗−𝜉^𝛼−1𝑑𝜉=−1𝛼𝜙_𝑗−𝜉^𝛼𝜙𝑗₀=1𝛼𝜙_𝑗^𝛼≤1𝛼,∀𝜙_𝑗∈𝐶₀,1. (49) Applying the integral mean theorem, we have ¹₀𝑖−𝑛+𝜙_𝑗−𝜉^𝛼−1=𝑑𝜉𝑖−𝑛+𝜙_𝑗−𝜃^𝛼−1()≤()0≤𝜃≤1𝑖−𝑛^𝛼−1𝜙1+^𝑗−𝜃𝑖−𝑛^𝛼−1≤(𝑖−𝑛)^𝛼−1𝜙1+^𝑗−1𝑖−𝑛^𝛼−1. (50) It follows from the inequalities 0<𝐶₀≤𝜙₁<𝜙₂≤1 and 0≤𝑛≤𝑖−1 that ¹₀𝑖−𝑛+𝜙_𝑗−𝜉^𝛼−1𝑑𝜉≤(𝑖−𝑛)^𝛼−1𝜙1+^𝑗−1𝑖−𝑛^𝛼−1≤(𝑖−𝑛)^𝛼−11+𝜙_𝑗−1^𝛼−1≤(𝑖−𝑛)^𝛼−1𝐶₀^𝛼−1,𝑗=1,2. (51)
Combining (45), (48), and (49) with (51) and defining 𝜖_𝑖=𝑆_𝑖⁽⁰⁾−𝑦𝑡_𝑖𝑆_𝑖⁽¹⁾−ℎ𝑦^′𝑡_𝑖,𝜖^′_𝑖=𝑆𝑡_𝑖+𝜙1,ℎ−𝑦𝑡_𝑖+𝜙1𝑆𝑡_𝑖+𝜙2,ℎ−𝑦𝑡_𝑖+𝜙2,𝑋_𝑖=max_{0≤𝑛≤𝑖+1}𝜖_𝑛, (52) where 𝑡_{𝑖+𝜙𝑗}=𝑡_𝑖+𝜙_𝑗ℎ, 𝑗=1,2, we can obtain 𝜖^′_𝑖≤ℎ^𝛼𝐿1Γ(𝛼)𝜆₁𝜖_𝑖+𝜆₂𝜖_𝑖+1𝛼+𝐶₁ℎ⁴Γ())ℎ)𝛼+1((𝑖+1^𝛼+1Γ(𝛼)^𝑖−1_𝑛=0𝐶₀^𝛼−1(𝑖−𝑛)^𝛼−1×𝜆₁𝜖_𝑛+𝜆₂𝜖_𝑛+1≤ℎ^𝛼𝐿1Γ(𝛼+1)𝜆𝑋_𝑖+1Γ(𝛼)^𝑖−1_𝑛=0𝐶₀^𝛼−1×(𝑖−𝑛)^𝛼−1𝜆𝑋_𝑛+𝐶₂ℎ⁴, (53) where 𝐶₁,𝐶₂ are some positive constants. If the matrix 𝐷(𝜙₁,𝜙₂) is invertible, then it follows from (47) that we have 𝑆⁽⁰⁾_𝑖+1−𝑦𝑡_𝑖+1𝑆⁽¹⁾_𝑖+1−ℎ𝑦^′𝑡_𝑖+1=−𝐷⁻¹𝜙₁,𝜙₂𝐶𝜙₁,𝜙₂𝑆_𝑖⁽⁰⁾−𝑦𝑡_𝑖𝑆_𝑖⁽¹⁾−ℎ𝑦^′𝑡_𝑖+𝐷⁻¹𝜙₁,𝜙₂𝑆𝑡_𝑖+ℎ𝜙₁,ℎ−𝑦𝑡_𝑖+ℎ𝜙₁𝑆𝑡_𝑖+ℎ𝜙₂,ℎ−𝑦𝑡_𝑖+ℎ𝜙₂+𝑂ℎ⁴. (54) Hence, 𝜖_𝑖+1≤𝐷⁻¹𝜙₁,𝜙₂𝐶𝜙₁,𝜙₂𝜖_𝑖+𝐷⁻¹𝜙₁,𝜙₂𝜖^′_𝑖+𝐶₃ℎ⁴≤𝐷⁻¹𝜙₁,𝜙₂𝐶𝜙₁,𝜙₂𝑋_𝑖+𝐷⁻¹𝜙₁,𝜙₂𝜖^′_𝑖+𝐶₃ℎ⁴, (55) where 𝐶₃ is an appropriate positive constant. Take ‖𝜖_ℓ+1‖=max_{0≤𝑛≤𝑖+1}‖𝜖_𝑛‖=𝑋_𝑖, 0≤ℓ≤𝑖. Then 𝑋_ℓ=max_{0≤𝑛≤ℓ+1}‖𝜖_𝑛‖=‖𝜖_ℓ+1‖, 𝑋_𝑖=𝑋_ℓ, and 𝑋_ℓ=𝜖_ℓ+1≤𝐷⁻¹𝜙₁,𝜙₂𝐶𝜙₁,𝜙₂𝜖_ℓ+𝐷⁻¹𝜙₁,𝜙₂𝜖^′_ℓ+𝐶₃ℎ⁴≤𝐷⁻¹𝜙₁,𝜙₂𝐶𝜙₁,𝜙₂𝑋_ℓ+𝐷⁻¹𝜙₁,𝜙₂𝜖^′_ℓ+𝐶₃ℎ⁴. (56)
Since ‖𝐷⁻¹(𝜙₁,𝜙₂)𝐶(𝜙₁,𝜙₂)‖<1, we have 1−𝐷⁻¹𝜙₁,𝜙₂𝐶𝜙₁,𝜙₂𝑋_ℓ≤𝐷⁻¹𝜙₁,𝜙₂𝜖^′_ℓ+𝐶₃ℎ⁴,𝑋_ℓ≤𝐶₁𝜖^′_ℓ+𝐶₄ℎ⁴, (57) where 𝐶₁=‖𝐷⁻¹(𝜙₁,𝜙₂)‖/(1−‖𝐷⁻¹(𝜙₁,𝜙₂)𝐶(𝜙₁,𝜙₂)‖), 𝐶₄ is an appropriate positive constant. According to (53), we obtain 𝑋_ℓ≤𝐶₁ℎ^𝛼𝐿1Γ(𝛼+1)𝜆𝑋_ℓ+1Γ(𝛼)×^ℓ−1_𝑛=0𝐶₀^𝛼−1(ℓ−𝑛)^𝛼−1𝜆𝑋_𝑛+𝐶₅ℎ⁴,𝐶1−¹ℎ^𝛼𝐿𝜆Γ()𝑋𝛼+1^ℓ≤𝐶₁ℎ^𝛼𝐿1Γ(𝛼)^ℓ−1_𝑛=0𝐶₀^𝛼−1(ℓ−𝑛)^𝛼−1𝜆𝑋_𝑛+𝐶₅ℎ⁴, (58) where 𝐶₅>0 is an appropriate constant. Obviously, there exist 𝐶∈(0,1) and ℎ₁>0 such that 𝐶₁ℎ^𝛼𝐿𝜆/Γ(𝛼+1)≤𝐶<1 for ℎ<ℎ₁. Thus, 𝑋_ℓ≤𝐶₁ℎ^𝛼𝐿𝜆𝐶11−Γ(𝛼)^ℓ−1_𝑛=0𝐶₀^𝛼−1(ℓ−𝑛)^𝛼−1𝑋_𝑛+𝐶₆ℎ⁴. (59) Let 𝐶=𝐶₁𝐶₀^𝛼−1𝜆/(1−𝐶)Γ(𝛼); then 𝜖_𝑖+1≤𝑋_ℓ≤𝐶ℎ^𝛼ℓ−1_𝑛=0(ℓ−𝑛)^𝛼−1𝑋_𝑛+𝐶₆ℎ⁴, (60) where 𝐶₆>0 is an appropriate constant. Applying Lemma 3, we get 𝜖_𝑖+1≤𝑋_ℓ≤𝐶₆ℎ⁴𝐸_𝛼𝐶Γ(𝛼)(ℓℎ)^𝛼≤𝐶₆ℎ⁴𝐸_𝛼𝐶Γ(𝛼)𝑇^𝛼. (61) By using the convergence of the Mittag-Leffler function 𝐸_𝛼(𝑧) [6], we have 𝑋_ℓ⟶0,𝜖_𝑖+1⟶0(ℎ⟶0). (62)
By means of (24), (47), we obtain sup_{0≤𝑡≤𝑇}𝑆(𝑡,ℎ)−𝑦(𝑡)≤max_{0≤𝑖≤𝑁−1}sup_{0≤𝑢<𝑣≤1}‖𝐶(𝑢,𝑣)‖𝜖_𝑖+sup_{0≤𝑢<𝑣≤1}‖𝐷(𝑢,𝑣)‖𝜖_𝑖+1+𝐶₇ℎ⁴≤𝜆₁+𝜆₂𝐸_𝛼𝐶Γ(𝛼)𝑇^𝛼𝐶₆ℎ⁴+𝐶₇ℎ⁴≤𝐶₇+𝜆₁+𝜆₂𝐶₆𝐸_𝛼𝐶Γ(𝛼)𝑇^𝛼ℎ⁴, (63) where 𝐶₇ is a positive constant. Consequently, we have lim_ℎ→0𝑆(𝑡,ℎ)=𝑦(𝑡). (64)

Remark 7. (1) If 𝜙₁≠𝜙₂, then 𝐷𝜙₁,𝜙₂=𝜙²₁2𝜙₁+1−𝜙²₁𝜙₁𝜙²₂2𝜙₂+1−𝜙²₂𝜙₂ (65) is invertible, and𝐷⁻¹𝜙₁,𝜙₂=−−1+𝜙₂𝜙²₁−𝜙₂+𝜙₁−1+𝜙₁𝜙²₂−𝜙₂+𝜙₁−−3+2𝜙₂𝜙²₁−𝜙₂+𝜙₁−3+2𝜙₁𝜙²₂−𝜙₂+𝜙₁. (66)
(2) Consider 𝐷⁻¹𝜙₁,𝜙₂𝐶𝜙₁,𝜙₂=−𝜙₁𝜙₂+𝜙₁+𝜙₂−1+𝜙₂−1+𝜙₁𝜙²₁𝜙²₂−−1+𝜙₂−1+𝜙₁𝜙₁𝜙₂2𝜙²₁−3𝜙₁+2𝜙₁𝜙₂−3𝜙₂+2𝜙²₂𝜙²₁𝜙²₂−𝜙₁𝜙₂−2𝜙₁+3−2𝜙₂𝜙₁𝜙₂. (67)

When ‖𝐷⁻¹(𝜙₁,𝜙₂)𝐶(𝜙₁,𝜙₂)‖<1, the range of values 𝜙₁,𝜙₂ can be obtained by using Algorithm 8, which can be shown in Figures 1 and 2, respectively.

(a)

(b)

Figure 1 

The range of values 𝜙1,𝜙2 when ‖𝐷−1𝐶‖1<1 (a) or ‖𝐷−1𝐶‖2<1 (b).

Figure 2 

The range of values 𝜙1,𝜙2 when ‖𝐷−1𝐶‖∞<1.

Algorithm 8.  For 0≤𝜙₁≤1 For 𝜙₁<𝜙₂≤1     if  ‖𝐷⁻¹(𝜙₁,𝜙₂)𝐶(𝜙₁,𝜙₂)‖<1       record 𝜙₁,𝜙₂;      end   end end plot(𝜙₁,𝜙₂).

Lemma 9. If the analytical solution of the problem (5) 𝑦(𝑡)∈𝐶⁴[0,𝑇], the spline functions 𝑆(𝑡,ℎ) and 𝑆(𝑡,ℎ) 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)If lim_ℎ→0𝑆(𝑡,ℎ)=𝑦(𝑡), then lim_ℎ→0𝑆(𝑡,ℎ)=𝑦(𝑡), where 𝑦(𝑡) is the analytical solution of the problem (5).(ii)Assume that 𝑦𝑆(𝑡,ℎ),𝑧𝑆(𝑡,ℎ) and 𝑦𝑆(𝑡,ℎ),𝑧𝑆(𝑡,ℎ) are the numerical solutions obtained by DSCM and ISCM for the problem ^𝐶₀𝐷^𝛼_𝑡𝑦(𝑡)(𝑡)[]𝑦(0)=𝑓𝑡,𝑦,𝑡∈𝐼=0,𝑇,0<𝛼≤1,=𝑦₀, (68)and the perturbed problem ^𝐶₀𝐷^𝛼_𝑡𝑧(𝑡)((𝑡[]𝑧(0)=𝑓𝑡,𝑧)),𝑡∈𝐼=0,𝑇,0<𝛼≤1,=𝑧₀, (69) where 𝑓 satisfies the Lipschitz condition (6), respectively. Then there exists the constant ℎ₀>0 such that, for all ℎ<ℎ₀, when |𝑦𝑆(𝑡_𝑛,ℎ)−𝑧𝑆(𝑡_𝑛,ℎ)|≤𝐶₁|𝑦₀−𝑧₀|, there exists the constant 𝐶>0 such that |𝑦𝑆(𝑡_𝑛,ℎ)−𝑧𝑆(𝑡_𝑛,ℎ)|≤𝐶|𝑦₀−𝑧₀|; here 𝐶₁>0 is a constant, 𝑛=1,2,…,𝑁.

Proof. Firstly, we give the proof of the conclusion (i). Applying DSCM, we have ^𝐶₀𝐷^𝛼_𝑡𝑆(𝑡,ℎ)_{𝑡=𝑡𝑖−1+𝜙ℎ}=𝑓𝑡_𝑖−1+𝜙ℎ,𝑆𝑡_𝑖−1+𝜙ℎ,ℎ,𝑖=1,2…,𝑁. (70) Since 𝑆(𝑡,ℎ)∈𝐶¹[0,𝑇], we obtain 𝑆𝑡_𝑖−1+𝜙ℎ,ℎ=𝑦₀+1Γ(𝛼)^𝑡𝑖−1₀^+𝜙ℎ𝑡_𝑖−1+𝜙ℎ−𝜏^𝐶₀^𝛼−1𝐷^𝛼_𝜏𝑆(𝜏,ℎ)𝑑𝜏=𝑦₀+1Γ(𝛼)^𝑖−1_𝑛=1^𝑡𝑛_𝑡𝑛−1𝑡_𝑖−1+𝜙ℎ−𝜏^𝐶₀^𝛼−1𝐷^𝛼_𝜏𝑆(𝜏,ℎ)𝑑𝜏+^𝑡𝑖−1_𝑡^+𝜙ℎ_𝑖−1𝑡_𝑖−1+𝜙ℎ−𝜏^𝐶₀^𝛼−1𝐷^𝛼_𝜏𝑆(𝜏,ℎ)𝑑𝜏. (71) It follows from Lemma 4 that 𝑆𝑡_𝑖−1+𝜙ℎ,ℎ=𝑦₀+1Γ(𝛼)^𝑖−1_𝑛=1^𝐶0𝐷^𝛼_𝜏𝑆(𝜏,ℎ)_{𝜏=𝑡𝑛−1+𝜙}×^𝑡𝑛_𝑡𝑛−1𝑡_𝑖−1+𝜙ℎ−𝜏^𝛼−1(ℎ)+𝑑𝜏+𝑜_𝐶0𝐷^𝛼_𝜏𝑆(𝜏,ℎ)_{𝜏=𝑡𝑖−1+𝜙ℎ}×^𝑡𝑖−1_𝑡^+𝜙ℎ_𝑖−1𝑡_𝑖−1+𝜙ℎ−𝜏^𝛼−1𝑑𝜏+𝑜(ℎ)=𝑦₀+1Γ(𝛼)^{𝑖−1𝐶0}_𝑛=1𝐷^𝛼_𝜏𝑆(𝜏,ℎ)_{𝜏=𝑡𝑛−1+𝜙}×^𝑡𝑛_𝑡𝑛−1𝑡_𝑖−1+𝜙ℎ−𝜏^𝛼−1+𝑑𝜏_𝐶0𝐷^𝛼_𝜏𝑆(𝜏,ℎ)_{𝜏=𝑡𝑖−1+𝜙ℎ}×^𝑡𝑖−1_𝑡^+𝜙ℎ_𝑖−1𝑡_𝑖−1+𝜙ℎ−𝜏^𝛼−1𝑑𝜏+𝑜(1), (72) where 𝑜(1)→0  as  ℎ→0.
By means of (70) and (72), we obtain 𝑆𝑡_𝑖−1+𝜙ℎ,ℎ=𝑦₀+1Γ(𝛼)^𝑖−1_𝑛=1𝑓𝑡_{𝑛−1+𝜙},𝑆𝑡_{𝑛−1+𝜙}×,ℎ_{𝑡𝑛𝑡𝑛−1}𝑡_𝑖−1+𝜙ℎ−𝜏^𝛼−1𝑡𝑑𝜏+𝑓_𝑖−1+𝜙ℎ,𝑆𝑡_𝑖−1+𝜙ℎ,ℎ×^𝑡𝑖−1_𝑡^+𝜙ℎ_𝑖−1𝑡_𝑖−1+𝜙ℎ−𝜏^𝛼−1𝑑𝜏+𝑜(1). (73) Applying Lemma 4, we get 𝑆𝑡_𝑖−1+𝜙ℎ,ℎ=𝑦₀+1Γ(𝛼)^𝑖−1_𝑛=1^𝑡𝑛_𝑡𝑛−1𝑡_𝑖−1+𝜙ℎ−𝜏^𝛼−1(((ℎ)+×𝑓𝜏,𝑆𝜏,ℎ))𝑑𝜏+𝑜_{𝑡𝑖−1𝑡+𝜙ℎ𝑖−1}𝑡_𝑖−1+𝜙ℎ−𝜏^𝛼−1×𝑓(𝜏,𝑆(𝜏,ℎ))𝑑𝜏+𝑜(ℎ)+𝑜(1)=𝑦₀+1Γ(𝛼)^𝑖−1_𝑛=1^𝑡𝑛_𝑡𝑛−1𝑡_𝑖−1+𝜙ℎ−𝜏^𝛼−1𝑓((+𝜏,𝑆𝜏,ℎ))𝑑𝜏_{𝑡𝑖−1𝑡+𝜙ℎ𝑖−1}𝑡_𝑖−1+𝜙ℎ−𝜏^𝛼−1×𝑓(𝜏,𝑆(𝜏,ℎ))𝑑𝜏+𝑜(1), (74) where there exists a constant 𝐶₈>0, 𝛾>0 such that |𝑜(1)|=𝐶₈ℎ^𝛾. Based on the definitions of 𝑆(𝑡,ℎ),𝑆(𝑡,ℎ), it follows that 𝑆𝑡_𝑖−1+𝜙_𝑗ℎ,ℎ−𝑆𝑡_𝑖−1+𝜙_𝑗≤1ℎ,ℎΓ(𝛼)^{𝑡𝑖−1+𝜙𝑗ℎ}_𝑡𝑖−1𝑡_𝑖−1+𝜙_𝑗ℎ−𝜏^𝛼−1×𝑓(𝜏,𝑆(𝜏,ℎ))−𝑓𝜏,𝑆()+𝜏,ℎ𝑑𝜏_{𝑖−1𝑛=1𝑡𝑛𝑡𝑛−1}𝑡_𝑖−1+𝜙_𝑗ℎ−𝜏^𝛼−1×𝑓(𝜏,𝑆(𝜏,ℎ))−𝑓𝜏,𝑆(+𝐶𝜏,ℎ)𝑑𝜏₈ℎ^𝛾≤1Γ(𝛼)^{𝑡𝑖−1+𝜙𝑗ℎ}_𝑡𝑖−1𝑡_𝑖−1+𝜙_𝑗ℎ−𝜏^𝛼−1×𝑓(𝜏,𝑆(𝜏,ℎ))−𝑓𝜏,𝑆()+𝜏,ℎ𝑑𝜏_{𝑖−1𝑛=1𝑡𝑛𝑡𝑛−1}𝑡_𝑖−1+𝜙_𝑗ℎ−𝜏^𝛼−1×𝑓(𝜏,𝑆(𝜏,ℎ))−𝑓𝜏,𝑆(𝜏,ℎ)𝑑𝜏+𝐶₈ℎ^𝛾≤𝐿Γ(𝛼)^{𝑡𝑖−1+𝜙𝑗ℎ}_𝑡𝑖−1𝑡_𝑖−1+𝜙_𝑗ℎ−𝜏^𝛼−1×𝑆(𝜏,ℎ)−𝑆(+𝜏,ℎ)𝑑𝜏_{𝑖−1𝑛=1𝑡𝑛𝑡𝑛−1}𝑡_𝑖−1+𝜙_𝑗ℎ−𝜏^𝛼−1×𝑆(𝜏,ℎ)−𝑆(𝜏,ℎ)𝑑𝜏+𝐶₈ℎ^𝛾,𝑗=1,2. (75)
Based on the similar proof to that of Lemma 6, it follows from Lemma 3 and the inequality ‖𝐷⁻¹(𝜙₁,𝜙₂)𝐶(𝜙₁,𝜙₂)‖<1 that 𝜖_𝑖≤𝐶₉ℎ^𝛾𝐸_𝛼𝐶Γ(𝛼)𝑇^𝛼,𝑖=1,2,…,𝑁, (76) where 𝜖_𝑖=𝑆_𝑖⁽⁰⁾−𝑆_𝑖⁽⁰⁾𝑆_𝑖⁽¹⁾−𝑆_𝑖⁽¹⁾, (77) and 𝐶 is defined in Lemma 6, 𝐶₉>0 is a constant. By means of the convergence of the Mittag-Leffler function 𝐸_𝛼(𝑧) [6], we get 𝜖_𝑖⟶0,ℎ⟶0;lim_ℎ→0𝑆(𝑡,ℎ)=lim_ℎ→0𝑆(𝑡,ℎ). (78) From Lemma 6, we have lim_ℎ→0𝑆(𝑡,ℎ)=𝑦(𝑡),lim_ℎ→0𝑆(𝑡,ℎ)=lim_ℎ→0𝑆(𝑡,ℎ)=𝑦(𝑡). (79)
Now, we give the proof of the conclusion (ii) of Lemma 6. According to the conclusion (i), we get 𝑦𝑆𝑡_𝑛,ℎ−𝑧𝑆𝑡_𝑛=𝑡,ℎ𝑦𝑆_𝑛,ℎ−𝑦𝑆𝑡_𝑛,ℎ+𝑦𝑆𝑡_𝑛,ℎ−𝑧𝑆𝑡_𝑛,ℎ+𝑧𝑆𝑡_𝑛,ℎ−𝑧𝑆𝑡_𝑛≤𝑡,ℎ𝑦𝑆_𝑛,ℎ−𝑦𝑆𝑡_𝑛+𝑦,ℎ𝑆𝑡_𝑛,ℎ−𝑧𝑆𝑡_𝑛+𝑧,ℎ𝑆𝑡_𝑛,ℎ−𝑧𝑆𝑡_𝑛,ℎ≤max_{0≤𝑖≤𝑁}𝑦𝐸𝑡_𝑖,ℎ+max_{0≤𝑖≤𝑁}𝑧𝐸𝑡_𝑖,ℎ+𝐶₁𝑦₀−𝑧₀,𝑛=1,2,…,𝑁, (80) where 𝑦𝐸(𝑡,ℎ)=𝑦𝑆(𝑡,ℎ)−𝑦𝑆(𝑡,ℎ), 𝑧𝐸(𝑡,ℎ)=𝑧𝑆(𝑡,ℎ)−𝑧𝑆(𝑡,ℎ).
According to the proof of the conclusion (i), there exists ℎ₀>0 such that, for all ℎ<ℎ₀, max_{0≤𝑖≤𝑁}𝑦𝐸𝑡_𝑖,ℎ+max_{0≤𝑖≤𝑁}𝑧𝐸𝑡_𝑖,ℎ≤𝑦₀−𝑧₀. (81) Then, 𝑦𝑆𝑡_𝑛,ℎ−𝑧𝑆𝑡_𝑛,ℎ≤1+𝐶₁𝑦₀−𝑧₀,𝑛=1,2,…,𝑁. (82) Obviously, take 𝐶=1+𝐶₁; hence, the conclusion (ii) is right.

Theorem 10 (convergence of DSCM). If the analytical solution of the problem (5) 𝑦(𝑡)∈𝐶⁴[0,𝑇] and the matrixes 𝐶(𝜙₁,𝜙₂) and 𝐷(𝜙₁,𝜙₂) of (14) satisfy that 𝐷⁻¹(𝜙₁,𝜙₂) exists, and 𝐷⁻¹𝜙₁,𝜙₂𝐶𝜙₁,𝜙₂<1, (83) then DSCM is convergent, and lim_ℎ→0𝑆(𝑡,ℎ)=𝑦(𝑡).

Proof. This result follows directly from Lemmas 6 and 9.

Now, we consider the stability of DSCM.

Theorem 11 (stability of DSCM). If the analytical solution of the problem (5) 𝑦(𝑡)∈𝐶⁴[0,𝑇] and the matrixes 𝐶(𝜙₁,𝜙₂) and 𝐷(𝜙₁,𝜙₂) in (14) satisfy that 𝐷⁻¹(𝜙₁,𝜙₂) exists and 𝐷⁻¹𝜙₁,𝜙₂𝐶𝜙₁,𝜙₂<1, (84) 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<ℎ<ℎ₁ (ℎ₁>0 is a constant), one has(i)𝜖_𝑖≤𝐶𝑦₀−𝑧₀,𝑖=1,2,…,𝑁, (85) where 𝜖_𝑖=𝑦𝑆_𝑖⁽⁰⁾−𝑧𝑆_𝑖⁽⁰⁾𝑦𝑆_𝑖⁽¹⁾−𝑧𝑆_𝑖⁽¹⁾, (86)𝐶>0 is a constant, {𝑦𝑆_𝑖⁽⁰⁾,𝑦𝑆_𝑖⁽¹⁾} are the numerical solutions of the problem (68), and {𝑧𝑆_𝑖⁽⁰⁾,𝑧𝑆_𝑖⁽¹⁾} are the numerical solutions of the perturbed problem (69), and(ii)sup_{0≤𝑡≤𝑇}𝑦𝑆(𝑡,ℎ)−𝑧𝑆(𝑡,ℎ)≤𝐶^′𝑦₀−𝑧₀, (87) where 𝐶^′>0 is a constant, 𝑦𝑆(𝑡,ℎ) is the numerical solution of the problem (68), and 𝑧𝑆(𝑡,ℎ) is the numerical solution of the problem (69).
Moreover, DSCM is stable.

Proof. According to Lemma 9, 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 𝑦𝑆(𝑡,ℎ)_{𝑡=𝑡𝑖+𝜙𝑗ℎ}=𝑦₀+1Γ(𝛼)^{𝑡𝑖+𝜙𝑗ℎ}₀𝑡_𝑖+𝜙_𝑗ℎ−𝜏^𝛼−1×𝑓𝜏,𝑦𝑆()𝑦𝜏,ℎ𝑑𝜏,𝑗=1,2,𝑆(0,ℎ)=𝑦₀,𝑧𝑆(𝑡,ℎ)_{𝑡=𝑡𝑖+𝜙𝑗ℎ}=𝑧₀+1Γ(𝛼)^{𝑡𝑖+𝜙𝑗ℎ}₀𝑡_𝑖+𝜙_𝑗ℎ−𝜏^𝛼−1×𝑓𝜏,𝑧𝑆(𝑧𝜏,ℎ)𝑑𝜏,𝑗=1,2,𝑆(0,ℎ)=𝑧₀. (88) The previous equations can be written as 𝑦𝑆𝑡_𝑖+𝜙_𝑗ℎ,ℎ=𝑦₀+1Γ(𝛼)^{𝑡𝑖+𝜙𝑗ℎ}_𝑡𝑖𝑡_𝑖+𝜙_𝑗ℎ−𝜏^𝛼−1𝑓𝜏,𝑦𝑆()+𝜏,ℎ𝑑𝜏_{𝑖𝑛=1𝑡𝑛𝑡𝑛−1}𝑡_𝑖+𝜙_𝑗ℎ−𝜏^𝛼−1𝑓𝜏,𝑦𝑆(𝜏,ℎ)𝑑𝜏,𝑗=1,2,𝑧𝑆𝑡_𝑖+𝜙_𝑗ℎ,ℎ=𝑧₀+1Γ(𝛼)^{𝑡𝑖+𝜙𝑗ℎ}_𝑡𝑖𝑡_𝑖+𝜙_𝑗ℎ−𝜏^𝛼−1𝑓𝜏,𝑧𝑆(+𝜏,ℎ)𝑑𝜏_{𝑖𝑛=1𝑡𝑛𝑡𝑛−1}𝑡_𝑖+𝜙_𝑗ℎ−𝜏^𝛼−1𝑓𝜏,𝑧𝑆(𝜏,ℎ)𝑑𝜏,𝑗=1,2. (89)(90) Thus, 𝑦𝑆𝑡_𝑖+𝜙_𝑗ℎ,ℎ−𝑧𝑆𝑡_𝑖+𝜙_𝑗=𝑦ℎ,ℎ₀−𝑧₀+1Γ(𝛼)×^{𝑡𝑖+𝜙𝑗ℎ}_𝑡𝑖𝑡_𝑖+𝜙_𝑗ℎ−𝜏^𝛼−1×𝑓𝜏,𝑦𝑆(𝜏,ℎ)−𝑓𝜏,𝑧𝑆(+𝜏,ℎ)𝑑𝜏_{𝑖𝑛=1𝑡𝑛𝑡𝑛−1}𝑡_𝑖+𝜙_𝑗ℎ−𝜏^𝛼−1×𝑓𝜏,𝑦𝑆(𝜏,ℎ)−𝑓𝜏,𝑧𝑆()=𝑦𝜏,ℎ𝑑𝜏₀−𝑧₀+ℎ^𝛼Γ(𝛼)×^𝜙𝑗₀𝜙_𝑗−𝜉^𝛼−1×𝑓𝑡_𝑖+𝜉ℎ,𝑦𝑆𝑡_𝑖𝑡+𝜉ℎ,ℎ−𝑓_𝑖+𝜉ℎ,𝑧𝑆𝑡_𝑖++𝜉ℎ,ℎ𝑑𝜉_{𝑖−1𝑛=010}𝑖+𝜙_𝑗−𝑛−𝜉^𝛼−1×𝑓𝑡_𝑛+𝜉ℎ,𝑦𝑆𝑡_𝑛𝑡+𝜉ℎ,ℎ−𝑓_𝑛+𝜉ℎ,𝑧𝑆𝑡_𝑛+𝜉ℎ,ℎ𝑑𝜉,𝑗=1,2. (91) Applying the Lipschitz condition, we get 𝑦𝑆𝑡_𝑖+𝜙_𝑗ℎ,ℎ−𝑧𝑆𝑡_𝑖+𝜙_𝑗≤𝑦ℎ,ℎ₀−𝑧₀+1Γ(𝛼)^{𝑡𝑖+𝜙𝑗ℎ}_𝑡𝑖𝑡_𝑖+𝜙_𝑗ℎ−𝜏^𝛼−1×𝑓𝜏,𝑦𝑆(𝜏,ℎ)−𝑓𝜏,𝑧𝑆()+𝜏,ℎ𝑑𝜏_{𝑖𝑛=1𝑡𝑛𝑡𝑛−1}𝑡_𝑖+𝜙_𝑗ℎ−𝜏^𝛼−1×𝑓𝜏,𝑦𝑆(𝜏,ℎ)−𝑓𝜏,𝑧𝑆()≤𝑦𝜏,ℎ𝑑𝜏₀−𝑧₀+1Γ(𝛼)^{𝑡𝑖+𝜙𝑗ℎ}_𝑡𝑖𝑡_𝑖+𝜙_𝑗ℎ−𝜏^𝛼−1×𝑓𝜏,𝑦𝑆(𝜏,ℎ)−𝑓𝜏,𝑧𝑆()+𝜏,ℎ𝑑𝜏_{𝑖𝑛=1𝑡𝑛𝑡𝑛−1}𝑡_𝑖+𝜙_𝑗ℎ−𝜏^𝛼−1×𝑓𝜏,𝑦𝑆(𝜏,ℎ)−𝑓𝜏,𝑧𝑆()≤𝑦𝜏,ℎ𝑑𝜏₀−𝑧₀+𝐿Γ(𝛼)^{𝑡𝑖+𝜙𝑗ℎ}_𝑡𝑖𝑡_𝑖+𝜙_𝑗ℎ−𝜏^𝛼−1×𝑦𝑆(𝜏,ℎ)−𝑧𝑆(+𝜏,ℎ)𝑑𝜏_{𝑖𝑛=1𝑡𝑛𝑡𝑛−1}𝑡_𝑖+𝜙_𝑗ℎ−𝜏^𝛼−1×𝑦𝑆(𝜏,ℎ)−𝑧𝑆(=𝑦𝜏,ℎ)𝑑𝜏₀−𝑧₀+ℎ^𝛼𝐿Γ(𝛼)^𝜙𝑗₀𝜙_𝑗−𝜏^𝛼−1×𝑦𝑆𝑡_𝑖+𝜏ℎ,ℎ−𝑧𝑆𝑡_𝑖++𝜏ℎ,ℎ𝑑𝜏_{𝑖−1𝑛=010}𝑖−𝑛+𝜙_𝑗−𝜏^𝛼−1×𝑦𝑆𝑡_𝑛+𝜏ℎ,ℎ−𝑧𝑆𝑡_𝑛+𝜏ℎ,ℎ𝑑𝜏, (92) where 𝑗=1,2. Note that the definitions of 𝑦𝑆(𝑡,ℎ),𝑧𝑆(𝑡,ℎ). Let 𝜆₁=sup_{0≤𝑢<𝑣≤1}‖𝐶(,𝜆𝑢,𝑣)‖₂=max_{0≤𝑢<𝑣≤1}‖𝐷(𝑢,𝑣)‖,𝜆=𝜆₁+𝜆₂,𝜖_𝑖=𝑦𝑆_𝑖⁽⁰⁾−𝑧𝑆_𝑖⁽⁰⁾𝑦𝑆_𝑖⁽¹⁾−𝑧𝑆_𝑖⁽¹⁾,𝑋_𝑖=max_{0≤𝑛≤𝑖+1}𝜖_𝑛,𝜖^′_𝑖=𝑦𝑆𝑡_𝑖+𝜙1,ℎ−𝑧𝑆𝑡_𝑖+𝜙1𝑦,ℎ𝑆𝑡_𝑖+𝜙2,ℎ−𝑧𝑆𝑡_𝑖+𝜙2,,ℎ (93) where 𝑡_{𝑖+𝜙𝑗}=𝑡_𝑖+𝜙_𝑗ℎ,𝑗=1,2.
It follows from (92), (93), (49), and (51) that 𝜖^′_𝑖≤𝑦₀−𝑧₀+ℎ^𝛼𝐿1Γ(𝛼)𝜆₁𝜖_𝑖+𝜆₂𝜖_𝑖+1𝛼+1Γ(𝛼)^𝑖−1_𝑛=0𝐶₀^𝛼−1(𝑖−𝑛)^𝛼−1×𝜆₁𝜖_𝑛+𝜆₂𝜖_𝑛+1≤𝑦₀−𝑧₀+ℎ^𝛼𝐿1Γ(𝛼+1)𝜆𝑋_𝑖+1Γ(𝛼)×^𝑖−1_𝑛=0𝐶₀^𝛼−1(𝑖−𝑛)^𝛼−1𝜆𝑋_𝑛≤𝑦₀−𝑧₀+ℎ^𝛼𝐿1Γ(𝛼+1)𝜆𝑋_𝑖+1Γ(𝛼)×^𝑖−1_𝑛=0𝐶₀^𝛼−1(𝑖−𝑛)^𝛼−1𝜆𝑋_𝑛. (94)
Take ‖𝜖_ℓ+1‖=max_{0≤𝑛≤𝑖+1}{‖𝜖_𝑛‖}=𝑋_𝑖; hence, 𝑋_𝑖=𝑋_ℓ,1≤ℓ≤𝑖. According to the fact that 𝐷(𝜙₁,𝜙₂) is invertible and ‖𝐷⁻¹(𝜙₁,𝜙₂)𝐶(𝜙₁,𝜙₂)‖<1, we have 𝑋_ℓ=𝜖_ℓ+1≤𝐷⁻¹𝜙₁,𝜙₂𝐶𝜙₁,𝜙₂𝜖_ℓ+𝐷⁻¹𝜙₁,𝜙₂𝜖^′_ℓ≤𝐷⁻¹𝜙₁,𝜙₂𝐶𝜙₁,𝜙₂𝑋_ℓ+𝐷⁻¹𝜙₁,𝜙₂𝜖^′_ℓ,1−𝐷⁻¹𝜙₁,𝜙₂𝐶𝜙₁,𝜙₂𝑋_ℓ≤𝐷⁻¹𝜙₁,𝜙₂𝜖^′_ℓ. (95)
Moreover, 𝑋_ℓ≤𝐶₁𝜖^′_ℓ, (96) where 𝐶₁=‖𝐷⁻¹(𝜙₁,𝜙₂)‖/(1−‖𝐷⁻¹(𝜙₁,𝜙₂)𝐶(𝜙₁,𝜙₂)‖). By means of (94), we have 𝑋_ℓ≤𝐶₁𝑦₀−𝑧₀+𝐶₁ℎ^𝛼𝐿1Γ(𝛼+1)𝜆𝑋_ℓ+1Γ(𝛼)^ℓ−1_𝑛=0𝐶₀^𝛼−1(ℓ−𝑛)^𝛼−1𝜆𝑋_𝑛,𝐶1−¹ℎ^𝛼𝐿𝜆Γ()𝑋𝛼+1^ℓ≤𝐶₁𝑦₀−𝑧₀+𝐶₁𝐿ℎ^𝛼1Γ(𝛼)^ℓ−1_𝑛=0𝐶₀^𝛼−1(ℓ−𝑛)^𝛼−1𝜆𝑋_𝑛. (97)(98)
Obviously, there exist 𝐶∈(0,1) and ℎ₁>0 such that 𝐶₁ℎ^𝛼𝐿𝜆/Γ(𝛼+1)≤𝐶<1 as ℎ<ℎ₁. Then 𝑋_ℓ≤𝐶₁𝐶𝑦1−⁰−𝑧₀+𝐶₁ℎ^𝛼𝐿Γ(𝛼)^ℓ−1_𝑛=0𝐶₀^𝛼−1(ℓ−𝑛)^𝛼−1𝑋_𝑛. (99) Let 𝐶=𝐶₁𝐶₀^𝛼−1𝐿𝜆/(1−𝐶)Γ(𝛼). We obtain 𝑋_ℓ≤𝐶₁𝐶𝑦1−⁰−𝑧₀+𝐶ℎ^𝛼ℓ−1_𝑛=0(ℓ−𝑛)^𝛼−1𝑋_𝑛. (100)
Applying Lemma 3 yields 𝑋_ℓ≤𝐶₁𝐶𝑦1−⁰−𝑧₀𝐸_𝛼𝐶Γ(𝛼)(ℓℎ)^𝛼≤𝐶₁𝐶𝑦1−⁰−𝑧₀𝐸_𝛼𝐶Γ(𝛼)𝑇^𝛼. (101) Using the convergence of Mittag-Leffler function 𝐸_𝛼(𝑧) and denoting 𝐶𝐶=¹𝐶𝐸1−^𝛼𝐶Γ(𝛼)𝑇^𝛼, (102) we obtain 𝜖_𝑖+1≤𝑋_𝑖=𝑋_ℓ≤𝐶𝑦₀−𝑧₀. (103)
Moreover, from the definitions of 𝑦𝑆(𝑡,ℎ) and 𝑧𝑆(𝑡,ℎ), we have 𝑦𝑆𝑡_𝑖+ℎ𝑢,ℎ−𝑧𝑆𝑡_𝑖𝑦+ℎ𝑢,ℎ𝑆𝑡_𝑖+ℎ𝑣,ℎ−𝑧𝑆𝑡_𝑖=+ℎ𝑣,ℎ𝑢²(2𝑢+1)𝑢²𝑣𝑣²(2𝑣+1)𝑣²𝑣𝑦𝑆_𝑖⁽⁰⁾−𝑧𝑆_𝑖⁽⁰⁾𝑦𝑆_𝑖⁽¹⁾−𝑧𝑆_𝑖⁽¹⁾+𝑢²(2𝑢+1)−𝑢²𝑢𝑣²(2𝑣+1)−𝑣²𝑣𝑦𝑆⁽⁰⁾_𝑖+1−𝑧𝑆⁽⁰⁾_𝑖+1𝑦𝑆⁽¹⁾_𝑖+1−𝑧𝑆⁽¹⁾_𝑖+1(𝑦=𝐶𝑢,𝑣)𝑆_𝑖⁽⁰⁾−𝑧𝑆_𝑖⁽⁰⁾𝑦𝑆_𝑖⁽¹⁾−𝑧𝑆_𝑖⁽¹⁾()𝑦+𝐷𝑢,𝑣𝑆⁽⁰⁾_𝑖+1−𝑧𝑆⁽⁰⁾_𝑖+1𝑦𝑆⁽¹⁾_𝑖+1−𝑧𝑆⁽¹⁾_𝑖+1,∀𝑢,𝑣∈[0,1]. (104) By means of (103), we get sup_{0≤𝑡≤𝑇}𝑦𝑆(𝑡,ℎ)−𝑧𝑆(𝑡,ℎ)≤max_{0≤𝑖≤𝑁−1}sup_{0≤𝑢<𝑣≤1}‖𝐶(𝑢,𝑣)‖𝜖_𝑖+sup_{0≤𝑢<𝑣≤1}‖𝐷(𝑢,𝑣)‖𝜖_𝑖+1≤𝜆₁+𝜆₂𝐶𝑦₀−𝑧₀≤𝐶^′𝑦₀−𝑧₀, (105) where 𝐶^′=(𝜆₁+𝜆₂)𝐶; 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 ‖error‖_𝐿1=ℎ^𝑁_𝑖=1𝑆𝑡_𝑖,ℎ−𝑦𝑡_𝑖,‖error‖_𝐿2=ℎ^𝑁_𝑖=1𝑆𝑡_𝑖,ℎ−𝑦𝑡_𝑖^21/2,order_𝐿1=log_1/2‖error‖_𝐿1(ℎ/2)‖error‖_𝐿1(ℎ),order_𝐿2=log_1/2‖error‖_𝐿2(ℎ/2)‖error‖_𝐿2(ℎ), (106) where ‖error‖(ℎ) is the error function which depends on the stepsize ℎ.

Example 1. Consider the initial value problem [10, 21] ^𝐶₀𝐷^𝛼_𝑡𝑦(𝑡)=−𝑦(𝑡)+𝑡²+2𝑡^2−𝛼Γ(),[]𝑦(0)3−𝛼𝑡∈𝐼=0,1,0<𝛼<1,=0. (107) The exact solution is 𝑦(𝑡)=𝑡².

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 ℎ=1/160; the absolute errors of DSCM, FWCM, and the method reported in [10] are shown in Table 1. When 𝛼=0.75, ℎ=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 (𝜙₁=0.85,𝜙₂=0.99) and the methods reported in [10, 21]; 𝛼=0.5.

Table 2 

The absolute errors of DSCM (𝜙₁=0.9,𝜙₂=0.98) and the methods reported in [10, 21]; 𝛼=0.75.

Figure 3 

The error between the numerical solutions and the exact solution of Example 1; 𝛼=0.75.

Example 2. Consider the initial value problem [12] ^𝐶₀𝐷^𝛼_𝑡𝑦(𝑡)=40320Γ()𝑡9−𝛼^8−𝛼Γ()−35+𝛼/2Γ()𝑡5−𝛼/2^4−𝛼/2+94Γ()+31+𝛼2𝑡^𝛼/2−𝑡⁴³−𝑦(𝑡)^3/2,[],𝑦(0)𝑡∈𝐼=0,1=0,0<𝛼≤1. (108) The exact solution is 𝑦(𝑡)=𝑡⁸−3𝑡^4+𝛼/2+(9/4)𝑡^𝛼. Obviously, 𝑦(𝑡)∉𝐶⁴[0,1]. In order to verify our theoretical results, we modify problem (108) as follows: ^𝐶₀𝐷^𝛼_𝑡𝑦(𝑡)=40320Γ()𝑡9−𝛼^8−𝛼Γ()−37+𝛼/2Γ()𝑡7−𝛼/2^6−𝛼/2+3Γ()𝑡325+𝛼⁴+32𝑡^2+𝛼/2−𝑡⁴³−𝑦(𝑡)^3/2,[],𝑦(0)𝑡∈𝐼=0,1=0,0<𝛼≤1. (109) The exact solution is 𝑦(𝑡)=𝑡⁸−3𝑡^6+𝛼/2+(9/4)𝑡^4+𝛼.

Tables 3 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 𝜙₁=0.9, 𝜙₂=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, 𝜙₁=0.88,𝜙₂=0.98.

Table 4 

The error between the numerical solutions and the exact solution of the problem (109); 𝛼=0.75, 𝜙₁=0.95,𝜙₂=0.99.

(a)

(b)

Figure 4 

The error between the numerical solutions and the exact solution of the problem (109); ℎ=0.05, 𝛼=0.000001.

(a)

(b)

Figure 5 

The error between the numerical solutions and the exact solution of the problem (109); 𝛼=0.9, ℎ=0.05.

Example 3. Consider the initial value problem ^𝐶₀𝐷^𝛼_𝑡𝑦(𝑡)+𝑡𝑦²(𝑡)(𝑡)[]𝑦(0)=𝑓,𝑡∈𝐼=0,𝑇,0<𝛼≤1,=𝜃₁ (110) and the perturbed problem ^𝐶₀𝐷^𝛼_𝑡𝑧(𝑡)+𝑡𝑧²(𝑡)(𝑡)[]𝑧(0)=𝑓,𝑡∈𝐼=0,𝑇,0<𝛼≤1,=𝜃₂. (111)

We use the cubic spline collocation method to solve the problems (110) and (111), respectively. Selecting 𝑓(𝑡)=𝑡sin(𝑡), 𝑇=20, 𝜃₁=0.1, 𝜃₂=0.2, we obtain the numerical results given in Figure 6. When 𝑓(𝑡)=𝑡sin(𝑡), 𝑇=20, 𝜃₁=2.0, 𝜃₂=−2.0, the numerical results are given in Figure 7. When 𝑓(𝑡)=𝑡𝑒^𝑡, 𝑇=20, 𝜃₁=1.0, 𝜃₂=2.0, the numerical results are shown in Figure 8.

(a)

(b)

Figure 6 

The numerical results of the problems (110) and (111); ℎ=0.5, 𝑇=20.0, 𝛼=0.8, 𝜃1=0.1, 𝜃2=0.2, 𝜙1=0.9, 𝜙2=0.98.

(a)

(b)

Figure 7 

The numerical results of the problems (110) and (111); ℎ=0.5, 𝛼=0.8, 𝜃1=2.0, 𝜃2=−2.0, 𝜙1=0.9, 𝜙2=0.98.

(a)

(b)

Figure 8 

The numerical results of the problems (110) and (111); ℎ=0.5, 𝛼=0.8, 𝜃1=1.0, 𝜃2=2.0, 𝜙1=0.9, 𝜙2=0.98.

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 𝑡 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.

Acknowledgments

This work is supported by projects from NSF of China (no. 11271311, no. 11226320, and no. 61144004), Program for Changjiang Scholars and Innovative Research Team in University of China (no. IRT1179), the Aid Program for Science and Technology Innovative Research Team in Higher Educational Institutions of Hunan Province of China, Huizhou Science and Technology Program (no. 20110103), and NSF of Huizhou University (no. 2012YB15).