Sage Journals: Discover world-class research

Abstract

The present contribution is concerned with the numerical solution of initial value problems for linear Caputo-type fractional differential equations. Global convergence estimates for a collocation method and its iterated version on special graded grids are derived. A numerical example confirming the theoretical results is also given.

Keywords

fractional differential equations Caputo fractional derivative collocation method iterated method graded grid

Introduction

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

Let $R = (- \infty, \infty)$ , $N = {1,2, \dots}$ . Let C^k(Ω) be the set of all k times (k ≥ 0) continuously differentiable functions on $Ω \subset R$ , C⁰(Ω) = C(Ω). In particular, by C[0, b] we denote the Banach space of continuous functions u on [0, b] with the typical norm ‖u‖_∞≔ max_0≤t≤b|u(t)|. Let α ∈ (0, 1). We consider the following initial value problem: $\begin{align} (D_{Cap}^{α} y) (t) + a (t) y (t) & = f (t), 0 \leq t \leq b, b > 0, \end{align}$ (1) $\begin{align} y (0) & = y_{0}, \end{align}$ (2)where $y_{0} \in R$ , a, f ∈ C[0, b] are some given continuous functions and $D_{Cap}^{α} y$ denotes the Caputo fractional derivative of order α of an unknown function y, defined as $(D_{Cap}^{α} y) (t) = \frac{1}{Γ (1 - α)} \frac{d}{d t} \int_{0}^{t} {(t - s)}^{- α} (y (s) - y (0)) d s, 0 < t \leq b .$

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

For given $b \in R, b > 0$ , $m \in N$ and 0 < ν < 1, by C^m,ν(0, b] we denote the set of continuous functions $u : [0, b] \to R$ which are m times continuously differentiable in (0, b] such that $| {| u}^{(i)} (t) | | \leq c t^{1 - ν - i}, t \in (0, b], i = 1, \dots, m,$ where c = c(u) is a positive constant independent of t. Equipped with the norm $‖ u ‖_{C^{m, ν} (0, b]} ≔ ‖ u ‖_{\infty} + \sum_{i = 1}^{m} \sup_{0 < t \leq b} t^{i - 1 + ν} |u^{(i)} (t)|, u \in C^{m, ν} (0, b],$ the set C^m,ν(0, b] becomes a Banach space. Note that $C^{q} [0, b] \subset C^{q, ν} (0, b] \subset C^{m, μ} (0, b] \subset C [0, b], q \geq m \geq 1, 0 < ν \leq μ < 1 .$

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

Theorem 1

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

Remark 1

Note that the statement y ∈ C^q,ν(0, b] for the solution y of (1)–(2) in Theorem 1 cannot, in general, be improved. For example, let us consider the following initial value problem: $(D_{Cap}^{\frac{1}{2}} y) (t) = \frac{\sqrt{π}}{2}, t \in [0,1]; y (0) = 1 .$

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

It is easy to see that problem (1)–(2) can be reformulated as an equivalent Volterra integral equation for y: $y (t) = y_{0} + \frac{1}{Γ (α)} \int_{0}^{t} {(t - s)}^{α - 1} [f (s) - a (s) y (s)] d s, 0 \leq t \leq b,$ or, in operator form, $y = T y + g,$ (3)where $\begin{align} (T y) (t) & = - \frac{1}{Γ (α)} \int_{0}^{t} {(t - s)}^{α - 1} a (s) y (s) d s, 0 \leq t \leq b, \end{align}$ (4) $\begin{align} g (t) & = y_{0} + \frac{1}{Γ (α)} \int_{0}^{t} {(t - s)}^{α - 1} f (s) d s, 0 \leq t \leq b . \end{align}$ (5)

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

Collocation method

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

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

Next, for a given integer k ≥ 0, we introduce $S_{k}^{(- 1)} (Π_{N}) ≔ \{v : v {| |}_{[t_{j - 1}, t_{j}]} \in π_{k}, j = 1, \dots, N\} .$

Here $v | |_{[t_{j - 1}, t_{j}]}$ is the restriction of $v : [0, b] \to R$ onto the subinterval [t_j−1, t_j] ⊂ [0, b] and π_k denotes the set of polynomials of degree not exceeding k. Note that the elements of $S_{k}^{(- 1)} (Π_{N})$ may have jump discontinuities at the interior points t₁, …, t_N−1 of the grid Π_N. We choose $m \in N$ points η₁, …, η_m in the interval [0, 1] so that $0 \leq η_{1} < \dots < η_{m} \leq 1$ (7)and define in every subinterval [t_j−1, t_j] ⊂ [0, b] the collocation points $t_{j p} = t_{j - 1} + η_{p} (t_{j} - t_{j - 1}), p = 1, \dots, m, j = 1, \dots, N .$ (8)

