Capturing Multiple Singularities with Spectral Accuracy for Multi-Term Fractional Differential Equations

1. Introduction

Fractional calculus has attracted considerable attention over the past few decades due to its intrinsic non-locality and memory effects, which provide a more accurate modeling framework than traditional integer-order calculus for a wide range of complex dynamical systems [1,2,3]. As the theoretical foundation has matured, fractional-order models have found applications in viscoelastic constitutive relations [4], bio-tissue electrical conductivity [5], financial dynamics [6], and advanced control theory [7].

In many multi-physics or multi-scale processes, multi-term fractional differential equations arise naturally. They offer a faithful description of systems with multiple relaxation times [4] and multi-rate anomalous diffusion in heterogeneous media. However, their numerical solution faces significant challenges. First, the solutions of multi-term fractional ordinary differential equations (FODEs) typically lack smoothness: due to the combined influence of several fractional derivatives with different orders, the solution may exhibit complicated initial layer behavior and multiple singularities [8]. Second, the convolutional non-locality of fractional derivatives leads to dense or densely-structured discretization matrices, resulting in computational and memory costs of at least $O\left(N^2\right)$ [9].

A more subtle but critical difficulty emerges when high-order algorithms, in particular spectral methods [10,11], are applied to multi-term FODEs. While spectral methods can achieve exponential convergence for single-term fractional equations, their direct extension to multi-term problems often yields disappointingly low accuracy [12,13,14]. The fundamental reason is that standard polynomial or trigonometric basis functions are ill-equipped to represent simultaneously the complex singular structures induced by multiple distinct fractional kernels.

To overcome these difficulties, various numerical strategies have been developed. Finite difference methods [2] and predictor-corrector schemes [8] are straightforward to implement but suffer from limited accuracy. Finite element methods can handle complex geometries but typically do not achieve high-order convergence. Spectral methods [11,15], leveraging the exponential convergence of global polynomial approximations, remain a mainstream choice for high-precision solutions of fractional differential equations. To remedy the loss of spectral accuracy caused by solution singularities, a promising direction is to construct basis functions that inherently capture the singular behavior. Specific proposals include using fractional orthogonal polynomials (e.g., Jacobi poly-fractonomials) [12,15,16] or explicitly incorporating singular functions into standard polynomial bases [12,17,18]. When the basis singularity exactly matches that of the exact solution, exponential convergence can in principle be restored.

In recent years, frame theory has emerged as a paradigm that transcends traditional bases. The seminal works by Adcock, Huybrechs, and co-authors [19,20,21,22] have systematically developed numerical frame theory, demonstrating that even when the resulting linear system is severely ill-posed, frame approximations can maintain stability and high accuracy through appropriate regularization. This flexibility offers a new avenue for addressing the singularity-induced ill-conditioning that plagues spectral methods for multi-term FODEs.

Motivated by the limitations of existing spectral approaches in handling multiple singularities and the flexibility offered by frame theory, this paper proposes a fractional polynomial frame collocation method (FPFCM) for solving general multi-term FODEs. Specifically, we construct a hybrid basis function frame, derive an exact recursive scheme for evaluating Caputo fractional derivatives within this hybrid space, and discretize the original differential equation via a collocation approach. To ensure numerical stability, we incorporate truncated singular value decomposition (TSVD) regularization. Through extensive numerical experiments on linear and nonlinear multi-term FODEs, including the fractional Bagley–Torvik equation, we validate the effectiveness, high accuracy, and superconvergence properties of the proposed method.

The main contributions of this work are as follows:

• We construct a fractional polynomial frame tailored to functions with multiple singularities and develop a corresponding collocation method for multi-term FODEs. Unlike existing spectral methods that rely on a single type of basis, the over-complete frame can adaptively capture several distinct singular exponents.
• We employ an oversampling strategy combined with TSVD regularization to stabilize the discrete system, which is severely ill-conditioned due to the near-linear dependence of the frame functions. This regularization preserves spectral accuracy while ensuring robustness.
• We provide an efficient algorithm for computing the Caputo derivative of the frame elements using a three-term recurrence relation, enabling fast assembly of the collocation matrix.
• We numerically demonstrate exponential convergence and, remarkably, a superconvergence phenomenon when a tunable parameter in the basis matches an exact singularity exponent. The performance is validated on the fractional Bagley–Torvik equation as well as on linear and nonlinear multi-term FODEs, showing clear advantages over standard polynomial spectral methods.

The remainder of this paper is organized as follows. Section 2 reviews essential preliminaries on the Caputo fractional derivative, frame theory, and the class of multi-term FODEs under consideration. Section 3 details the construction of the fractional polynomial frame and presents approximation results. Section 4 describes the fractional polynomial frame collocation method, the TSVD regularization strategy, and the fast, stable evaluation of fractional derivatives of frame functions. Section 5 presents numerical examples demonstrating the effectiveness and superiority of the proposed method, including the fractional Bagley–Torvik equation, linear and nonlinear multi-term FODEs. Finally, Section 6 concludes the paper.

2. Preliminaries

We first give the definition of Caputo fractional derivative and some of its properties, which is used extensively in diverse fields. Then, we present the frame concept and some examples of frame. At last, we introduce the concerned problem.

 Definition 1

([3,9]). Let $f\left(t\right)$ be sufficiently smooth on $[0,T]$, $\gamma >0,n-1<\gamma \le n$, and $n\in \mathbb{N}$. Then the γ-th order Caputo fractional derivative of $f\left(t\right)$ is defined as:

$$D^{\gamma }f\left(t\right)={\displaystyle \frac{1}{\mathrm{\Gamma }(n-\gamma )}}{\int }_0^t{\displaystyle \frac{f^{\left(n\right)}\left(\tau \right)}{{(t-\tau )}^{\gamma -n+1}}}\mathrm{d}\tau ,$$

where $\mathrm{\Gamma }(·)$ is the Euler’s Gamma function, defined as

$$\mathrm{\Gamma }\left(z\right)={\int }_0^{\infty }{t}^{z-1}{e}^{-t}\mathrm{d}t.$$

Some properties of the Caputo fractional derivative are listed as follows [9,23]:

• If A is a constant, $D^{\gamma }A=0.$
• For a polynomial function $t^m$, there holds

  $$D^{\gamma }{t}^m=\left\{\begin{array}{c}0,m\in \mathbb{N}\cup \left\{0\right\}\mathrm{and}m<\lceil \gamma \rceil \hfill \\ {\displaystyle \frac{\mathrm{\Gamma }(m+1)}{\mathrm{\Gamma }(m+1-\gamma )}}{t}^{m-\gamma },m\in \mathbb{N}\cup \left\{0\right\}\mathrm{and}m\ge \lceil \gamma \rceil \mathrm{or}m\notin \mathbb{N}\mathrm{and}m>\lfloor \gamma \rfloor ,\hfill \end{array}\right.$$

  where $\lceil \gamma \rceil$ denotes the smallest integer greater than or equal to $\gamma$ (ceiling function), and $\lfloor \gamma \rfloor$ denotes the largest integer less than or equal to $\gamma$ (floor function).
• If $\beta >\gamma -1$, then

  $$D^{\gamma }\left(t^{\beta }lnt\right)={\displaystyle \frac{\mathrm{\Gamma }(\beta +1)}{\mathrm{\Gamma }(\beta -\gamma +1)}}{t}^{\beta -\gamma }\left(lnt+\psi (\beta +1)-\psi (\beta -\gamma +1)\right),$$

  (1)

  where $\psi \left(z\right):={\displaystyle \frac{{\mathrm{\Gamma }}^{\prime }\left(z\right)}{\mathrm{\Gamma }\left(z\right)}}$ is the digamma function.
• The Caputo fractional derivative possesses the property of linear operation:

  $$D^{\gamma }\left(\sum _{i=1}^n{c}_i{f}_i\left(t\right)\right)=\sum _{i=1}^n{c}_i{D}^{\gamma }{f}_i\left(t\right),$$

  where $c_i(i=1,...,n)$ are constant numbers.

As a generalization of orthonormal basis, frame is flexible and easy to construct, which breaks the strict constraints of linear independence and uniqueness of representation imposed by traditional bases, allowing redundancy. The definition of frame in a Hilbert space is as follows.

 Definition 2

([19,24]). Let $\mathcal{H}$ be a Hilbert space with inner product $(·,·)$. If a set $\mathrm{\Phi }={\left\{{\varphi }_n\right\}}_{n\in I}$ satisfies the frame condition:

$${A\parallel f\parallel }^2\le \sum _{n\in I}{\left|(f,{\varphi }_n)\right|}^2\le B{\parallel f\parallel }^2,\forall f\in \mathcal{H},$$

where $A,B>0$ are frame bounds, then Φ is called a frame of $\mathcal{H}$. If $A=B$, the frame is tight.

To demonstrate the specific applications of the frames in numerical approximation, three typical frame construction examples [19] are presented below:

• Fourier frames for complex geometries: Assume that $\mathrm{\Omega }\subset {(-1,1)}^d$ is a complex geometric domain (not a simple rectangle/cube), by restricting the orthonormal Fourier basis on ${(-1,1)}^d$ to $\mathrm{\Omega }$, one has

  $$\mathrm{\Phi }={\left\{{\varphi }_{\mathbf{n}}\right\}}_{\mathbf{n}\in {\mathbb{Z}}^d},{\varphi }_{\mathbf{n}}\left(\mathbf{t}\right)=2^{-d/2}{e}^{i\pi \mathbf{n}·\mathbf{t}},\mathbf{t}\in \mathrm{\Omega }.$$
  This system is a tight frame of $L^2\left(\mathrm{\Omega }\right)$ and is linearly independent, with frame bounds $A=B=1$. The Fourier frame is employed to numerically solve problem on complex geometric domain [19,25].
• Augmented Fourier basis: In the Hilbert space $L^2(-1,1)$, by adding a finite number of polynomials to the basis of the standard Fourier basis, one has

  $$\mathrm{\Phi }={\left\{{\phi }_n\right\}}_{n\in \mathbb{Z}}\cup {\left\{{\psi }_k\right\}}_{k=1}^K,{\phi }_n\left(t\right)={\displaystyle \frac{1}{\sqrt{2}}}{e}^{i\pi nt},$$

  where ${\psi }_k=\sqrt{k+1/2}{P}_k$, and $P_k$ is the k-th order Legendre polynomial. $\mathrm{\Phi }$ forms a frame for $L^2(-1,1)$, with frame bounds $A=1,B=2$, and it is linearly independent. The frame aims to approximate the nonperiodic and oscillatory functions.
• Polynomials plus modified polynomials: Assume that $w\in L^{\infty }(-1,1)$ may be singular, oscillatory, or possess other features that make approximation difficult. Let ${\left\{{\phi }_n\right\}}_{n\in \mathbb{N}}$ be the orthonormal basis of Legendre polynomials, the system defined as

  $$\mathrm{\Phi }={\left\{{\phi }_n\right\}}_{n\in \mathbb{N}}\cup {\left\{{\psi }_n\right\}}_{n\in \mathbb{N}},{\psi }_n\left(t\right)=w\left(t\right){\phi }_n\left(t\right)$$

  gives rise to a frame for the space $L^2\left([-1,1]\right)$. The frame bounds are

  $$A=1+{\mathrm{essinf}}_{t\in (-1,1)}{\left|w\left(t\right)\right|}^2,B=1+{\mathrm{esssup}}_{t\in (-1,1)}{\left|w\left(t\right)\right|}^2.$$
  This frame is linearly dependent if and only if w is a rational function of two polynomials. We are interested in this frame since the fractional derivative brings a weak endpoint singularity like $w\left(t\right)={(1+t)}^{\delta }$.

Our method is designed for solving multi-term FODEs, following the last typical frame described above. To fix the idea, let $0<{\gamma }_1<\dots <{\gamma }_S$, we consider the Cauchy problem of the multi-term FODE:

$$F\left(t,y\left(t\right),D^{{\gamma }_1}y\left(t\right),\dots ,D^{{\gamma }_S}y\left(t\right)\right)=0,t\in (0,T]$$

(2)

with initial value conditions:

$$y^{\left(k\right)}\left(0\right)=a_k,k=0,1,\dots ,⌈{\gamma }_S⌉-1,$$

(3)

where ${\gamma }_S$ denotes the order of the fractional derivative, and F is a known (linear or nonlinear) function of $t,y$ and the fractional-order derivatives $D^{{\gamma }_k}y\left(t\right)(1\le k\le S)$. In practice, our method is easy to extend to the boundary value problems and even the FODE system.

3. Fractional Polynomial Frame

3.1. Shifted Fractional Jacobi and Legendre Polynomials

Let ${\mathbb{P}}_M\left(\Lambda \right)$ be the collection of all algebraic polynomials on the interval $\Lambda$ of degree at most M. For $\alpha ,\beta >-1$, the classical Jacobi polynomial system ${\left\{P_i^{(\alpha ,\beta )}\left(x\right)\right\}}_{i=0}^{\infty }$ forms an orthogonal basis of the space $L_{{\omega }^{(\alpha ,\beta )}}^2(-1,1)$ with the inner product as

$${(u,v)}_{{\omega }^{(\alpha ,\beta )}}:={\int }_{-1}^{1}u\left(x\right)v\left(x\right){(1-x)}^{\alpha }{(1+x)}^{\beta }\mathrm{d}x.$$

The weight function is ${\omega }^{(\alpha ,\beta )}\left(x\right)={(1-x)}^{\alpha }{(1+x)}^{\beta }$ and the induced norm is denoted by ${\parallel ·\parallel }_{{\omega }^{(\alpha ,\beta )}}$. We will drop the subscript when $\alpha =\beta =0$ since ${\omega }^{(0,0)}=1$.

For $\alpha ,\beta >-1$, there holds some important properties of the Jacobi polynomial $\left\{P_i^{\alpha ,\beta }\left(x\right)\right\}$, that are presented for readers in Appendix A. The classical Legendre polynomial is the special case $\alpha =\beta =0$ of the Jacobi polynomial $P_i^{(\alpha ,\beta )}\left(x\right)$, denote by $P_i\left(x\right):=P_i^{(0,0)}\left(x\right)$.

The family of shifted fractional Jacobi polynomials is defined for $\delta \ge 0$, as in [26]:

$${\varphi }_i^{(\delta ,\alpha ,\beta )}\left(t\right):=t^{\delta }{P}_i^{(\alpha ,\beta )}\left({\displaystyle \frac{2}{T}}t-1\right),i\in {\mathbb{N}}_0,t\in [0,T].$$

(4)

Some properties on the shifted fractional Jacobi polynomials are very useful.

 Theorem 1. 

For $\alpha ,\beta >-1,\delta \ge 0,t\in (0,T]$ and $i\in {\mathbb{Z}}_0$, there hold

(i)

$t{\varphi }_i^{(\delta ,\alpha ,\beta )}\left(t\right)={\varphi }_i^{(\delta +1,\alpha ,\beta )}\left(t\right).$

(ii)

The first-order derivative:

$${\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left[{\varphi }_i^{(\delta ,\alpha ,\beta )}\left(t\right)\right]=\delta {\varphi }_i^{(\delta -1,\alpha ,\beta )}\left(t\right)+{\displaystyle \frac{i+\alpha +\beta +1}{T}}{\varphi }_{i-1}^{(\delta ,\alpha +1,\beta +1)}\left(t\right),$$

and the k-order derivative:

$$\begin{array}{cc}\hfill {\displaystyle \frac{{\mathrm{d}}^k}{\mathrm{d}{t}^k}}\left[{\varphi }_i^{(\delta ,\alpha ,\beta )}\left(t\right)\right]=& \sum _{j=0}^k{d}_{j,i}^{\delta ,k}{\varphi }_{i-k+j}^{(\delta -j,\alpha +k-j,\beta +k-j)}\left(t\right),\hfill \\ \hfill d_{j,i}^{\delta ,k}=& {\displaystyle \frac{\mathrm{\Gamma }(k+1)\mathrm{\Gamma }(\delta +1)\mathrm{\Gamma }(i+k-j+\alpha +\beta +1)}{T^{k-j}\mathrm{\Gamma }(j+1)\mathrm{\Gamma }(k-j+1)\mathrm{\Gamma }(\delta -j+1)\mathrm{\Gamma }(i+\alpha +\beta +1)}}\hfill \end{array}$$

(5)

with ${\varphi }_i^{(\delta ,\alpha ,\beta )}\left(t\right)=0$ if $i<0$.

(iii)

The recurrence relation:

$$\begin{array}{cc}\hfill {\varphi }_{i+1}^{(\delta ,\alpha ,\beta )}\left(t\right)& ={\displaystyle \frac{2A_i^{\alpha ,\beta }}{T}}{\varphi }_i^{(\delta +1,\alpha ,\beta )}\left(t\right)-(A_i^{\alpha ,\beta }+B_i^{\alpha ,\beta }){\varphi }_i^{(\delta ,\alpha ,\beta )}\left(t\right)-C_i^{\alpha ,\beta }{\varphi }_{i-1}^{(\delta ,\alpha ,\beta )}\left(t\right),i\ge 1,\hfill \\ \hfill {\varphi }_1^{(\delta ,\alpha ,\beta )}\left(t\right)& ={\displaystyle \frac{\alpha +\beta +2}{T}}{t}^{\delta +1}-(\beta +1)t^{\delta },{\varphi }_0^{(\delta ,\alpha ,\beta )}\left(t\right)=t^{\delta },\hfill \end{array}$$

(6)

with $A_i^{\alpha ,\beta },B_i^{\alpha ,\beta },C_i^{\alpha ,\beta }$ as in (A2).

(iv)

The integral relation:

$${\displaystyle \frac{1}{\mathrm{\Gamma }\left(\gamma \right)}}{\int }_0^t{\displaystyle \frac{{\varphi }_i^{(\delta ,\alpha ,\delta )}\left(\tau \right)}{{(t-\tau )}^{1-\gamma }}}\mathrm{d}\tau ={\displaystyle \frac{\mathrm{\Gamma }(i+\delta +1)}{\mathrm{\Gamma }(i+\delta +\gamma +1)}}{\varphi }_i^{(\delta +\gamma ,\alpha -\gamma ,\delta +\gamma )}\left(t\right).$$

(7)

Proof From the definition (4), the equality in (i) holds clearly. Making use of (A3) in Appendix A, the derivative formula in (ii) can be obtained by direct calculation.

As for the recurrence relation in (iii), it can be derived from the definition (4), (i) and the three-term recurrence relation (A1) in Appendix A.

We now show the integral relation (7). By transformation $s={\displaystyle \frac{2}{T}}\tau -1$ and $x={\displaystyle \frac{2}{T}}t-1$, one has ${\varphi }_i^{(\delta ,\alpha ,\delta )}\left(\tau \right)={\left({\displaystyle \frac{T}{2}}\right)}^{\delta }{(1+s)}^{\delta }{P}_i^{(\alpha ,\delta )}\left(s\right)$ and ${(t-\tau )}^{1-\gamma }={\left({\displaystyle \frac{T}{2}}\right)}^{1-\gamma }{(x-s)}^{1-\gamma }$, then

$${\displaystyle \frac{1}{\mathrm{\Gamma }\left(\gamma \right)}}{\int }_0^t{\displaystyle \frac{{\varphi }_i^{(\delta ,\alpha ,\delta )}\left(\tau \right)}{{(t-\tau )}^{1-\gamma }}}\mathrm{d}\tau ={\left({\displaystyle \frac{T}{2}}\right)}^{\delta +\gamma }{\displaystyle \frac{1}{\mathrm{\Gamma }\left(\gamma \right)}}{\int }_{-1}^x{\displaystyle \frac{{(1+s)}^{\delta }{P}_i^{(\alpha ,\delta )}\left(s\right)}{{(x-s)}^{1-\gamma }}}\mathrm{d}s.$$

Hence, from (A6) in Appendix A, we obtain

$${\displaystyle \frac{1}{\mathrm{\Gamma }\left(\gamma \right)}}{\int }_0^t{\displaystyle \frac{{\varphi }_i^{(\delta ,\alpha ,\delta )}\left(\tau \right)}{{(t-\tau )}^{1-\gamma }}}\mathrm{d}\tau ={\left({\displaystyle \frac{T}{2}}\right)}^{\delta +\gamma }{\displaystyle \frac{\mathrm{\Gamma }(i+\delta +1)}{\mathrm{\Gamma }(i+\delta +\gamma +1)}}{(1+x)}^{\delta +\gamma }{P}_i^{(\alpha -\gamma ,\delta +\gamma )}\left(x\right).$$

This is the right equality (7). □

We will use the normalized shifted fractional Legendre polynomial for $\delta \ge 0$, defined as:

$${\varphi }_i^{\left(\delta \right)}\left(t\right):=\sqrt{{\displaystyle \frac{2i+1}{T}}}{t}^{\delta }{P}_i\left({\displaystyle \frac{2}{T}}t-1\right),i\in {\mathbb{N}}_0,t\in [0,T].$$

(8)

When $\delta =0$ occurs, the shifted fractional Legendre polynomial ${\varphi }_i^{\left(0\right)}\left(t\right)$ reduces to the shifted normalized Legendre polynomial.

The shifted fractional Legendre polynomial ${\varphi }_i^{\left(\delta \right)}\left(t\right)$ inherits the three-term recurrence relation of the classical Legendre polynomial, which reads as:

$$\begin{array}{cc}\hfill & {\varphi }_0^{\left(\delta \right)}\left(t\right)={\displaystyle \frac{1}{\sqrt{T}}}{t}^{\delta },{\varphi }_1^{\left(\delta \right)}\left(t\right)=\sqrt{{\displaystyle \frac{3}{T}}}{t}^{\delta }\left({\displaystyle \frac{2}{T}}t-1\right),\hfill \\ \hfill & {\varphi }_{i+1}^{\left(\delta \right)}\left(t\right)={\displaystyle \frac{\sqrt{(2i+1)(2i+3)}}{i+1}}\left({\displaystyle \frac{2}{T}}t-1\right){\varphi }_i^{\left(\delta \right)}\left(t\right)-{\displaystyle \frac{i}{i+1}}\sqrt{{\displaystyle \frac{2i+3}{2i-1}}}{\varphi }_{i-1}^{\left(\delta \right)}\left(t\right),i\ge 1.\hfill \end{array}$$

(9)

Additionally, the orthogonality property of the shifted fractional Legendre polynomial ${\varphi }_i^{\left(\delta \right)}\left(t\right)$ is given by

$${\int }_0^T{t}^{-2\delta }{\varphi }_i^{\left(\delta \right)}\left(t\right){\varphi }_j^{\left(\delta \right)}\left(t\right)\mathrm{d}t={\displaystyle \frac{\sqrt{(2i+1)(2j+1)}}{T}}{\int }_0^T{P}_i\left({\displaystyle \frac{2}{T}}t-1\right)P_j\left({\displaystyle \frac{2}{T}}t-1\right)\mathrm{d}t={\delta }_{ij}.$$

(10)

Hence, ${\left\{{\varphi }_i^{\left(\delta \right)}\left(t\right)\right\}}_{i=0}^{\infty }$ forms an orthonomal basis of $L_{t^{-2\delta }}^2(0,T)$.

3.2. Fractional Polynomial Frame and Approximation

Let $0={\delta }_1<{\delta }_2<\dots <{\delta }_K\in [0,+\infty )$ be K distinct nonnegative real numbers. We will use the set

$$\mathrm{\Phi }=\left\{{\varphi }_i^{\left({\delta }_k\right)}\left(t\right),i=0,\dots ,\infty ,k=1,\dots ,K\right\},$$

(11)

which is the collection of K families of the shifted fractional Legendre polynomials.

 Theorem 2.

Φ is a frame of $L^2(0,T)$.

Proof Noting that ${\delta }_1=0$. Since ${\left\{{\varphi }_i^{\left(0\right)}\right\}}_{i=0}^{\infty }$ is a complete orthogonal basis of $L^2(0,T)$, then for every $f\in L^2(0,T)$, we have

$${\parallel f\parallel }_{L^2}^2=\sum _{i=0}^{\infty }{\left|(f,{\varphi }_i^{\left(0\right)})\right|}^2.$$

Hence, we obtain a frame lower bound $A=1$.

Since that ${\varphi }_i^{\left({\delta }_k\right)}=t^{{\delta }_k}{\varphi }_i^{\left(0\right)}$, and $t^{{\delta }_k}$ is positive, increasing, bounded by $T^{{\delta }_k}$ on $[0,T]$, we have

$$|(f,{\varphi }_i^{\left({\delta }_k\right)}){|}^2=|(f,t^{{\delta }_k}{\varphi }_i^{\left(0\right)}){|}^2\le T^{2{\delta }_k}{\left|(f,{\varphi }_i^{\left(0\right)})\right|}^2.$$

Summing up all above inequalities, we have

$$\sum _{k=1}^K\sum _{i=0}^{\infty }|(f,{\varphi }_i^{\left({\delta }_k\right)}){|}^2\le B\sum _{i=0}^{\infty }|(f,{\varphi }_i^{\left(0\right)}){|}^2=B{\parallel f\parallel }_{L^2}^2$$

with $B={\sum }_{k=1}^K{T}^{2{\delta }_k}$. Now we give a frame upper bound. □

For discretization, we have to consider the truncated frame. Let $M_1,\dots ,M_K$ be K positive integer numbers. Denote

$$M=\sum _{k=1}^K\left(M_k+1\right).$$

The truncated frame is

$${\mathrm{\Phi }}_M=\left\{{\varphi }_i^{\left({\delta }_k\right)}\left(t\right),i=0,\dots ,M_k,k=1,\dots ,K\right\},$$

(12)

Hence, the dimension of space $\mathrm{span}\left\{{\mathrm{\Phi }}_M\right\}$ is M.

 Theorem 3.

${\mathrm{\Phi }}_M$ is a frame of $\mathrm{span}\left\{{\mathrm{\Phi }}_M\right\}$.

 Proof. 

This is a direct result of Proposition 1.1.2 in [24]. □

We sample $N+1$ points ${\left\{t_j\right\}}_{j=0}^N$ from $[0,T]$, referred as collocation points. Frame approximation of function $u\left(t\right)$ by ${\mathrm{\Phi }}_M$ leads to the least square problem

$$u_N\left(t\right)=\arg \underset{\varphi \in \mathrm{span}\left\{{\mathrm{\Phi }}_M\right\}}{min}{\displaystyle \frac{1}{N+1}}\sum _{j=0}^N{|u\left(t_j\right)-\varphi \left(t_j\right)|}^2.$$

Equivalently, the coefficients $\mathbf{c}=\left(c_i\right),i=1,\dots ,M$ of $u_N$ are the solution of the algebraic least-squares problem

$$\mathbf{c}=\arg \underset{\mathbf{x}\in {\mathbb{R}}^M}{min}\parallel \mathbf{A}\mathbf{x}-\mathbf{b}\parallel$$

with

$$\begin{array}{cc}\hfill \mathbf{A}& ={\left\{{\displaystyle \frac{1}{\sqrt{N+1}}}{\varphi }_m\left(t_j\right)\right\}}_{j=0,\dots ,N;m=1,\dots ,M}\in {\mathbb{R}}^{(N+1)\times M}\hfill \\ \hfill \mathbf{b}& ={\left\{{\displaystyle \frac{1}{\sqrt{N+1}}}u\left(t_j\right)\right\}}_{j=0,\dots ,N}\in {\mathbb{R}}^{N+1}.\hfill \end{array}$$

Theoretically, using frames to compute a best approximation requires solving an ill-conditioned linear system due to redundancy, which poses challenges for achieving high accuracy [19]. We address this by combining oversampling with truncated singular value decomposition (TSVD) regularization [20], which proves highly effective (see numerical results).

The oversampling technique requires $N>M$, with the ratio $N/M$ defined as the sampling ratio [20,21]. Thus, when the collocation condition satisfies $N+1>M$ (i.e., the number of collocation points exceeds the number of basis functions), the resulting system is overdetermined. We solve this system using the truncated singular value decomposition method, as described in Section 4.2.

The approximation errors of the frame ${\mathrm{\Phi }}_M$ in the general case remain unclear, mainly due to the complexity introduced by the singularity of the approximated function, leaving us unable to characterize the situation precisely. Fortunately, some results are known for special cases.

If $K=1$ and ${\delta }_1=0$, then $\mathrm{span}\left\{{\mathrm{\Phi }}_M\right\}={\mathbb{P}}_M$. In this case, the orthogonal projection ${\mathcal{P}}_M:L^2(-1,1)\to {\mathbb{P}}_M\left([-1,1]\right)$ converges spectrally fast, as stated in the following theorem.

 Theorem 4

([11]). If $u\in B_{0,0}^m(-1,1)$ and $0\le l\le m$, then

$${∥{\displaystyle \frac{{\mathrm{d}}^l}{\mathrm{d}{x}^l}}\left({\mathcal{P}}_{M}u-u\right)∥}_{{\omega }^{(l,l)}}\le C{M}^{l-m}{∥{\displaystyle \frac{{\mathrm{d}}^m}{\mathrm{d}{x}^m}}u∥}_{{\omega }^{(m,m)}},$$

where C is a generic positive constant independent of $M,u$, and the non-uniformly (or anisotropic) Jacobi-weighted Sobolev space is defined as

$$B_{\alpha ,\beta }^m(-1,1):=\left\{u:{\displaystyle \frac{{\mathrm{d}}^k}{\mathrm{d}{x}^k}}u\in L_{{\omega }^{(\alpha +k,\beta +k)}}^2(-1,1),0\le k\le m\right\},$$

with inner product

$${(u,v)}_{B_{\alpha ,\beta }^m}=\sum _{k=0}^m{\left({\displaystyle \frac{{\mathrm{d}}^k}{\mathrm{d}{x}^k}}u,{\displaystyle \frac{{\mathrm{d}}^k}{\mathrm{d}{x}^k}}v\right)}_{{\omega }^{(\alpha +k,\beta +k)}}.$$

If $K=2,{\delta }_1=0,{\delta }_2=\delta >0$, the frame ${\mathrm{\Phi }}_M$ is expected to approximate the function class as

$$S^{\delta }:=\left\{y\left(t\right)=g\left(t\right)+t^{\delta }h\left(t\right):g\left(t\right),h\left(t\right)\in C^{\infty }[0,T]\right\}.$$

 Theorem 5

(Proposition 5.9, in [19]). Let $y\left(t\right)\in S^{\delta }$, then there exists a coefficient vector $\mathbf{c}={(c_0^1,\dots ,c_{M_1}^1,c_0^2,\dots ,c_{M_2}^2)}^T\in {\mathbb{R}}^M$ such that ${\parallel \mathbf{c}\parallel }_2\le {\parallel g\parallel }_{L^2}+{\parallel h\parallel }_{L^2}$ and

$${∥y-\left[\sum _{j=0}^{M_1}{c}_j^1{\varphi }_j^{\left(0\right)}+\sum _{j=0}^{M_2}{c}_j^2{\varphi }_j^{\left(\delta \right)}\right]∥}_{L^{\infty }}\le {∥g-{\mathcal{P}}_{M_1}g∥}_{L^{\infty }}+T^{\delta }{∥h-{\mathcal{P}}_{M_2}h∥}_{L^{\infty }}.$$

Moreover, the numerical approximation $y_h={\mathcal{P}}_{N,M}^{ϵ}y$ (see Section 4.2) obtained by the frame collocation method with TSVD regularization satisfies

$${∥y-y_h∥}_{L^{\infty }}\le {∥g-{\mathcal{P}}_{M_1}g∥}_{L^{\infty }}+T^{\delta }{∥h-{\mathcal{P}}_{M_2}h∥}_{L^{\infty }}+\sqrt{ϵ}\left({\parallel g\parallel }_{L^2}+{\parallel h\parallel }_{L^2}\right),$$

where ϵ is the regularization parameter in TSVD.

4. Frame Collocation Method for Multi-Term FODEs

4.1. Fractional Polynomial Frame Collocation Scheme

Under consideration of initial value problem, we fix $t_0=0$ and ${\left\{t_j\right\}}_{j=1}^N\in (0,T]$ for collocation points, that can be chosen uniformly distributed, randomly distributed, or others like Gauss-type points. The frame collocation method for multi-term FODEs (2) reads as: to find $y_h\in \mathrm{span}\left\{{\mathrm{\Phi }}_M\right\}$ such that

$$\begin{array}{cc}\hfill & F\left(t,y_h\left(t\right),D^{{\gamma }_1}{y}_h\left(t\right),\dots ,D^{{\gamma }_S}{y}_h\left(t\right)\right)=0,t=t_j,j=1,\dots ,N\hfill \\ \hfill & y_h^{\left(k\right)}\left(t_0\right)=a_k,k=0,\dots ,\lceil {\gamma }_S\rceil -1.\hfill \end{array}$$

(13)

We rewrite the scheme (13) as a least square problem

$$\begin{array}{cc}\hfill y_h=\arg \underset{\varphi \in \mathrm{span}\left\{{\mathrm{\Phi }}_M\right\}}{min}(& {\displaystyle \frac{1}{N}}\sum _{j=1}^N{\left|F(t_j,\varphi \left(t_j\right),D^{{\gamma }_1}\varphi \left(t_j\right),\dots ,D^{{\gamma }_S}\varphi \left(t_j\right))\right|}^2\hfill \\ \hfill & +{\displaystyle \frac{1}{\lceil {\gamma }_S\rceil }}\sum _{k=0}^{\lceil {\gamma }_S\rceil -1}{|{\varphi }^{\left(k\right)}\left(t_0\right)-a_k|}^2).\hfill \end{array}$$

(14)

Let the approximate solution $y_h\left(t\right)$ be expressed as

$$y_h\left(t\right)=\sum _{k=1}^K\sum _{i=0}^{M_k}{c}_i^k{\varphi }_i^{\left({\delta }_k\right)}\left(t\right),t\in [0,T],$$

where $c_i^k$ are unknown coefficients to be determined. Denote

$$\begin{array}{cc}\hfill & \mathbf{t}={\left[t_1,\dots ,t_N\right]}^T,\hfill \\ \hfill & \mathbf{c}={\left[c_0^1,\dots ,c_{M_1}^1,\dots ,c_0^K,\dots ,c_{M_K}^K\right]}^T,\hfill \\ \hfill & \Psi \left(t\right)=\left[{\varphi }_0^{\left({\delta }_1\right)}\left(t\right),\dots ,{\varphi }_{M_1}^{\left({\delta }_1\right)}\left(t\right),\dots ,{\varphi }_0^{\left({\delta }_K\right)}\left(t\right),\dots ,{\varphi }_{M_K}^{\left({\delta }_K\right)}\left(t\right)\right].\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill & {\mathbf{y}}_h=y_h\left(\mathbf{t}\right),{\mathbf{y}}_h^{\left({\gamma }_k\right)}=y_h^{\left({\gamma }_k\right)}\left(\mathbf{t}\right),{\mathrm{\Phi }}_0={\left[\Psi \left(\mathbf{t}\right)\right]}_{N\times M},\hfill \\ \hfill & {\mathbf{D}}^{{\gamma }_k}={\left[D^{{\gamma }_k}\Psi \left(\mathbf{t}\right)\right]}_{N\times M},{\mathbf{D}}^{\left(k\right)}\Psi \left(0\right)={\left[{\Psi }^{\left(k\right)}\left(0\right)\right]}_{1\times M}(\mathbf{D}={\mathbf{D}}^{\left(1\right)}),\hfill \end{array}$$

then we have

$${\mathbf{y}}_h={\mathrm{\Phi }}_0\mathbf{c},{\mathbf{y}}_h^{\left({\gamma }_k\right)}={\mathbf{D}}^{{\gamma }_k}\mathbf{c},y_h^{\left(k\right)}\left(0\right)={\mathbf{D}}^{\left(k\right)}\Psi \left(0\right)\mathbf{c}.$$

For the initial condition, we denote

$$\mathbf{B}={\left[\Psi \left(0\right);\mathbf{D}\Psi \left(0\right);\dots ;{\mathbf{D}}^{(⌈{\gamma }_S⌉-1)}\Psi \left(0\right)\right]}_{\lceil {\gamma }_S\rceil \times M},$$

and $\mathbf{a}={\left[a_0,a_1,\dots ,a_{\lceil {\gamma }_S\rceil -1}\right]}^T$. Then initial condition gives $\mathbf{B}\mathbf{c}=\mathbf{a}.$

We introduce a nonlinear mapping G : ${\mathbb{R}}^M\to {\mathbb{R}}^{N+\lceil {\gamma }_S\rceil }$ as:

$$G\left(\mathbf{c}\right)=\left(\begin{array}{c}F\left(\mathbf{t},{\mathrm{\Phi }}_0\mathbf{c},{\mathbf{D}}^{{\gamma }_1}\mathbf{c},\dots ,{\mathbf{D}}^{{\gamma }_S}\mathbf{c}\right)\\ \mathbf{B}\mathbf{c}-\mathbf{a}\end{array}\right).$$

Now scheme (14) is equivalently the algebraic least-squares problem

$$\mathbf{c}=\arg \underset{\mathbf{x}\in {\mathbb{R}}^M}{min}\parallel G\left(\mathbf{x}\right)\parallel .$$

(15)

 Remark 1.

According to [20] , the number $N+1$ of the collocation points should be larger than M, empirically, ${\displaystyle \frac{N}{M}}\approx 2\mathrm{or}3$ (it is refered to as the sampling ratio).

We stress that the system (13) is overdetermined. In fact, there are $N+\lceil {\gamma }_S\rceil$ equations for M unknowns, and we always have $N+\lceil {\gamma }_S\rceil >N>M$, as stated in Remark 1. Hence, an efficient technique must be employed to solve the system (13), even in the linear case.

4.2. Truncated Singular Value Decomposition (TSVD)

We employ truncated singular value decomposition (TSVD) regularization [25], which solves the minimization problem

$$\underset{\mathbf{c}\in {\mathbb{R}}^M}{min}{\parallel \mathbf{A}\mathbf{c}-\mathbf{b}\parallel }_2^2,$$

(16)

where $\mathbf{A}\in {\mathbb{R}}^{N\times M}$ and $\mathbf{b}\in {\mathbb{R}}^N$ are given. Typically $N>M$, and the linear system $\mathbf{A}\mathbf{c}=\mathbf{b}$ is ill-posed.

Let the singular value decomposition of $\mathbf{A}$ be

$$\mathbf{A}=U\Sigma V^T,$$

where $U\in {\mathbb{R}}^{N\times N},V\in {\mathbb{R}}^{M\times M}$ are orthogonal matrices, and

$$\Sigma =diag\left({\sigma }_1,{\sigma }_2,\dots ,{\sigma }_M\right)\in {\mathbb{R}}^{N\times M}$$

contains the singular values ordered as

$${\sigma }_1\ge {\sigma }_2\ge \dots \ge {\sigma }_M.$$

Given a truncation threshold $ϵ>0$ (in our numerical experiments we set $ϵ={10}^{-14}$), we define the regularized pseudoinverse ${\Sigma }_{ϵ}^{†}$ whose diagonal entries are

$${\left({\Sigma }_{ϵ}^{†}\right)}_{ii}=\left\{\begin{array}{cc}{\displaystyle \frac{1}{{\sigma }_i}},\hfill & {\sigma }_i>ϵ\hfill \\ 0,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

The TSVD-regularized least-squares solution of (16) is then

$${\mathbf{c}}_{ϵ}={\mathbf{A}}^{†}\mathbf{b},\mathrm{where}{\mathbf{A}}^{†}=V{\Sigma }_{ϵ}^{†}{U}^T.$$

The corresponding approximate solution on $[0,T]$ is reconstructed as

$$y_h\left(t\right)=\sum _{j=1}^M{\left({\mathbf{c}}_{ϵ}\right)}_j{\Psi }_j\left(t\right).$$

Thus, the entire frame-based collocation method combined with TSVD regularization defines a direct mapping

$${\mathcal{P}}_{N,M}^{ϵ}:C\left([0,T]\right)⟶\mathrm{span}\left\{{\mathrm{\Phi }}_M\right\},y⟼y_h=\sum _{j=1}^M{\left({\mathbf{c}}_{ϵ}\right)}_j{\Psi }_j\left(t\right).$$

4.3. Evaluation of the Caputo Derivative of the Frame

A key step to implement the proposed frame collocation method is how to effectively evaluate the Caputo derivative of the frame basis function ${\varphi }_i^{\left(\delta \right)}\left(t\right)$ defined in (8) in Section 3.1.

Noting that

$${\varphi }_i^{\left(\delta \right)}\left(t\right)=\sqrt{{\displaystyle \frac{2i+1}{T}}}{\varphi }_i^{(\delta ,0,0)}\left(t\right),$$

we can evaluate $D^{\gamma }{\varphi }_i^{\left(\delta \right)}\left(t\right)$ by making use of Theorem 1.

For $0<\mu <1,i,k\in {\mathbb{N}}_0,\delta \ge 0$, let

$${\mathrm{I}}_{i,k}^{(\delta ,\mu )}={\mathrm{I}}_{i,k}^{(\delta ,\mu )}\left(t\right):={\displaystyle \frac{1}{\mathrm{\Gamma }\left(\mu \right)}}{\int }_0^t{\displaystyle \frac{{\varphi }_i^{(\delta ,k,k)}\left(\tau \right)}{{(t-\tau )}^{1-\mu }}}\mathrm{d}\tau ,$$

(17)

and

$${\widehat{\mathrm{I}}}_{i,k}^{(\delta ,\mu )}={\widehat{\mathrm{I}}}_{i,k}^{(\delta ,\mu )}\left(t\right):={\displaystyle \frac{1}{\mathrm{\Gamma }\left(\mu \right)}}{\int }_0^t{\displaystyle \frac{{\varphi }_i^{(\delta ,k,k)}\left(\tau \right)}{{(t-\tau )}^{-\mu }}}\mathrm{d}\tau .$$

(18)

The following recurrence relation is similar to Theorem 4.1 in [17].

 Theorem 6.

For $t\in [0,T]$ and $i\ge 1$, then ${\mathrm{I}}_{i,k}^{(\delta ,\mu )}$ and ${\widehat{\mathrm{I}}}_{i,k}^{(\delta ,\mu )}$ satisfy,

$${\mathrm{I}}_{i+1,k}^{(\delta ,\mu )}=A_i^{k,k}\left({\displaystyle \frac{2}{T}}t-1\right){\mathrm{I}}_{i,k}^{(\delta ,\mu )}-C_i^{k,k}{\mathrm{I}}_{i-1,k}^{(\delta ,\mu )}-{\displaystyle \frac{2}{T}}{A}_i^{k,k}{\widehat{\mathrm{I}}}_{i,k}^{(\delta ,\mu )},$$

(19)

and

$${\widehat{\mathrm{I}}}_{i+1,k}^{(\delta ,\mu )}={\tilde{A}}_i^{k,k}{\widehat{\mathrm{I}}}_{i,k}^{(\delta ,\mu )}+{\tilde{C}}_i^{k,k}{\widehat{\mathrm{I}}}_{i-1,k}^{(\delta ,\mu )}+t\left({\tilde{a}}_i^{k,k}{\mathrm{I}}_{i-1,k}^{(\delta ,\mu )}+{\tilde{c}}_i^{k,k}{\mathrm{I}}_{i+1,k}^{(\delta ,\mu )}\right)$$

(20)

with the starting terms

$$\begin{array}{cc}\hfill {\mathrm{I}}_{0,k}^{(\delta ,\mu )}=& {\displaystyle \frac{\mathrm{\Gamma }(\delta +1)}{\mathrm{\Gamma }(\delta +\mu +1)}}{t}^{\delta +\mu },{\mathrm{I}}_{1,k}^{(\delta ,\mu )}=(k+1)\left({\displaystyle \frac{2}{T}}{\displaystyle \frac{\delta +1}{\delta +\mu +1}}t-1\right){\mathrm{I}}_{0,k}^{(\delta ,\mu )},\hfill \\ \hfill {\widehat{\mathrm{I}}}_{0,k}^{(\delta ,\mu )}=& {\displaystyle \frac{\mu \mathrm{\Gamma }(\delta +1)}{\mathrm{\Gamma }(\delta +\mu +2)}}{t}^{\delta +\mu +1},{\widehat{\mathrm{I}}}_{1,k}^{(\delta ,\mu )}=(k+1)\left({\displaystyle \frac{2}{T}}{\displaystyle \frac{\delta +1}{\delta +\mu +2}}t-1\right){\widehat{\mathrm{I}}}_{0,k}^{(\delta ,\mu )}.\hfill \end{array}$$

where the coefficients are

$$\begin{array}{cc}\hfill {\tilde{A}}_i^{k,k}=& -{\displaystyle \frac{T^2{A}_i^{k,k}}{T^2-4(\delta +\mu +1)A_i^{k,k}{c}_i^{k,k}}},{\tilde{C}}_i^{k,k}=-{\displaystyle \frac{T^2{C}_i^{k,k}-4(\delta +\mu +1)A_i^{k,k}{a}_i^{k,k}}{T^2-4(\delta +\mu +1)A_i^{k,k}{c}_i^{k,k}}},\hfill \\ \hfill {\tilde{a}}_i^{k,k}=& {\displaystyle \frac{4\mu A_i^{k,k}{c}_i^{k,k}}{T^2-4(\delta +\mu +1)A_i^{k,k}{c}_i^{k,k}}},{\tilde{c}}_i^{k,k}={\displaystyle \frac{4\mu A_i^{k,k}{a}_i^{k,k}}{T^2-4(\delta +\mu +1)A_i^{k,k}{c}_i^{k,k}}},\hfill \end{array}$$

$A_i^{k,k},C_i^{k,k}$and $a_i^{k,k},c_i^{k,k}$ are given as in (A2) and (A5), respectively.

 Proof. 

From Theorem 1 with $B_i^{k,k}=0$, we have

$${\varphi }_{i+1}^{(\delta ,k,k)}\left(\tau \right)=A_i^{k,k}\left({\displaystyle \frac{2}{T}}t-1\right){\varphi }_i^{(\delta ,k,k)}\left(\tau \right)-C_i^{k,k}{\varphi }_{i-1}^{(\delta ,k,k)}\left(\tau \right)-{\displaystyle \frac{2}{T}}{A}_i^{k,k}(t-\tau ){\varphi }_i^{(\delta ,k,k)}\left(\tau \right).$$

Then the first relation (19) is obtained by the definition (17).

Next we will show the second relation (19). From the recurrence relation (6), we obtain

$${\widehat{\mathrm{I}}}_{i+1,k}^{(\delta ,\mu )}=-A_i^{k,k}{\widehat{\mathrm{I}}}_{i,k}^{(\delta ,\mu )}-C_i^{k,k}{\widehat{\mathrm{I}}}_{i-1,k}^{(\delta ,\mu )}+{\displaystyle \frac{2}{T}}{A}_i^{k,k}{\widehat{\mathrm{I}}}_{i,k}^{(\delta +1,\mu )}.$$

(21)

Let $h\left(x\right)=a_i^{k,k}{P}_{i-1}^{(k,k)}\left(x\right)+c_i^{k,k}{P}_{i+1}^{(k,k)}\left(x\right),$ then from (A4) we have $P_i^{(k,k)}\left({\displaystyle \frac{2}{T}}t-1\right)={\displaystyle \frac{2}{T}}{h}^{\prime }\left({\displaystyle \frac{2}{T}}t-1\right)$ and

$$\begin{array}{cc}\hfill {\widehat{\mathrm{I}}}_{i,k}^{(\delta +1,\mu )}=& {\displaystyle \frac{1}{\mathrm{\Gamma }\left(\mu \right)}}{\int }_0^t{\displaystyle \frac{{\tau }^{\delta +1}{P}_i^{(k,k)}\left(\frac{2}{T}\tau -1\right)}{{(t-\tau )}^{-\mu }}}\mathrm{d}\tau ={\displaystyle \frac{2}{T}}{\displaystyle \frac{1}{\mathrm{\Gamma }\left(\mu \right)}}{\int }_0^t{\displaystyle \frac{{\tau }^{\delta +1}{h}^{\prime }\left(\frac{2}{T}\tau -1\right)}{{(t-\tau )}^{-\mu }}}\mathrm{d}\tau \hfill \\ \hfill =& {\displaystyle \frac{2}{T}}{\displaystyle \frac{1}{\mathrm{\Gamma }\left(\mu \right)}}{\int }_0^t\left[\mu t{(t-\tau )}^{\mu -1}-(\delta +\mu +1){(t-\tau )}^{\mu }\right]{\tau }^{\delta }h\left({\displaystyle \frac{2}{T}}\tau -1\right)\mathrm{d}\tau \hfill \\ \hfill =& {\displaystyle \frac{2\mu t}{T}}{\displaystyle \frac{1}{\mathrm{\Gamma }\left(\mu \right)}}{\int }_0^t{\displaystyle \frac{{\tau }^{\delta }h\left(\frac{2}{T}\tau -1\right)}{{(t-\tau )}^{1-\mu }}}\mathrm{d}\tau +{\displaystyle \frac{2(\delta +\mu +1)}{T}}{\displaystyle \frac{1}{\mathrm{\Gamma }\left(\mu \right)}}{\int }_0^t{\displaystyle \frac{{\tau }^{\delta }h\left(\frac{2}{T}\tau -1\right)}{{(t-\tau )}^{-\mu }}}\mathrm{d}\tau \hfill \\ \hfill =& {\displaystyle \frac{2\mu t}{T}}\left[a_i^{k,k}{\mathrm{I}}_{i-1,k}^{(\delta ,\mu )}+c_i^{k,k}{\mathrm{I}}_{i+1,k}^{(\delta ,\mu )}\right]+{\displaystyle \frac{2(\delta +\mu +1)}{T}}\left[a_i^{k,k}{\widehat{\mathrm{I}}}_{i-1,k}^{(\delta ,\mu )}+c_i^{k,k}{\widehat{\mathrm{I}}}_{i+1,k}^{(\delta ,\mu )}\right].\hfill \end{array}$$

Hence, the second relation (19) is derived by merging the above equality and (21).

The starting terms can be derived by directly calculating. □

According to the equation (5) in Theorem 1,

$$\begin{array}{cc}\hfill D^{\gamma }{\varphi }_i^{\left(\delta \right)}\left(t\right)=& \sqrt{{\displaystyle \frac{2i+1}{T}}}{\displaystyle \frac{1}{\mathrm{\Gamma }(n-\gamma )}}\sum _{j=0}^n{d}_{j,i}^{\delta ,n}{\int }_0^t{\displaystyle \frac{{\varphi }_{i-n+j}^{(\delta -j,n-j,n-j)}\left(\tau \right)}{{(t-\tau )}^{\gamma -n+1}}}\mathrm{d}\tau \hfill \\ \hfill =& \sqrt{{\displaystyle \frac{2i+1}{T}}}\sum _{j=0}^n{d}_{j,i}^{\delta ,n}{\mathrm{I}}_{i-n+j,n-j}^{(\delta -j,n-\gamma )},\hfill \end{array}$$

(22)

where $n=\lceil \gamma \rceil$. Consequently, equation (22) together with Theorem 6 provides an efficient algorithm for evaluating $D^{\gamma }{\varphi }_i^{\left(\delta \right)}\left(t\right)$.

An interesting case of $\delta =0$ reduces to only one term in the summation of equation (22), that is,

$$D^{\gamma }{\varphi }_i^{\left(0\right)}\left(t\right)=\sqrt{{\displaystyle \frac{2i+1}{T}}}{\displaystyle \frac{(n+i)!}{T^{n}i!}}{\mathrm{I}}_{i-n,n}^{(0,n-\gamma )}.$$

 Remark 2.

An alternative way to evaluate $D^{\gamma }{\varphi }_i^{\left(\delta \right)}\left(t\right)$ is motivated by the integral relation (7). As shown in Appendix A (see (A7) and (A8)), there exists a transformation between different indices of the Jacobi polynomials, and a similar transformation holds for the shifted fractional Jacobi polynomial ${\varphi }_i^{(\delta ,\alpha ,\beta )}\left(t\right)$. Thus, we can transform the shifted fractional Jacobi polynomial in equation (22) into one with appropriate indices, so that the integral relation (7) can be utilized.

5. Numerical Experiments

To demonstrate the feasibility of the fractional polynomial frame collocation method (FPFCM), we implement the algorithm described in Section 4 using MATLAB to solve several FODEs, including the fractional Bagley–Torvik equations as well as linear and nonlinear multi-term FODEs.

To measure the accuracy, we employ the absolute maximum error defined as

$$E=\underset{0\le i\le 1000}{max}\left\{|y\left(t_i\right)-y_h\left(t_i\right)|\right\},t_i={\displaystyle \frac{iT}{1000}},i=0,1,\dots ,1000,$$

where $y\left(t\right)$ and $y_h\left(t\right)$ are the exact and numerical solutions of the problem, respectively.

Since no further information is available to guide the selection of collocation points, our numerical experiments employ uniformly distributed collocation points, that is,

$$t_j=j{\displaystyle \frac{T}{N}},j=0,\dots ,N.$$

In fact, other selection strategies (such as random or Gauss-type points) do not introduce additional difficulties in the fractional polynomial frame collocation method.

5.1. The Fractional Bagley–Torvik Equation

The fractional Bagley–Torvik equation is presented to describe the motion of a rigid plate immersed in a Newtonian fluid of infinite extent [27,28]. Since the Bagley–Torvik equation plays a vital role in many problems of applied science and engineering, a large number of researchers show interest in its solution such as in [23,29,30,31,32].

The fractional Bagley-Torvik equation reads

$$A{y}^{\prime \prime }\left(t\right)+B{D}^{\gamma }y\left(t\right)+Cy\left(t\right)=f\left(t\right),0<t\le T.$$

(23)

where $A,B,C$ are constant numbers and $f\left(t\right)$ is the source term. Equation (23) may equipped with initial value condition as

$$y\left(0\right)=a_0,y^{\prime }\left(0\right)=a_1.$$

(24)

 Example 1.

Consider the initial value problem of equation (23) with $\gamma =3/2,A=B=C=1,T=1.$ For the exact solution $y\left(t\right)$, we consider the following three cases:

C11.

$y\left(t\right)=8t^{3.5}+{\displaystyle \frac{1}{8}}{t}^4-t^3$.

C12.

$y\left(t\right)=t^2{E}_{0.5,3}\left({\displaystyle \frac{10t^{0.5}}{13}}\right)$.

C13.

$y\left(t\right)=e^t+t^{2.5}sin\left(t\right)$.

Here in C12, $E_{\alpha ,\beta }(·)$ is the two-parameter Mittag–Leffler function, defined as

$$E_{\alpha ,\beta }\left(z\right)=\sum _{k=0}^{\infty }{\displaystyle \frac{z^k}{\mathrm{\Gamma }(\alpha k+\beta )}}.$$

The initial values $a_0,a_1$ are determined by the exact solution $y\left(t\right)$. It is clear that the initial value conditions are $a_0=a_1=0$ for the first two cases, and $a_0=a_1=1$ for the C13.

The source terms $f\left(t\right)$ for the three cases are specified as

$$\begin{array}{cc}\hfill C11:f\left(t\right)=& \left(70-{\displaystyle \frac{8}{\sqrt{\pi }}}\right)t^{1.5}+\left(1.5+{\displaystyle \frac{105\sqrt{\pi }}{4}}\right)t^2+{\displaystyle \frac{8}{5\sqrt{\pi }}}{t}^{2.5}-6t-t^3+8t^{3.5}+{\displaystyle \frac{1}{8}}{t}^4,\hfill \\ \hfill C12:f\left(t\right)=& E_{0.5,1}\left({\displaystyle \frac{10t^{0.5}}{13}}\right)+t^{0.5}{E}_{0.5,1.5}\left({\displaystyle \frac{10t^{0.5}}{13}}\right)+t^2{E}_{0.5,3}\left({\displaystyle \frac{10t^{0.5}}{13}}\right),\hfill \\ \hfill C13:f\left(t\right)=& 2\sum _{n=0}^{\infty }{\displaystyle \frac{t^n}{n!}}+\sum _{k=0}^{\infty }{\displaystyle \frac{{(-1)}^k}{(2k+1)!}}{\displaystyle \frac{\mathrm{\Gamma }(2k+4.5)}{\mathrm{\Gamma }(2k+3)}}{t}^{2k+2}+\sum _{k=0}^{\infty }{\displaystyle \frac{t^{2k+0.5}}{\mathrm{\Gamma }(2k+1.5)}}\hfill \\ \hfill & +\sum _{k=0}^{\infty }\left[{\displaystyle \frac{1}{\mathrm{\Gamma }(2k+2.5)}}+{\displaystyle \frac{{(-1)}^k(10k+8.75)}{(2k+1)!}}\right]t^{2k+1.5}.\hfill \end{array}$$

In this test, the sampling ratio is taken as 2, which means $N=2M=2{\sum }_{k=1}^K(M_k+1)$.

For the case C11, the frame is specified by taking $K=2$. Since the exact solution takes form $y\left(t\right)=g\left(t\right)+t^{3.5}h\left(t\right)$ with $g\left(t\right),h\left(t\right)\in C^{\infty }[0,T]$, we choose ${\delta }_1=0,{\delta }_2=1.5\&2.5\&3.5$. In this case, $g\left(t\right)\in {\mathbb{P}}_4$, it is enough to take $M_1=4$. We test for different $M_2=0,1,2,3$. For the regularization parameter $ϵ$ in TSVD, we take $ϵ={10}^{-13}$ in the case.

The maximum absolute errors of the case C11 for $M_2=0,1,2,3$ are presented in Table 1.

It is shown from Table 1 that the FPFCM achieves almost machine precision only using a few basis functions from the frame.

For the case C12, the problem is solved by direct piecewise polynomial collocation methods [32]. Since the exact solution is low regularity, it is challenging to achieve high-order convergence rate by most of the existing numerical methods. The FPFCM is employed to solve this problem with $K=2,M_1=M_2$ and ${\delta }_1=0,{\delta }_2=1.5\&2.5$. For the regularization parameter $ϵ$ in TSVD, we take $ϵ={10}^{-13}$ in the case. We list the maximum absolute errors in Table 2, together with the errors obtained by direct piecewise polynomial collocation methods [32] for comparison.

From Table 2, it is clear that the proposed method achieves “spectral accuracy”, whilst the direct piecewise polynomial collocation method in [32] has convergence rate about $1.5\sim 3.0$.

For the case C13, the FPFCM is employed to solve this problem with $K=2,M_1=M_2$ and ${\delta }_1=0,{\delta }_2=1.5\&2.5\&3.5$. For the regularization parameter $ϵ$ in TSVD, we take $ϵ={10}^{-14}$ in the case.

We list the maximum absolute errors in Table 3, it is clearly shown that the convergence rate is “spectral accuracy” for this case.

5.2. The Multi-Term Linear FODE

Consider a general multi-term linear FODE with constant or variable coefficients as

$$\sum _{k=1}^S{A}_k\left(t\right)D^{{\gamma }_k}y\left(t\right)=f\left(t\right),0<t\le T.$$

(25)

Equation (25) may be equipped with some boundary or initial conditions for well-posedness, which will be specified in the following numerical examples. The existence, uniqueness, convergence of the solution for multi-order fractional differential equation are investigated, such as in [33,34].

 Example 2.

Consider equation (25) of the initial value problem with $S=4,T=1$, orders, coefficents, initial value conditions and exact solution as the following two cases:

C21.

Orders and coefficents:

$$({\gamma }_1,{\gamma }_2,{\gamma }_3,{\gamma }_4)=(3,2.5,1.5,0),(A_1,A_2,A_3,A_4)=(4/3,-15/4,-7/3,-12/5).$$

The initial value conditions: $y\left(0\right)=y^{\prime }\left(0\right)=y^{\left(2\right)}\left(0\right)=0.$

The exact solution: $y\left(t\right)=8t^{3.5}+{\displaystyle \frac{1}{8}}{t}^4-t^3$.

C22.

Orders and coefficents:

$$({\gamma }_1,{\gamma }_2,{\gamma }_3,{\gamma }_4)=(2,0,1,1/2),(A_1,A_2,A_3,A_4)=(1,1,-2,1).$$

The initial value conditions: $y\left(0\right)=y^{\prime }\left(0\right)=0.$

The exact solution: $y\left(t\right)=t^7-t^2$.

For the case C21, Wang and Liang solved the same problem by direct piecewise polynomial collocation method in [32]. We employ the FPFCM to solve this problem with $K=2,M_1=4$ and ${\delta }_1=0,{\delta }_2=2.5\&3.5$. We take regularization parameter $ϵ={10}^{-14}$ in TSVD. Table 4 lists the maximum absolute errors together with the errors obtained by direct piecewise polynomial collocation method [32] for comparison.

From Table 4, our method achieves almost machine precision only using a few basis functions, whilst the direct piecewise polynomial collocation method in [32] has finite convergence rate about $2.8$. It is shown our method is superior.

For the case C22, Han et. al. solved the same problem by a method of the shifted Chebyshev polynomials fractional differential operator matrix with error correction in [35]. Since the exact solution belongs to ${\mathbb{P}}_7$, it is enough to take $K=1,{\delta }_1=0$ and $M=M_1=7(N=2M)$ in the FPFCM. Table 5 lists the maximum absolute errors together with the errors obtained in [35] for comparison.

As shown in Table 5, our method achieves nearly the same accuracy as the method in [35]. In fact, our method converges exponentially for the solution, be it smooth or involving multiple singularities.

 Example 3.

Consider the boundary value problem in equation (25) with $S=3$, $T=1$, $0<\gamma <1$. The orders and coefficients for the two cases are:

C31.

Orders and coefficents: $({\gamma }_1,{\gamma }_2,{\gamma }_3)=(2,\gamma ,0),\mathrm{and}(A_1,A_2,A_3)=(1,t^2-1,-sin\left(2\pi t\right)).$

The boundary value conditions: $y\left(0\right)=2,y\left(1\right)+y^{\prime }\left(1\right)=1$.

The exact solution: $y\left(t\right)=t^{\gamma }(t^2-t^3)+2$.

C32.

Orders and coefficents: $({\gamma }_1,{\gamma }_2,{\gamma }_3)=(\gamma +1,\gamma ,0),\mathrm{and}(A_1,A_2,A_3)=(1,t,t^2).$

The boundary value conditions: $y\left(0\right)=0,y\left(1\right)=1$.

The exact solution: $y\left(t\right)=t^{3-\gamma /3}+t^{1.3}lnt$.

The source terms $f\left(t\right)$ of equation (25) for the two cases are specified as

$$\begin{array}{cc}\hfill & C31:f\left(t\right)=(2+\gamma )(1+\gamma )t^{\gamma }-(3+\gamma )(2+\gamma )t^{1+\gamma }+(t^2-1)\left[{\displaystyle \frac{\mathrm{\Gamma }(3+\gamma )}{2}}{t}^2-{\displaystyle \frac{\mathrm{\Gamma }(4+\gamma )}{6}}{t}^3\right]\hfill \\ \hfill & -sin\left(2\pi t\right)(t^{2+\gamma }-t^{3+\gamma }+2).\hfill \\ \hfill & C32:f\left(t\right)={\displaystyle \frac{\mathrm{\Gamma }(4-\gamma /3)}{\mathrm{\Gamma }(3-4\gamma /3)}}{t}^{2-4\gamma /3}+{\displaystyle \frac{\mathrm{\Gamma }\left(2.3\right)}{\mathrm{\Gamma }(1.3-\gamma )}}{t}^{0.3-\gamma }\left(lnt+\psi \left(2.3\right)-\psi (1.3-\gamma )\right)\hfill \\ \hfill & +{\displaystyle \frac{\mathrm{\Gamma }(4-\gamma /3)}{\mathrm{\Gamma }(4-4\gamma /3)}}{t}^{4-4\gamma /3}+{\displaystyle \frac{\mathrm{\Gamma }\left(2.3\right)}{\mathrm{\Gamma }(2.3-\gamma )}}{t}^{2.3-\gamma }\left(lnt+\psi \left(2.3\right)-\psi (2.3-\gamma )\right)\hfill \\ \hfill & +t^{5-\gamma /3}+t^{3.3}lnt,\hfill \end{array}$$

where $\psi \left(z\right):={\displaystyle \frac{{\mathrm{\Gamma }}^{\prime }\left(z\right)}{\mathrm{\Gamma }\left(z\right)}}$ is the digamma function.

For the case C31, Mary et. al. solved the same problem by a method of the finite difference method combined with L1 Scheme (FDM-L1) in [36]. We employ the FPFCM with $K=2,N=2M,{\delta }_1=0,{\delta }_2=\gamma \&(\gamma +1)\&(\gamma +2),M_1=0$ to solve it. We take regularization parameter $ϵ={10}^{-14}$ in TSVD. Table 6 and Table 7 list the maximum absolute errors together with the errors obtained by the finite difference method combined with L1 Scheme (FDM-L1) in [36] for comparison.

From Table 6Table 7, it is shown that the proposed method achieves superior accuracy with a few basis functions if appropricate parameters ${\left\{{\delta }_k\right\}}_{k=1}^K$ are choosed.

For the case C32, the exact solution exhibits distinct types of singularities at the left endpoint $t=0$. We employ the FPFCM with $K=2,N=2M,{\delta }_1=\left(1-{\displaystyle \frac{\gamma }{3}}\right)\&\left(2-{\displaystyle \frac{\gamma }{3}}\right)\&\left(1-{\displaystyle \frac{\gamma }{3}}\right),{\delta }_2=1.3,M_1=M_2$ to solve it. We take regularization parameter $ϵ={10}^{-14}$ in TSVD. Table 8, Table 9 and Table 10 list the maximum absolute errors.

For the special case, the proposed method is observed to have a finite convergence rate, indicating that the constructed frame fails to capture the singularity of $lnt$. Despite this limitation, the method provides a relatively satisfactory computational precision with a small number of basis functions, as shown in Table 8, Table 9 and Table 10.

5.3. The Nonlinear Two-Term FODEs

In the subsection, we consider the following nonlinear multi-term FODEs.

$$D^{\gamma }y\left(t\right)=F(t,y\left(t\right),D^{{\gamma }_1}y\left(t\right)),0<t\le T.$$

(26)

We limit $0\le {\gamma }_1<\gamma \le 1$, so only one initial value condition is required, that is, $y\left(0\right)=a_0$.

 Example 4.

Consider equation (26) with $T=1$ and

C41.

$F(t,y,D^{{\gamma }_1}y)=f\left(t\right)-y^{{\displaystyle \frac{3}{2}}}$ with

$$f\left(t\right)={\displaystyle \frac{40320}{\mathrm{\Gamma }(9-\gamma )}}{t}^{8-\gamma }-{\displaystyle \frac{3\mathrm{\Gamma }(5+\gamma /2)}{\mathrm{\Gamma }(5-\gamma /2)}}{t}^{4-\gamma /2}+{\left({\displaystyle \frac{3}{2}}{t}^{\gamma /2}-t^4\right)}^3+{\displaystyle \frac{9}{4}}\mathrm{\Gamma }(\gamma +1).$$

The initial value is $a_0=0$. The exact solution [37,38] is

$$y\left(t\right)=t^8-3t^{4+\gamma /2}+{\displaystyle \frac{9}{4}}{t}^{\gamma }.$$

C42.

$F(t,y,D^{{\gamma }_1}y)=f\left(t\right)+{\displaystyle \frac{{\gamma }_1{t}^{{\gamma }_1}}{5}}{D}^{{\gamma }_1}y+{\displaystyle \frac{t}{10}}{y}^3$ with

$$f\left(t\right)=\left\{\begin{array}{cc}{\displaystyle \frac{\mathrm{\Gamma }(2\gamma +1)}{\mathrm{\Gamma }(\gamma +1)}}{t}^{\gamma }-{\displaystyle \frac{t}{10}}{\left(t^{2\gamma }+1\right)}^3,\hfill & {\gamma }_1=0\hfill \\ {\displaystyle \frac{\mathrm{\Gamma }(2\gamma +1)}{\mathrm{\Gamma }(\gamma +1)}}{t}^{\gamma }+{\displaystyle \frac{\mathrm{\Gamma }(2{\gamma }_1+1)}{\mathrm{\Gamma }(2{\gamma }_1-\gamma +1)}}{t}^{2{\gamma }_1-\gamma }-{\displaystyle \frac{{\gamma }_1}{5}}{\displaystyle \frac{\mathrm{\Gamma }(2\gamma +1)}{\mathrm{\Gamma }(2\gamma -{\gamma }_1+1)}}{t}^{2\gamma }\hfill & \\ -{\displaystyle \frac{{\gamma }_1}{5}}{\displaystyle \frac{\mathrm{\Gamma }(2{\gamma }_1+1)}{\mathrm{\Gamma }({\gamma }_1+1)}}{t}^{2{\gamma }_1}-{\displaystyle \frac{t}{10}}{\left(t^{2\gamma }+t^{2{\gamma }_1}\right)}^3,\hfill & {\gamma }_1>0.\hfill \end{array}\right.$$

The initial value is $a_0=0$ if ${\gamma }_1>0$ and $a_0=1$ if ${\gamma }_1=0$. The exact solution is

$$y\left(t\right)=t^{2\gamma }+t^{2{\gamma }_1}.$$

For C41, the exact solution contains terms $t^8$, $t^{4+\gamma /2}$, and $t^{\gamma }$, which exhibit different singularity exponents at $t=0$. We employ FPFCM with $K=3$, $N=2M$, ${\delta }_1=0$, ${\delta }_2=\gamma$, ${\delta }_3=\left(1+{\displaystyle \frac{\gamma }{2}}\right)\&\left(2+{\displaystyle \frac{\gamma }{2}}\right)\&\left(3+{\displaystyle \frac{\gamma }{2}}\right)$, with $M_1=8,M_2=0$. We give the absolute errors in Table 11, Table 12 and Table 13.

From Table 11, Table 12 and Table 13, it is shown that the proposed method achieves high precision with a small number of basis functions. In particular, when ${\delta }_3=1+{\displaystyle \frac{\gamma }{2}}$ and $M=14$ ($M_3=3$), the absolute errors reach the order of ${10}^{-13}$ for all $\gamma$ values. Similarly, for ${\delta }_3=2+{\displaystyle \frac{\gamma }{2}}$ and $3+{\displaystyle \frac{\gamma }{2}}$, errors drop to ${10}^{-14}$–${10}^{-15}$ as M increases. These results demonstrate that choosing fractional basis parameters matching the singularity exponents of the exact solution leads to superior performance.

For C42, the exact solution is $y\left(t\right)=t^{2\gamma }+t^{2{\gamma }_1}$, which has two distinct singularity exponents at $t=0$. We employ the FPFCM with $K=2$, ${\delta }_1=2\gamma$, ${\delta }_2=2{\gamma }_1$, and $M_1=0$. Table 14 lists the maximum absolute errors for $\gamma =0.9$ and various ${\gamma }_1$.

It is shown that the proposed method achieves excellent accuracy with a few basis functions.

 Remark 3.

Some interesting facts about the convergence behavior of the proposed method are confirmed by numerical tests, which we summarize as follows:

• The exponential convergence of the proposed method with the right parameters is demonstrated for both smooth solutions and low-regularity solutions of the type $t^{\delta }$. Nevertheless, the convergence rate is degraded when the solution involves other kinds of singularities, such as $lnt$ (see case C32 of Example 3).
• The choice of the parameters of the proposed method is crucial to the convergence rate. In practice, the optimal parameters are those that match the singularity of the solution.
• The data of the FODEs, including coefficients and source terms, have little effect on the convergence behavior of the proposed method.
• If little information is available about the regularity of the solution, one can take $K=2$, ${\delta }_1=0$, and ${\delta }_2=0.5$; then the FODE problem can be solved effectively in most cases.
• The optimal sampling ratio remains an open problem.

6. Conclusions

Due to their enhanced capability to accurately describe systems characterized by multiple relaxation times and multi-rate anomalous diffusion in heterogeneous media, multi-term fractional differential equations (FODEs) have found extensive applications in multi-physics coupling and multi-scale processes. Nevertheless, their numerical solution presents considerable challenges, primarily arising from the inherent lack of smoothness in the solutions of multi-term FODEs.

In the present study, a high-precision numerical method is successfully developed for solving multi-term fractional ordinary differential equations whose solutions exhibit multiple singularities, through the integration of frame theory with the collocation method. The proposed approach employs shifted fractional Legendre polynomials to construct a hybrid basis function frame, which flexibly accommodates the singular structure of the solution. Consequently, the rapid convergence characteristic of spectral methods is retained, while the accuracy limitations of conventional approaches in handling singular solutions are effectively overcome. To mitigate the severe ill-conditioning frequently encountered in discretized systems, truncated singular value decomposition (TSVD) regularization is introduced, thereby substantially enhancing the numerical stability of the algorithm.

Numerical experiments confirm that the proposed method achieves superior convergence performance in the numerical solution of multi-term fractional ODEs, exhibiting superconvergence particularly when the singularity exponent coincides with the basis function parameter $\delta$. The method is not only applicable to problems exhibiting pronounced singular structures but also efficient for smooth solutions, thereby demonstrating broad applicability and practical utility. It thus constitutes a reliable and efficient numerical tool for high-precision simulation of multi-term fractional differential equations with complex singular structures.

It is also noted that the proposed method can be readily extended to solve fractional differential equations of variable-order of a more general form [26]. Future work will be directed toward the application of the proposed method to fractional differential equations in higher dimensions.