Approximations $y_{N} \in S_{m - 1}^{(- 1)} (Π_{N})$ to the solution y of (3) can be found by assuming that the following conditions hold: $y_{N} (t_{j p}) = (T y_{N}) (t_{j p}) + g (t_{j p}), p = 1, \dots, m, j = 1, \dots, N,$ (9)where T and g are defined by (4) and (5), respectively. The collocation conditions (9) have an operator equation representation $y_{N} = P_{N} T y_{N} + P_{N} g,$ (10)where $P_{N} = P_{N, m} : C [0, b] \to S_{m - 1}^{(- 1)} (Π_{N})$ is defined by $(P_{N} v) (t_{j p}) = v (t_{j p}), p = 1, \dots, m, j = 1, \dots, N, v \in C [0, b] .$

If η₁ = 0, then by (P_Nv)(t_j1) we denote the right limit $\lim_{t \to t_{j - 1}^{+}} (P_{N} v) (t)$ . If η_m = 1, then by (P_Nv)(t_jm) we denote the left limit $\lim_{t \to t_{j}^{-}} (P_{N} v) (t)$ .

The collocation conditions (9) form a system of equations whose exact form is determined by the choice of a basis in the space $S_{m - 1}^{(- 1)} (Π_{N})$ . We will use the Lagrange fundamental polynomial representation: $y_{N} (t) = \sum_{l = 1}^{N} \sum_{k = 1}^{m} c_{l k} φ_{l k} (t), t \in [0, b],$ (11)where $φ_{l k} (t) = \{\begin{cases} 0 & f o r t \notin [t_{l - 1}, t_{l}], \\ \prod_{i = 1, i \neq k}^{m} \frac{t - t_{l i}}{t_{l k} - t_{l i}} f o r & t \in [t_{l - 1}, t_{l}], \end{cases} k = 1, \dots, m, l = 1, \dots, N .$

As usual, we take $\prod_{i = 1, i \neq k}^{m} (t - t_{l i}) / (t_{l k} - t_{l i}) = 1$ if m = 1. Note that y_N(t_lk) = c_lk, k = 1, …, m, l = 1, …, N, and using representation (11) the collocation conditions (9) form the following system of linear algebraic equations with respect to the unknown values {c_jp}: $c_{j p} = \sum_{l = 1}^{N} \sum_{k = 1}^{m} (T φ_{l k}) (t_{j p}) c_{l k} + g (t_{j p}), p = 1, \dots, m, j = 1, \dots, N .$

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

Theorem 2

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

Then there exists an integer N₀ such that for all N ≥ N₀ equation (10) possesses a unique solution $y_{N} \in S_{m - 1}^{(- 1)} (Π_{N})$ and $\sup_{0 \leq t \leq b} | y (t) - y_{N} (t) | \to 0 a s N \to \infty,$ where y is the solution of (1)–(2).

(ii) If, in addition to (i), we assume that a, f ∈ C^m,μ(0, b], where 0 < μ < 1, then for all N ≥ N₀ the error estimate $\sup_{0 \leq t \leq b} | y (t) - y_{N} (t) | \leq c \{\begin{cases} N^{- r (1 - ν)} & f o r 1 \leq r < \frac{m}{1 - ν}, \\ N^{- m} & f o r r \geq \frac{m}{1 - ν}, \end{cases}$ (12)holds. Here ν = max{1 − α, μ}, r ≥ 1 is the grading parameter of the grid (6) and c is a constant independent of N.

Iterated numerical method

Theorem 2 shows that if the functions a and f are smooth enough, and if the grading parameter r is chosen to be sufficiently large, the expected convergence order of the proposed numerical method is of order O(N^−m). This convergence order can be improved by using an iterated method, assuming a little more smoothness of a and f, together with a more precise choice of collocation parameters η₁, …, η_m. More precisely, let y_N be the solution of (10), with N ≥ N₀ (see Theorem 2). Denote $y_{N}^{i t} ≔ T y_{N} + g, N \geq N_{0},$ (13)where T and g are defined by (4) and (5), respectively. Note that $y_{N}^{i t} \in C [0, b]$ . Based on Theorems 1 and 2, and using some ideas of^3,7,14 we can prove the following result, which generalizes and refines the corresponding error estimates in²⁰.

Theorem 3

Let α ∈ (0, 1) and $m \in N$ . Assume that a ∈ C^m+1[0, b] and f ∈ C^m+1,μ(0, b], where 0 < μ < 1. Let $N \in N$ and assume that the grid points (6) and the collocation parameters η₁, …, η_m satisfying (7) are used. Moreover, assume that the collocation parameters are chosen so that the quadrature approximation $\int_{0}^{1} F (x) d x \approx \sum_{k = 1}^{m} w_{k} F (η_{k})$ (14)with appropriate weights {w_k} is exact for all polynomials F of degree m.

Then there exists $N_{0} \in N$ so that for N ≥ N₀ the approximation $y_{N}^{i t}$ defined by (13) is unique and the error estimate $‖ y - y_{N}^{i t} ‖_{\infty} \leq c \{\begin{cases} N^{- r (1 + α - ν)} & f o r 1 \leq r < \frac{m + α}{1 + α - ν}, \\ N^{- m - α} & f o r r \geq \frac{m + α}{1 + α - ν}, \end{cases}$ (15)holds. Here y is the exact solution of problem (1)–(2), ν = max{1 − α, μ}, r ≥ 1 is the grading parameter of the grid (6) and c is a constant independent of N.

Numerical example

Consider the following initial value problem: $(D_{Cap}^{\frac{1}{2}} y) (t) + (1 - t) y (t) = t^{\frac{1}{2}} - t^{\frac{3}{2}} - 2 t + 2 + Γ (\frac{3}{2}), 0 \leq t \leq 1; y (0) = 2 .$ (16)

We see that (16) is a special problem of (1)–(2) with $α = \frac{1}{2}, a (t) = 1 - t, f (t) = t^{\frac{1}{2}} - t^{\frac{3}{2}} - 2 t + 2 + Γ (\frac{3}{2}), y_{0} = 2, b = 1 .$

Note that a ∈ C^m[0, 1] and f ∈ C^m,1/2(0, 1] for all $m \in N$ . Therefore by Theorem 1 it follows that y ∈ C^m,ν(0, b], where $ν = \max \{1 - \frac{1}{2}, \frac{1}{2}\} = \frac{1}{2} .$ (17)

Special case of collocation parameters

We apply the collocation method (see Section 2) and its iterated version (see Section 3) with m = 2, using as collocation parameters the shifted Gauss-Legendre points $η_{1} = \frac{3 - \sqrt{3}}{6}, η_{2} = 1 - η_{1},$ which satisfy the conditions set for collocation parameters in Theorem 3. It follows from Theorem 2 that for sufficiently large $N \in N$ we have $\sup_{0 \leq t \leq b} | y (t) - y_{N} (t) | \leq c_{0} \{\begin{cases} N^{- 0.5 r} & i f 1 \leq r < 4, \\ N^{- 2} & i f r \geq 4, \end{cases}$ (18)where $y (t) = t^{\frac{1}{2}} + 2, 0 \leq t \leq 1,$ is the exact solution of problem (16) and c₀ is a positive constant independent of N. Similarly, it follows from Theorem 3 that $‖ y - y_{N}^{i t} ‖_{\infty} \leq c_{1} \{\begin{cases} N^{- r} & i f 1 \leq r < \frac{5}{2}, \\ N^{- 2.5} & i f r \geq \frac{5}{2}, \end{cases}$ (19)where c₁ is a positive constant independent of N.

In Table 1 some results of numerical experiments for different values of the parameter r are presented. The actual numerical error ɛ_N for sup_0≤t≤b|y(t) − y_N(t)| and

ε_{N}^{i t}

for

‖ y - y_{N}^{i t} ‖_{\infty}

are calculated as follows:

\begin{array}{l} ε_{N} & ≔ \max_{j = 1, \dots, N} \max_{k = 0, \dots, 10} | y (τ_{j k}) - y_{N} (τ_{j k}) |, \\ ε_{N}^{i t} & ≔ \max_{j = 1, \dots, N} \max_{k = 0, \dots, 10} | y (τ_{j k}) - y_{N}^{i t} (τ_{j k}) |, \end{array}

where

τ_{j k} ≔ t_{j - 1} + k (t_{j} - t_{j - 1}) / 10, k = 0, \dots, 10, j = 1, \dots, N,

with t_j defined by (6). The ratios

ϱ_{N} ≔ \frac{ε_{N / 2}}{ε_{N}}, ϱ_{N}^{i t} ≔ \frac{ε_{N / 2}^{i t}}{ε_{N}^{i t}},

characterizing the observed convergence rates, are also presented. Due to (18), the ratios ϱ_N for the non-iterated collocation method for r = 1, r = 2 and r = 4 ought to be approximately 2^0.5 ≈ 1.414, 2¹ = 2 and 2² = 4, respectively. Due to (19), the ratios

ϱ_{N}^{i t}

for the iterated method for r = 1, r = 2 and r = 4 ought to be approximately 2¹ = 2, 2² = 4 and 2^2.5 = 5.667, respectively. All these ratios are also given in the corresponding rows of Table 1. As we can see from Table 1 the numerical results are in good accord with the theoretical estimates given by Theorem 2 and Theorem 3 (the estimates (12) and (15)).

Table 1.

Numerical results for shifted Gauss points.

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

General case of collocation parameters

We also briefly highlight the case when the chosen collocation parameters do not satisfy the quadrature condition set in Theorem 3, by applying both methods (for m = 2) with collocation parameters $η_{1} = 0.1, η_{2} = 0.9 .$

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

Table 2.

Numerical results for non-Gauss points.

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

Statements and declarations

Footnotes

Acknowledgments

We thank the reviewers for their remarks and comments.

Conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research,authorship,and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research,authorship,and/or publication of this article: This research was supported by Estonian Research Council grant PRG864.

References

Baleanu

Diethelm

Trujillo

. Fractional calculus. Models and numerical methods. Singapore: World Scientific Publishing Co. Pte. Ltd, 2016.

Brugnano

Burrage

, et al. A spectrally accurate step-by-step method for the numerical solution of fractional differential equations. J Sci Comput 2024; 99: 48.

Brunner

. Collocation methods for Volterra integral and related functional differential equations. Cambridge: Cambridge University Press, 2004.

Cardone

Conte

Paternoster

. Stability of two-step spline collocation methods for initial value problems for fractional differential equations. Commun Nonlinear Sci Numer Simul 2022; 115: 106726.

Diethelm

. The analysis of fractional differential equations. Berlin: Springer, 2010.

Ford

Morgado

Rebelo

. A nonpolynomial collocation method for fractional terminal value problems. J Comput Appl Math 2015; 275: 392–402.

Ford

Pedas

Vikerpuur

. High order approximations of solutions to initial value problems for linear fractional integro-differential equations. Fract Calc Appl Anal 2023; 26: 2069–2100.

Garrappa

. Numerical solution of fractional differential equations: a survey and a software tutorial. Math 2018; 6: 16.

Kolk

Pedas

Tamme

. Modified spline collocation for linear fractional differential equations. J Comput Appl Math 2015; 283: 28–40.

10.

Kolk

Pedas

Tamme

. Smoothing transformation and spline collocation for linear fractional boundary value problems. Appl Math Comput 2016; 283: 234–250.

11.

Liang

Stynes

. Collocation methods for general Caputo two-point boundary value problems. J Sci Comput 2018; 76: 390–425.

12.

Lillemäe

Pedas

Vikerpuur

. Central part interpolation schemes for a class of fractional initial value problems. Acta Commentationes Univ Tartuensis Math 2022; 26: 161–178.

13.

Lillemäe

Pedas

Vikerpuur

. Central part interpolation schemes for fractional differential equations. Appl Numer Math 2024; 200: 318–330.

14.

Pedas

Tamme

. Spline collocation methods for linear multi-term fractional differential equations. J Comput Appl Math 2011; 236: 167–176.

15.

Pedas

Vikerpuur

. Spline collocation for multi-term fractional integro-differential equations with weakly singular kernels. Fractal Fract 2021; 5: 90.

16.

Podlubny

. Fractional differential equations. Cambridge, MA: Academic Press, 1999.

17.

Timofejeva

Navickas

Telksnys

, et al. An operator-based scheme for the numerical integration of FDEs. Mathematics 2021; 9: 1372.

18.

Vainikko

. Multidimensional weakly singular integral equations. Berlin: Springer, 1993.

19.

Vikerpuur

. Two collocation type methods for fractional differential equations with non-local boundary conditions. Math Model Anal 2017; 22: 654–670.

20.

Wang

Liang

. Superconvergence and postprocessing of collocation methods for fractional differential equations. J Sci Comput 2023; 97: 29.