A Unified Gegenbauer Framework for Arbitrary-Order Riemann-Liouville Fractional Integrals with Tunable Endpoint Regularity

1. Introduction

Fractional calculus provides a powerful framework for modeling systems exhibiting memory, nonlocal interactions, and power-law dynamics [1,2]. Central to this theory is the RLFI of order $\alpha >0$, defined by

$${\mathfrak{I}}_t^{\alpha }f={\displaystyle \frac{1}{\Gamma \left(\alpha \right)}}{\mathcal{I}}_t^{\left(\tau \right)}\left({(t-\tau )}^{\alpha -1}f\right).$$

(1)

While fractional orders in ${\mathbf{\Omega }}_1^{\circ }$ have been extensively studied due to their prevalence in physical applications [42,43,44], emerging problems in applied electromagnetics [45], structural dynamics [46], intelligent infrastructure, electrical engineering, nonlinear dynamics, biomechanics, control systems, thermal and fluid sciences, and biomedical signal analysis increasingly require integration orders beyond this interval [47,48].

1.1. Challenges in Arbitrary-Order RLFI Computation

Despite its closed-form definition (1), the RLFI rarely admits an analytic expression for general f and arbitrary $\alpha >0$. Even for smooth or elementary f, the nonlocal kernel ${(t-\tau )}^{\alpha -1}$ couples all past values of f, defying symbolic integration except in special cases (e.g., monomials or exponentials with integer $\alpha$). Moreover, for non-integer $\alpha >1$, the kernel belongs only to $C^{\lfloor \alpha -1\rfloor }\left({\mathbf{\Omega }}_t\right)$, which precludes the use of classical higher-order Newton–Cotes or Gauss rules without transformation. In practical modeling–such as anomalous transport, viscoelasticity, and control–where f is often given numerically or as part of a coupled system, fast, stable, and repeated evaluation of ${\mathfrak{I}}_t^{\alpha }f$ becomes essential.

1.2. Overview of Existing Numerical Approaches

This need has spurred decades of research into numerical quadrature methods. Foundational work in polynomial-based and operational matrix methods for fractional calculus includes early contributions from orthogonal polynomial systems [52] as well as more recent extensions and refinements [50,51,54]. More recent developments include operational matrices using Charlier polynomials [24], fractional-order Chelyshkov functions [41,49], and SG-based GBFA schemes [4] that exhibit super-exponential convergence for $0<\alpha <1$—a result we later extend and refine in Section A.3, showing that full super-convergence of the classical GBFA formulation (including the quadrature component) occurs precisely when $1/\alpha \in {\mathbb{Z}}^{+}$. Concurrently, spline interpolation schemes (linear to quintic) offered flexible local approximations with rigorous error estimates [25,26,53]. Composite and trapezoidal-type quadrature rules for fractional integrals, frequently combined with Romberg extrapolation, have been developed and analyzed in recent years [27,28]. In parallel, a variety of projection and operational matrix techniques–employing block-pulse functions, hybrid polynomial systems, and Galerkin formulations–have been developed over the past decade to improve computational efficiency [30,31,32,33,34,41]. Advanced quadrature techniques for weakly singular fractional integrals continue to be an active research area. These include variable transformations such as the sigmoidal Sidi–Laurie maps [56,57], Sinc-based collocation [37,38,39], and modern adaptations in Sobolev settings [35,36], all aimed at improving convergence by regularizing endpoint behavior. Spectral and wavelet bases, such as Mittag-Leffler wavelets, have further enhanced adaptivity and accuracy [40,55].

1.3. Limitations of Current Methods

While these existing alternatives provide valuable tools for specific regimes of $\alpha$, they generally lack one or more of the following attributes essential for high-performance, repeated RLFI evaluation: (i) fixed-node, non-adaptive structure compatible with global solvers; (ii) precomputable integration matrices enabling low online cost; (iii) generally offers rapid (often spectral) convergence for smooth f across all $\alpha >0$; and (iv) tunable endpoint regularity to control quadrature accuracy. The GBFA framework developed in this work–extending the earlier approach of Elgindy [4]–uniquely satisfies all four properties simultaneously, as demonstrated later in the paper.

The GBFA method [4] employs the transformation

$$\tau =t(1-y^{1/\alpha }),y\in {\mathbf{\Omega }}_1,$$

(2)

to obtain the exact identity

(3)

Although the kernel ${(t-\tau )}^{\alpha -1}$ of (1) is nonsingular on ${\mathbf{\Omega }}_t$ for $\alpha \ge 1$, the composed integrand $f\left(t(1-y^{1/\alpha })\right)$ is only Hölder-continuous of order $1/\alpha$ at $y=0$, with all derivatives blowing up as $y\to 0^{+}$ (as we demonstrate later via Faà di Bruno’s formula). This loss of endpoint smoothness degrades Gauss-type quadrature convergence despite the original kernel being $C^{\lfloor \alpha -1\rfloor }$. Adaptive quadrature can compensate but is computationally expensive, non-reproducible across evaluations, and incompatible with fixed-grid solvers (e.g., collocation, operational matrices, time-stepping schemes). Hence, a unified, non-adaptive, highly accurate framework operating on fixed nodes, yielding precomputable FSGIMs, maintaining high algebraic, spectral, or even super-exponential convergence rates whenever the endpoint regularity permits, and integrating seamlessly into global approximation architectures remains compelling.

1.4. The GBFA Method and Its Original Scope

Within this landscape, Elgindy [4] introduced the GBFA method, which achieves super-exponential convergence for the RLFI when $1/\alpha \in {\mathbb{Z}}^{+}$ via SG interpolation and precomputable FSGIMs, while exhibiting algebraic convergence for other fractional orders $0<\alpha <1$; cf. Section A.3. The method hinges on transformation (2), yielding identity (3), and exhibits: (i) super-exponential convergence for analytic f and $1/\alpha \in {\mathbb{Z}}^{+}$; (ii) precomputable FSGIMs; and (iii) tunable SG parameters $\lambda$, ${\lambda }_q$ for optimal stability and accuracy. Although the original presentation restricted $\alpha$ to ${\mathbf{\Omega }}_1^{\circ }$, transformation (2) is mathematically valid ${\forall }_a\alpha >0$ as we show later in Section 3. This paper demonstrates that the entire GBFA framework–including SG interpolation, FSGIM construction, and error analysis–extends naturally to arbitrary $\alpha >0$ without algorithmic modification. Additionally, we introduce and analyze a family of generalized $\varphi \left(\alpha \right)$-transformations that restore $C^r$ endpoint regularity for $\alpha >1$ (for any $r\in {\mathbb{Z}}^{+}$), thereby enabling high-order quadrature convergence across all $\alpha >0$. We further note that $C^r$ endpoint regularity is likewise attainable for $0<\alpha <1$ whenever $1/\alpha \notin {\mathbb{Z}}^{+}{\forall }_{sl}{n}_q$, as demonstrated later. This insight significantly broadens GBFA’s applicability to problems requiring $\alpha \ge 1$ or even integer-order repeated integration.

Existing alternatives for arbitrary-order RLFIs–such as adaptive quadrature [5], operational matrix techniques [6,7], and wavelet-based methods [8,9]–often trade accuracy for computational cost or implementation complexity. In contrast, the generalized GBFA preserves spectral accuracy, enables efficient FSGIM reuse, and, as shown herein, attains near machine-precision across $\alpha >0$ with minimal parameter tuning.

1.5. Contributions of This Work

The main contributions of this work include:

(i)

A proof that the GBFA method of [4], originally developed for $0<\alpha <1$, requires no algorithmic modification to extend to arbitrary $\alpha >0$;

(ii)

Extension of the exact transformed identity (3)—originally derived in [4] for $0<\alpha <1$—and all associated theoretical results (interpolation, quadrature, FSGIM construction, and error analysis) to any arbitrary $\alpha >0$;

(iii)

Numerical validation of excellent convergence rates and robust performance over a wide range of $\alpha >0$;

(iv)

Introduction of a novel $\varphi \left(\alpha \right)$-generalized transformation for $\alpha >1$ that enforces tunable $C^r$ endpoint regularity (for any prescribed $r\ge 0$) in the transformed integrand, dramatically accelerating quadrature convergence and enabling near machine-precision accuracy with modest quadrature orders through appropriate choice of r and spectral parameters $\lambda ,{\lambda }_q$; moreover, we show that $C^r$ endpoint regularity is likewise attainable for $0<\alpha <1$ whenever $1/\alpha \notin {\mathbb{Z}}^{+}{\forall }_{sl}{n}_q$, as demonstrated later.

(v)

Provision of a unified framework which delivers an exact algebraic substitution that precisely cancels the kernel ${(t-\tau )}^{\alpha -1}$ for $0<\alpha <1$ (classical case) while restoring controllable endpoint smoothness for $\alpha \ge 1$ (generalized case), and likewise for $0<\alpha <1$ whenever $1/\alpha \notin {\mathbb{Z}}^{+}$; this is distinct from approximate sigmoidal regularization techniques such as Sidi–Laurie transformations, which are primarily designed for weakly singular kernels when $0<\alpha <1$.

(vi)

A refined error analysis that establishes a complete convergence hierarchy for the ITE across all smoothness classes of f (globally analytic, Gevrey-regular with index $1<\mu <2$, general $C^{\infty }$, and finitely smooth $C^k$), clarifying the precise regularity conditions required for super-exponential, log-stretched-exponential, qualitative, or algebraic decay rates—extending and sharpening the bounds in [4], which focused primarily on analytic data.

(vii)

A sharp characterization of the endpoint regularity of the transformed integrand under Elgindy’s classical map, establishing that full super-exponential convergence of the GBFA method (including exact quadrature for polynomial data) occurs iff $1/\alpha \in {\mathbb{Z}}^{+}$; otherwise, the integrand belongs only to a Hölder space $C^{K,\theta }\left({\mathbf{\Omega }}_1\right)$ with $K=\lfloor 1/\alpha \rfloor$ and $\theta =1/\alpha -K$, limiting quadrature to algebraic asymptotic rates. In particular, for $0<\alpha <1$ with $1/\alpha \notin {\mathbb{Z}}^{+}$, one may still enforce an error of order $\mathcal{O}\left(n^{-r}\right)$ for any $r>0$${\forall }_{sl}{n}_q$. Crucially, we clarify via Remark A9 that for real-analytic f, the pre-asymptotic quadrature error often decays dramatically faster than this worst-case rate due to exponentially decaying coefficients in the fractional power series expansion of the transformed integrand, reconciling observed near-exponential convergence in practical computations even when $1/\alpha \notin {\mathbb{Z}}^{+}$.

1.6. Organization of the Paper

The remainder of the paper is organized as follows: Section 2 reviews the GBFA framework of [4]. Section 3 proves this framework is valid ${\forall }_a\alpha >0$, analyzes endpoint regularity, and presents convergence results for analytic, Gevrey, $C^{\infty }$, and $C^k$ data, along with quadrature error and $\lambda$-dependence. Section 4 introduces the $\varphi \left(\alpha \right)$-generalized transformation, and constructs weighted FSGIMs. Section 5 analyzes complexity. Section 6 derives the smoothness condition. Section 7 validates both transformations for $\alpha >0$, compares GBFA against some methods and solvers available in the literature. Section 8 summarizes implications, guidelines, trade-offs, and future directions.

2. Background: GBFA Method

This section briefly summarizes the GBFA method as presented in [4], highlighting the components that naturally generalize to arbitrary $\alpha >0$.

2.1. SG Polynomials, the RLFI Transformation, and FSGIMs

As detailed in [4, Section 2], the SG polynomials ${\widehat{G}}_n^{\left(\lambda \right)}\left(t\right)$, defined for $t\in {\mathbf{\Omega }}_1$ and $\lambda >-1/2$, form an orthogonal basis on ${\mathbf{\Omega }}_1$ with respect to the weight function ${(t-t^2)}^{\lambda -1/2}$. For $f\in L^2\left({\mathbf{\Omega }}_1\right)$ and SGG nodes ${\widehat{\mathbb{G}}}_n^{\lambda }=\left\{{\widehat{t}}_{n,0:n}^{\lambda }\right\}$, the GBFA interpolant is:

(4)

where ${\mathcal{L}}_{0:n}^{\lambda }\left(t\right)$ are the Lagrange basis polynomials in SG form. The core innovation in [4] is Elgindy’s transformation (2), which, for $0<\alpha <1$, yields (3). The GBFA approximation then becomes:

(5)

The approximation formula (5) shows that the RLFI reduces, after the transformation (2), to evaluating integrals of the form

(6)

We refer to (6) as the SG–transformed integrals. They arise naturally from substituting the SG interpolant of f into the transformed RLFI and represent the contribution of each SG basis function to the fractional integral.

Because the mapping $y\mapsto t(1-y^{1/\alpha })$ is a $C^{\infty }$-diffeomorphism on the open interval ${\mathbf{\Omega }}_1^{\circ }$ for all $\alpha >0$, the integrands in (6) are smooth functions of y on ${\mathbf{\Omega }}_1^{\circ }$, even when the original RLFI kernel ${(t-\tau )}^{\alpha -1}$ is weakly singular for $0<\alpha <1$. This interior regularity supports high-order SG quadrature convergence. However, for $\alpha >1$, the mapping is only Hölder-continuous of order $1/\alpha$ at the endpoint $y=0$, with all derivatives blowing up as $y\to 0^{+}$ (as shown later in Section 3.6). This endpoint behavior reduces the global regularity of the transformed integrand on ${\mathbf{\Omega }}_1$, affecting the quadrature error component but not the interpolation error (which depends only on the smoothness of f in the original variable t).

In practice, the integrals (6) are approximated using high-order Gauss quadratures such as the $(n_q,{\lambda }_q)$-GBFA quadrature rule [4]. This introduces a quadrature truncation error that we will analyze later in Section 3.5 for arbitrary $\alpha >0$. The explicit definition (6) is therefore essential for understanding both the construction of the FSGIMs and the total error structure of the GBFA method.

As developed in [4, Section 3], the $\alpha$th-order FSGIM ${}^E\mathbf{Q}_n^{\alpha }$ enables efficient computation:

$${\mathfrak{I}}_{{\mathbf{z}}_M}^{\alpha }f\approx {}^E\mathbf{Q}_n^{\alpha }{f}_{0:n},$$

(7)

for evaluation points ${\mathbf{z}}_M\in {\mathbf{\Omega }}_1^{M+1}$. These matrices are precomputable and depend only on $\alpha$, n, $\lambda$, and quadrature parameters; see [4, Eqs. (18) and (20)].

Motivated by the structural properties reviewed above, we now turn to one of the central theoretical contributions of this work.

3. Generalization to Arbitrary Fractional Orders

With the essential machinery now in place, we demonstrate in this section how the GBFA framework naturally extends to arbitrary $\alpha >0$ without modification.

3.1. Mathematical Validity of Elgindy’s Transformation

The transformation (2) requires careful examination of the exponent $1/\alpha$. For $\alpha >0$, the mapping (2) is bijective on ${\mathbf{\Omega }}_t^{\circ }\leftrightarrow {\mathbf{\Omega }}_1^{\circ }$ with Jacobian (A46) (see A.2). Substituting into the RLFI definition yields:

(8)

which directly institutes (3). Notice here the precise cancellation of the exponents of y: the term $y^{(\alpha -1)/\alpha }$ arising from the kernel ${(t-\tau )}^{\alpha -1}$ is perfectly neutralized by the term $y^{(1-\alpha )/\alpha }$ appearing in the Jacobian $d\tau /dy$. This identity holds for all $\alpha >0$; thus, the core GBFA framework—including the approximation formula (5), the FSGIM construction from [4, Section 2], and the evaluation form (7)—extends directly to arbitrary positive orders without modification. Since the transformation (2) is exact for all $\alpha >0$, the algorithm for generating ${}^E\mathbf{Q}_n^{\alpha }$ requires no $\alpha$-dependent changes. However, for $\alpha >1$, the endpoint regularity of the transformed integrand is reduced (becoming only Hölder-continuous of order $1/\alpha$ at $y=0$), which affects the quadrature truncation error component but leaves the interpolation error convergence rate intact as we prove later in Section 3.6.

With the transformation rigorously justified, we can now analyze how the error structure behaves under arbitrary fractional orders.

3.2. Error Analysis

The error analysis developed in [4, Section 4] carries over to the present framework without any change in the proof structure or asymptotic rates; only the constants in the bounds acquire an implicit dependence on $\alpha$ through the composed integrand $f\left(t(1-y^{1/\alpha })\right)$. The corresponding convergence results are given below in Theorems 1–5, covering analytic, Gevrey-regular, general $C^{\infty }$-smooth functions, and finitely smooth data with algebraic convergence rates. In particular, Theorems 1 and 2 quantify super-exponential and log-stretched-exponential decay for analytic and Gevrey-regular inputs, Theorems 3 and 4 ensure qualitative convergence for $C^{\infty }$ data, and Theorem 5 establishes algebraic convergence of order $\mathcal{O}\left(n^{-k}\right)$ for functions in $C^k\left({\mathbf{\Omega }}_1\right)$.

Theorem 1 

(Super-Exponential Convergence for Globally Analytic f). Let $\alpha >0$ and $t\in {\mathbf{\Omega }}_1$. Suppose that f is real-analytic on ${\mathbf{\Omega }}_1$ and that there exist constants $C_f,B_f>0$ such that

$$\parallel f^{\left(k\right)}{\parallel }_{L^{\infty }\left({\mathbf{\Omega }}_1\right)}\le C_f{B}_f^k,\forall k\in {\mathbb{Z}}_0^{+}.$$

(9)

Let ${\mathfrak{I}}_t^{\alpha }f$ be approximated via the GBFA method using SG interpolation of degree n. Then $\exists \xi =\xi \left(t\right)\in {\mathbf{\Omega }}_1^{\circ }$ such that the ITE admits the exact representation

(10)

where ${\widehat{K}}_n^{\lambda }$ is given by [4, Eq. (24)]. Moreover, the GBFA approximation converges to ${\mathfrak{I}}_t^{\alpha }f$ at a super-exponential rate in n; more precisely, there exists a constant $c>0$ independent of n such that

$$\left|{\mathfrak{I}}_t^{\alpha }f-{\mathfrak{I}}_t^{\alpha }{I}_{n}f\right|=O\left(e^{-cnlnn}\right),n\to \infty .$$

(11)

Proof. 

The exact error representation (10) follows directly from [4, Theorem 1], which remains valid for any $\alpha >0$ since the transformation (2) is a $C^{\infty }$-diffeomorphism between ${\mathbf{\Omega }}_1^{\circ }$ and ${\mathbf{\Omega }}_t^{\circ }$${\forall }_a$$\alpha >0$; see Corollary A2. To obtain an explicit asymptotic bound, we invoke the general estimate derived in [4, Theorem 2], which yields

$$\left|{\mathfrak{I}}_t^{\alpha }f-{\mathfrak{I}}_t^{\alpha }{I}_{n}f\right|\stackrel{\sim }{<}{\mathcal{A}}_{n+1}{\vartheta }_{\alpha ,\lambda }{\left(\frac{e}{4}\right)}^n{n}^{-{\displaystyle \frac{3}{2}}-n+\alpha }{{\rm Y}}_{{\sigma }^{\lambda }}\left(n\right),$$

(12)

where

$${\mathcal{A}}_{n+1}:={\parallel f^{(n+1)}\parallel }_{L^{\infty }\left({\mathbf{\Omega }}_1\right)},$$

$$\begin{array}{c}{\vartheta }_{\alpha ,\lambda }=\frac{1}{4\pi }{e}^{\alpha +1}{\alpha }^{-{\displaystyle \frac{1}{2}}-\alpha }{\left(1+2\lambda \right)}^{-{\displaystyle \frac{1}{2}}-2\lambda }\left[1+{\displaystyle \frac{1}{1620{\left(1+\lambda \right)}^5}}\right]\times \\ {\left[{\displaystyle \frac{\left(1+\lambda \right)csch\left(\frac{1}{1+2\lambda }\right)}{1+2\lambda }}\right]}^{{\displaystyle \frac{1}{2}}+\lambda }{\left[\alpha sinh\left({\displaystyle \frac{1}{\alpha }}\right)\right]}^{-\alpha /2}{\left[\left(1+\lambda \right)sinh\left({\displaystyle \frac{1}{1+\lambda }}\right)\right]}^{{\displaystyle \frac{1+\lambda }{2}}},\end{array}$$

(13)

and

$${{\rm Y}}_{{\sigma }^{\lambda }}\left(n\right)=\left\{\begin{array}{cc}1,\hfill & \lambda \ge 0,\hfill \\ {\sigma }^{\lambda }{n}^{-\lambda },\hfill & -{\displaystyle \frac{1}{2}}<\lambda <0,\hfill \end{array}\right.$$

with ${\sigma }^{\lambda }>1$. By assumption (9), we have the growth bound

$${\mathcal{A}}_{n+1}\le C_f{B}_f^{n+1}.$$

(14)

Substituting (14) into (12) gives

$$\left|{\mathfrak{I}}_t^{\alpha }f-{\mathfrak{I}}_t^{\alpha }{I}_{n}f\right|\stackrel{\sim }{<}{C}_f{B}_f^{n+1}{\left(\frac{e}{4}\right)}^n{n}^{-{\displaystyle \frac{3}{2}}-n+\lambda }{{\rm Y}}_{{\sigma }^{\lambda }}\left(n\right).$$

(15)

Collecting the n-dependent factors, we may rewrite (15) as

$$\left|{\mathfrak{I}}_t^{\alpha }f-{\mathfrak{I}}_t^{\alpha }{I}_{n}f\right|\le C_1{\left({\displaystyle \frac{e{B}_f}{4}}\right)}^n{n}^{-{\displaystyle \frac{3}{2}}-n+\lambda }{{\rm Y}}_{{\sigma }^{\lambda }}\left(n\right),$$

(16)

for some constant $C_1>0$ independent of n. Since ${{\rm Y}}_{{\sigma }^{\lambda }}\left(n\right)$ is at most algebraic in n, there exists $p\in \mathbb{R}$ and $C_2>0$ such that

$${{\rm Y}}_{{\sigma }^{\lambda }}\left(n\right)\le C_2{n}^p,n\ge 1.$$

Thus,

$$\begin{array}{cc}\hfill \left|{\mathfrak{I}}_t^{\alpha }f-{\mathfrak{I}}_t^{\alpha }{I}_{n}f\right|& \le C_3{\left({\displaystyle \frac{e{B}_f}{4}}\right)}^n{n}^{-n+q},\hfill \\ \hfill & =C_{3}exp\left(nln(e{B}_f/4)-nlnn+qlnn\right).\hfill \end{array}$$

(17)

${\forall }_s$$C_3>0$ and $q\in \mathbb{R}$ independent of n. As $n\to \infty$, the dominant term in the exponent is $-nlnn$. Hence, there exists $c>0$ and $N\in {\mathbb{Z}}^{+}$ such that

$$nln(e{B}_f/4)-nlnn+qlnn\le -cnlnn,\forall n\ge N,$$

which implies

$${\left({\displaystyle \frac{e{B}_f}{4}}\right)}^n{n}^{-n+q}\le e^{-cnlnn},n\ge N.$$

This proves that the GBFA error decays super-exponentially in n under the growth condition (9).    □

The following theorem is tailored for smooth, but non-globally analytic function f; see, for example, [13] for a detailed discussion of Gevrey-regularity. By log-stretched-exponential convergence here we mean stretched-exponential with logarithmic enhancement; i.e., faster than stretched-exponential but slower than exponential.

Theorem 2 

(Log-Stretched-Exponential Convergence for Gevrey Regular f). Let $\alpha >0$ and $t\in {\mathbf{\Omega }}_1$. Suppose that f is $C^{\infty }$-smooth on ${\mathbf{\Omega }}_1$ and that there exist constants $C_f,B_f>0$ and $1<\mu <2$ such that

$$\parallel f^{\left(k\right)}{\parallel }_{L^{\infty }\left({\mathbf{\Omega }}_1\right)}\le C_f{B}_f^k{(k!)}^{\mu -1},\forall k\in {\mathbb{Z}}_0^{+}.$$

(18)

Let ${\mathfrak{I}}_t^{\alpha }f$ be approximated via the GBFA method using SG interpolation of degree n. Then $\exists \xi =\xi \left(t\right)\in {\mathbf{\Omega }}_1^{\circ }$ such that the ITE admits the exact representation (10). Moreover, there exists a constant $c>0$ independent of n such that

$$\left|{\mathfrak{I}}_t^{\alpha }f-{\mathfrak{I}}_t^{\alpha }{I}_{n}f\right|=O\left(e^{-c{n}^{1/\mu }lnn}\right),n\to \infty .$$

(19)

In particular, the GBFA approximation converges faster than any algebraic rate in n for Gevrey-regular f of order μ.

Proof. 

As in Theorem 1, the exact error representation (10) follows from [4, Theorem 1] and the diffeomorphism property of the transformation (2). Consider again the general GBFA error estimate (12). Under the Gevrey growth condition (18), we have

$${\mathcal{A}}_{n+1}\le C_f{B}_f^{n+1}{\left((n+1)!\right)}^{\mu -1}.$$

(20)

Using Stirling’s approximation $(n+1)!\sim \sqrt{2\pi (n+1)}{\left({\displaystyle \frac{n+1}{e}}\right)}^{n+1}$, we obtain

$${\left((n+1)!\right)}^{\mu -1}\sim {\left(2\pi \right)}^{{\displaystyle \frac{\mu -1}{2}}}{(n+1)}^{{\displaystyle \frac{\mu -1}{2}}}{\left({\displaystyle \frac{n+1}{e}}\right)}^{(\mu -1)(n+1)}.$$

(21)

Substituting (20) and (21) into (12), and absorbing ${\left(2\pi \right)}^{(\mu -1)/2}$ and the factor $B_f$ into a generic constant, we obtain ${\forall }_{sl}$n:

$$\begin{array}{cc}\hfill \left|{\mathfrak{I}}_t^{\alpha }f-{\mathfrak{I}}_t^{\alpha }{I}_{n}f\right|& \le C_1{\left({\displaystyle \frac{e{B}_f}{4}}\right)}^n{\left({\displaystyle \frac{n}{e}}\right)}^{(\mu -1)n}{n}^{\lambda +{\displaystyle \frac{3}{2}}\mu -n-3}{{\rm Y}}_{{\sigma }^{\lambda }}\left(n\right),\hfill \end{array}$$

(22)

${\forall }_s{C}_1>0$ independent of n. Since ${{\rm Y}}_{{\sigma }^{\lambda }}\left(n\right)$ grows at most algebraically in n, there exist $q\in \mathbb{R}$ and $C_2>0$ such that

$${{\rm Y}}_{{\sigma }^{\lambda }}\left(n\right)\le C_2{n}^q,n\ge 1.$$

Thus, from (22) we obtain

$$\begin{array}{cc}\hfill \left|{\mathfrak{I}}_t^{\alpha }f-{\mathfrak{I}}_t^{\alpha }{I}_{n}f\right|& \le C_3{\left({\displaystyle \frac{e{B}_f}{4}}\right)}^n{\left({\displaystyle \frac{n}{e}}\right)}^{(\mu -1)n}{n}^{q^{\prime }-n},\hfill \end{array}$$

(23)

${\forall }_s{C}_3>0$ and $q^{\prime }\in \mathbb{R}$ independent of n. Writing the right-hand side of (23) in exponential form,

$$\begin{array}{cc}\hfill & {\left({\displaystyle \frac{e{B}_f}{4}}\right)}^n{\left({\displaystyle \frac{n}{e}}\right)}^{(\mu -1)n}{n}^{q^{\prime }-n}\hfill \\ \hfill & =exp\left((\mu -2)nlnn+n\left(ln(e{B}_f/4)-(\mu -1)\right)+q^{\prime }lnn\right).\hfill \end{array}$$

For $1<\mu <2$, the coefficient $\mu -2$ is negative, hence the leading term is

$$(\mu -2)nlnn=-(2-\mu )nlnn,$$

which dominates the lower-order contributions $\mathcal{O}\left(n\right)$ and $q^{\prime }lnn$ as $n\to \infty$. Therefore, there exist constants $c_1>0$, $C_4>0$, and $N\in {\mathbb{Z}}^{+}$ such that

$${\left({\displaystyle \frac{e{B}_f}{4}}\right)}^n{\left({\displaystyle \frac{n}{e}}\right)}^{(\mu -1)n}{n}^{q^{\prime }-n}\le C_{4}exp\left(-c_{1}nlnn\right),n\ge N.$$

Since $nlnn\ge n^{1/\mu }lnn$${\forall }_s$$n\ge 2$ and $1<\mu <2$, it follows that

$$exp\left(-c_{1}nlnn\right)\le exp\left(-c{n}^{1/\mu }lnn\right),{\forall }_{sl}n,$$

${\forall }_{s}c>0$, which proves (19).    □

Theorem 3 

(Qualitative Convergence for Higher-Order Gevrey-Regular f). Let $\alpha >0$ and $t\in {\mathbf{\Omega }}_1$. Suppose that f is $C^{\infty }$-smooth on ${\mathbf{\Omega }}_1$ and that there exist constants $C_f,B_f>0$ and $\mu \ge 2$ such that

$$\parallel f^{\left(k\right)}{\parallel }_{L^{\infty }\left({\mathbf{\Omega }}_1\right)}\le C_f{B}_f^k{(k!)}^{\mu -1},\forall k\in {\mathbb{Z}}_0^{+}.$$

(24)

Let ${\mathfrak{I}}_t^{\alpha }f$ be approximated via the GBFA method using SG interpolation of degree n. Then the GBFA approximation converges pointwise to the exact RLFI, i.e.,

$$\underset{n\to \infty }{lim}\left({\mathfrak{I}}_t^{\alpha }f-{\mathfrak{I}}_t^{\alpha }{I}_{n}f\right)=0.$$

(25)

Moreover, the convergence is uniform in t on any compact subset $K\subset {\mathbf{\Omega }}_1^{\circ }$.

Proof. 

The Gevrey condition (24) with $\mu \ge 2$ implies that $f\in C^{\infty }\left({\mathbf{\Omega }}_1\right)$. For such smooth functions, spectral interpolation at Gauss-type nodes yields uniform convergence. Specifically, for SG interpolation at the SGG nodes ${\widehat{\mathbb{G}}}_n^{\lambda }$, there exists a sequence ${\left\{{\eta }_n\right\}}_{n=1}^{\infty }$ with ${\eta }_n\to 0$ such that

$$\parallel f-I_n{f\parallel }_{L^{\infty }\left({\mathbf{\Omega }}_1\right)}\le {\eta }_n,\forall n\in {\mathbb{Z}}^{+}.$$

(26)

This follows from standard results in polynomial approximation theory; see, e.g., the convergence properties of interpolation in orthogonal polynomial bases at Gauss quadrature nodes [10,11,12,58].

Next, the RLFI ${\mathfrak{I}}_t^{\alpha }(·)$ defines a bounded linear operator on $C\left({\mathbf{\Omega }}_1\right)$. For any $g\in C\left({\mathbf{\Omega }}_1\right)$,

$$\left|{\mathfrak{I}}_t^{\alpha }g\right|\le {\displaystyle \frac{1}{\Gamma \left(\alpha \right)}}{\int }_0^t{(t-\tau )}^{\alpha -1}\left|g\left(\tau \right)\right|d\tau \le {\displaystyle \frac{t^{\alpha }}{\Gamma (\alpha +1)}}{\parallel g\parallel }_{L^{\infty }\left({\mathbf{\Omega }}_1\right)},t\in {\mathbf{\Omega }}_1.$$

(27)

Applying this bound with $g=f-I_{n}f$ and using (26), we obtain for every $t\in {\mathbf{\Omega }}_1$:

$$\begin{array}{cc}\hfill \left|{\mathfrak{I}}_t^{\alpha }f-{\mathfrak{I}}_t^{\alpha }{I}_{n}f\right|& =\left|{\mathfrak{I}}_t^{\alpha }(f-I_{n}f)\right|\le {\displaystyle \frac{t^{\alpha }}{\Gamma (\alpha +1)}}{\parallel f-I_{n}f\parallel }_{L^{\infty }\left({\mathbf{\Omega }}_1\right)}\hfill \\ \hfill & \le {\displaystyle \frac{t^{\alpha }}{\Gamma (\alpha +1)}}{\eta }_n.\hfill \end{array}$$

(28)

Consequently, for any fixed $t\in {\mathbf{\Omega }}_1^{\circ }$, since ${\eta }_n\to 0$, we have the pointwise convergence (25).

For the uniform convergence on compact subsets, let $K\subset {\mathbf{\Omega }}_1^{\circ }$ be compact. Then ${sup}_{t\in K}{t}^{\alpha }<\infty$, and we obtain

$$\underset{t\in K}{sup}\left|{\mathfrak{I}}_t^{\alpha }f-{\mathfrak{I}}_t^{\alpha }{I}_{n}f\right|\le {\displaystyle \frac{1}{\Gamma (\alpha +1)}}\underset{t\in K}{sup}{t}^{\alpha }{\eta }_n⟶0\mathrm{as}n\to \infty ,$$

which establishes uniform convergence on K.    □

Theorem 4 

(Qualitative Convergence for Smooth Data). Let $\alpha >0$ and $t\in {\mathbf{\Omega }}_1$. Suppose that $f\in C^{\infty }\left({\mathbf{\Omega }}_1\right)$ and that, for each $n\in {\mathbb{Z}}^{+}$, $I_{n}f$ denotes the SG interpolant of f at the SGG nodes ${\widehat{\mathbb{G}}}_n^{\lambda }\subset {\mathbf{\Omega }}_1$ as in (4). Assume that the SG interpolation operator satisfies

$$\underset{n\to \infty }{lim}{\parallel f-I_{n}f\parallel }_{L^{\infty }\left({\mathbf{\Omega }}_1\right)}=0.$$

(29)

Let ${\mathfrak{I}}_t^{\alpha }f$ be approximated via the GBFA method using SG interpolation of degree n. Then (25) is satisfied, and the convergence is uniform in t on any compact subset $K\subset {\mathbf{\Omega }}_1^{\circ }$.

Proof. 

By assumption $f\in C^{\infty }\left({\mathbf{\Omega }}_1\right)$ and (29) holds, i.e., the SG interpolants $I_{n}f$ at the SGG nodes converge uniformly to f on ${\mathbf{\Omega }}_1$. This is precisely the interpolation property used in the proof of Theorem 3, where the Gevrey condition (24) was invoked only to guarantee smoothness and the uniform convergence of SG interpolation. The remaining step in the proof of Theorem 3 relies on two ingredients: (i) the bounded linearity of the RLFI ${\mathfrak{I}}_t^{\alpha }(·)$ on $C\left({\mathbf{\Omega }}_1\right)$, and (ii) the uniform convergence of $I_{n}f$ to f in $L^{\infty }\left({\mathbf{\Omega }}_1\right)$. These ingredients are entirely independent of the Gevrey growth condition and require only (29). Consequently, the same argument shows that (25) holds true for each fixed $t\in {\mathbf{\Omega }}_1^{\circ }$, and that the convergence is uniform in t on any compact subset $K\subset {\mathbf{\Omega }}_1^{\circ }$. This establishes the proof.    □

Remark 1 

(Minimal Versus Structured Regularity). Theorem 4 identifies minimal smooth-data assumptions for qualitative convergence of the GBFA method–namely: (i) f is $C^{\infty }$-smooth on ${\mathbf{\Omega }}_1$, and (ii) the SG interpolants converge uniformly to f as in (29). In fact, the proof of convergence fundamentally relies only on the interpolation property (29) and the boundedness of the RLFI operator, and does not preclude convergence in weaker function spaces. By contrast, Theorem 3 concerns the more restrictive class of Gevrey-regular functions with index $\mu \ge 2$. Although this theorem also yields qualitative convergence, it plays a different conceptual role: it completes the Gevrey hierarchy established in Theorems 2–3 and shows that, as the Gevrey index crosses the threshold $\mu =2$, the GBFA method transitions from quantitative (stretched-exponential) convergence to purely qualitative convergence. Thus, Theorem 3 serves to characterize the behavior of the GBFA method across the full spectrum of Gevrey regularity, whereas Theorem 4 isolates the weakest smooth-data conditions under which convergence still holds in this framework.

3.3. Dependence of the GBFA ITE on the Parameter $\lambda$

The dependence of the GBFA ITE on the SG parameter $\lambda$ arises through several components of the general estimate (12). The dominant $\lambda$-dependence is contained in the constant ${\vartheta }_{\alpha ,\lambda }$, which, as shown in [4], is unimodal in $\lambda$, attaining its maximum near ${\lambda }^{*}\approx -0.1351$ and then decaying exponentially as $\lambda$ increases. This exponential decay is the primary mechanism governing the behavior of the error with respect to $\lambda$. In contrast, the remaining $\lambda$-dependent factors in (12), namely the algebraic term $n^{-{\displaystyle \frac{3}{2}}-n+\lambda }$ and, for $-1/2<\lambda <0$, the auxiliary factor ${{\rm Y}}_{{\sigma }^{\lambda }}\left(n\right)={\sigma }^{\lambda }{n}^{-\lambda }$, contribute only algebraic growth or decay in n. For fixed n, these algebraic contributions are dominated by the exponential decay of ${\vartheta }_{\alpha ,\lambda }$, and therefore do not alter the qualitative trend that the GBFA ITE decreases as $\lambda$ increases beyond ${\lambda }^{*}$.

For moderate fractional orders $\alpha \in {\mathbf{\Omega }}_{1,2}^{\circ }$, this unimodal structure of ${\vartheta }_{\alpha ,\lambda }$ explains why increasing $\lambda$ improves the ITE. For larger $\alpha$, the dominant factor ${\alpha }^{-{\displaystyle \frac{1}{2}}-\alpha }$ in (13) decays super-exponentially, suppressing the $\lambda$-dependence of all terms in (12). This accounts for the numerical observation in Section 7 that, for $\alpha >2$, the GBFA error becomes nearly invariant across all tested values of $\lambda$. Although the analysis in [4] focused on $0<\alpha <1$, the structure of the $\lambda$-dependent terms in (12) remains unchanged for all $\alpha >0$, and the qualitative behavior of the ITE with respect to $\lambda$ therefore extends naturally to arbitrary fractional orders.

Figure 1 extends the numerical study of [4, Fig. 3] from the classical regime $0<\alpha <1$ to fractional orders $\alpha >1$. The results exhibit the same qualitative structure observed previously: for moderate values of $\alpha$ (e.g., $\alpha \in {\mathbf{\Omega }}_{1,2}^{\circ }$), the GBFA error retains a clear dependence on the SG parameter $\lambda$, with smaller errors obtained for increasing $\lambda >{\lambda }^{*}\approx -0.1351$, in agreement with the unimodal behavior of ${\vartheta }_{\alpha ,\lambda }$ described earlier in this section. As $\alpha$ increases further, the dominant factor ${\alpha }^{-{\displaystyle \frac{1}{2}}-\alpha }$ in the constant ${\vartheta }_{\alpha ,\lambda }$ decays super-exponentially, effectively suppressing the $\lambda$-dependence of the entire bound (12). This explains the near-coincidence of all curves for increasing $\alpha$ values in the classical transformation, where the total GBFA error is dominated by the quadrature component. In this regime, the super-exponential decay of the factor ${\alpha }^{-{\displaystyle \frac{1}{2}}-\alpha }$ in ${\vartheta }_{\alpha ,\lambda }$ suppresses the underlying $\lambda$-dependence, so the method appears only weakly sensitive to $\lambda$. Thus, Figure 1 demonstrates that the parameter-sensitivity patterns reported in [4] for $0<\alpha <1$ persist for $\alpha >1$, but gradually diminish as $\alpha$ grows, in full agreement with the theoretical predictions.

Remark 2. 

Although the theoretical bound (12) indicates that the dependence on λ should diminish as α increases, numerical experiments using the $\varphi \left(\alpha \right)$-generalized transformation introduced in Section 4 exhibit a more nuanced trend. For this generalized mapping, moderate values of λ (typically $\lambda \in {\mathbf{\Omega }}_1$) continue to deliver excellent accuracy for large α. In contrast, extreme values near the theoretical lower limit $\lambda =-0.5^{+}$ show noticeably degraded performance as α grows, effectively widening the practical λ-sensitivity range. These observations suggest that, in applications of the generalized transformation, λ should be selected from ${\mathbf{\Omega }}_1$ to avoid the blow-up region near $-0.5^{+}$, while ${\lambda }_q$ should be fixed at $0.5$ to fully exploit the $C^r$ endpoint regularization, as demonstrated in Section 7.3.

3.4. Finite Smoothness Case

The previous results focused on globally analytic functions, Gevrey-regular data, and $C^{\infty }$-smooth functions. In many practical applications, however, the input f may possess only finite smoothness, say $f\in C^k\left({\mathbf{\Omega }}_1\right)$ for some finite $k\in {\mathbb{Z}}^{+}$. In this case, one expects algebraic—rather than (super-)exponential—convergence, in line with classical results from spectral and pseudospectral approximation theory; see, for example, [11]. The following theorem formalizes this behavior for the GBFA approximation of the RLFI.

Theorem 5 

(Algebraic convergence for $C^k$ data). Let $\alpha >0$ and $t\in {\mathbf{\Omega }}_1$. Suppose that $f\in C^k\left({\mathbf{\Omega }}_1\right)$ for some integer $k\ge 1$, and that its kth derivative satisfies

$$\parallel f^{\left(k\right)}{\parallel }_{L^{\infty }\left({\mathbf{\Omega }}_1\right)}\le C_k,$$

(30)

for some constant $C_k>0$. For each $n\in {\mathbb{Z}}^{+}$, let $I_{n}f$ denote the SG interpolant of f at the SGG nodes ${\widehat{\mathbb{G}}}_n^{\lambda }\subset {\mathbf{\Omega }}_1$ as in (4). Then there exist constants $C_{\mathrm{app}},C_{\alpha ,k}>0$, independent of n, such that

$$\parallel f-I_n{f\parallel }_{L^{\infty }\left({\mathbf{\Omega }}_1\right)}\le C_{\mathrm{app}}{n}^{-k},\forall n\in {\mathbb{Z}}^{+},$$

(31)

and the GBFA approximation of the RLFI satisfies the algebraic error bound

$$\left|{\mathfrak{I}}_t^{\alpha }f-{\mathfrak{I}}_t^{\alpha }{I}_{n}f\right|\le C_{\alpha ,k}{n}^{-k},\forall n\in {\mathbb{Z}}^{+}.$$

(32)

In particular, the GBFA approximation of ${\mathfrak{I}}_t^{\alpha }f$ converges to the exact RLFI at an algebraic rate $O\left(n^{-k}\right)$ as $n\to \infty$ for every fixed $t\in {\mathbf{\Omega }}_1$.

Proof. 

By assumption $f\in C^k\left({\mathbf{\Omega }}_1\right)$ with bounded kth derivative (30). For polynomial approximation on ${\mathbf{\Omega }}_1$ with respect to an orthogonal basis associated with a Jacobi-type weight (in particular, SG polynomials), classical Jackson-type estimates and interpolation results imply that there exists a constant $C_{\mathrm{app}}>0$, independent of n, such that the SG interpolation error at the SGG nodes satisfies the algebraic bound (31) for every fixed $\lambda >-1/2$ [14,15].

To transfer the interpolation error to the RLFI, we use the boundedness of the RLFI as a linear operator on $C\left({\mathbf{\Omega }}_1\right)$. Using (28) and (31), we obtain

$$\left|{\mathfrak{I}}_t^{\alpha }f-{\mathfrak{I}}_t^{\alpha }{I}_{n}f\right|\le {\displaystyle \frac{t^{\alpha }}{\Gamma (\alpha +1)}}{C}_{\mathrm{app}}{n}^{-k}.$$

Setting

$$C_{\alpha ,k}:={\displaystyle \frac{C_{\mathrm{app}}}{\Gamma (\alpha +1)}}\underset{t\in {\mathbf{\Omega }}_1}{sup}{t}^{\alpha },$$

which is finite since $t^{\alpha }\le 1$ on ${\mathbf{\Omega }}_1$, yields (32).    □

Remark 3 

(Convergence Hierarchy for Smooth Functions). The GBFA method exhibits a classical smoothness-dependent convergence hierarchy in terms of its ITE:

• Globally analytic f ($\mu =1$): Super-exponential convergence $\mathcal{O}\left(e^{-cnlnn}\right)$ as shown by Theorem 1.
• Gevrey-regular f with $1<\mu <2$: Log-stretched-exponential convergence $\mathcal{O}\left(e^{-c{n}^{1/\mu }lnn}\right)$ as shown by Theorem 2. This naturally interpolates between algebraic and super exponential rates.
• Gevrey-regular f with $\mu \ge 2$ (and general $C^{\infty }$ data): Qualitative convergence guaranteed but no explicit rate from current analysis (Theorems 3 and 4). The restriction $1<\mu <2$ in Theorem 2 is essential because for $\mu \ge 2$, the leading term $(\mu -2)nlnn$ in the error exponent becomes non-negative, requiring different analysis techniques.
• Finite smoothness $f\in C^k\left({\mathbf{\Omega }}_1\right)$: Algebraic convergence $\mathcal{O}\left(n^{-k}\right)$ as shown by Theorem 5.

This hierarchy demonstrates that the GBFA method preserves the classical spectral property: convergence accelerates with increased smoothness of f, from algebraic for finitely differentiable functions to super-exponential for analytic functions. The RLFI operator does not alter this fundamental behavior but merely introduces a bounded scaling factor dependent on α and t.

3.5. Quadrature Truncation Error for GBFA-Based RLFI

In the GBFA framework for $\alpha \in {\mathbf{\Omega }}_1^{\circ }$, the SG–transformed integrals (6) appearing in the RLFI approximation are evaluated numerically using the $(n_q,{\lambda }_q)$–GBFA quadrature; cf. the SGIRV construction in [16, Algorithms 6–7]. This numerical step introduces a quadrature truncation error that combines with the interpolation error to determine the total accuracy of the method. The quadrature truncation error for SG–based RLFI integrands is rigorously characterized in [4, Theorems 3–4]. These results provide a closed-form expression for the quadrature truncation error in terms of the $(n_q+1)$th SG derivative of the composed integrand, and sharp asymptotic bounds describing its decay with respect to $n_q$, including the regimes $j\sim n_q$ and $j\gg n_q$. In particular, [4, Theorem 3] gives the exact quadrature error representation for each SG mode j, while [4, Theorem 4] establishes its asymptotic decay rate and the influence of the quadrature index ${\lambda }_q$. Finally, the quadrature error produced by the $(n_q,{\lambda }_q)$–GBFA rule combines with the SG interpolation error analyzed in Theorems 1–5 to yield the total truncation error of the GBFA method. A complete asymptotic bound for this combined error is already established in [4, Theorem 5]. For analytic functions (Theorem 1), the regularity assumptions required by [4, Theorem 5] are fully satisfied, so the total GBFA error inherits the super-exponential decay in n coming from the ITE together with the fast decay in $n_q$ associated with the $(n_q,{\lambda }_q)$–GBFA quadrature. In particular, as shown in [4, Theorem 4], the quadrature error exhibits a two-regime structure: it vanishes identically whenever $j<n_q+1$, and for $j\ge n_q+1$ its dominant terms decay exponentially in $n_q$ (up to algebraic prefactors depending only on $\lambda$ and ${\lambda }_q$) under the smallness conditions on $\frac{2t{\eta }^{1/\alpha -1}}{\alpha }$ and $\frac{t{\eta }^{1/\alpha -1}}{\alpha }$ specified in [4, Eqs. (41) and (43)], where $\eta \in {\mathbf{\Omega }}_1^{\circ }$. Consequently, for analytic f, the total truncation error of the GBFA method is typically dominated by the interpolation component and therefore decays super-exponentially in n, while the quadrature component either vanishes or becomes negligible for sufficiently large $n_q$ in the admissible parameter regime. Thus, the overall accuracy of the GBFA method is governed jointly by the smoothness class of f and the chosen interpolation and quadrature parameters: the interpolation component decays at the rate dictated by the regularity of f, whereas the quadrature component is either zero or decays at least exponentially in $n_q$ when the transformed integrand is sufficiently smooth (e.g., when $1/\alpha \in {\mathbb{Z}}^{+}$); otherwise, it exhibits algebraic convergence due to limited endpoint regularity.

Remark 4 

(On the Smoothness Assumptions in the GBFA Error Analysis [4]). The truncation–error analysis in the GBFA work [4] (Theorems 1–5 therein) is carried out under the single assumption that $f\in C^{n+1}\left({\mathbf{\Omega }}_1\right)$, since the Lagrange interpolation remainder formula requires the existence and boundedness of $f^{(n+1)}$. No distinction is made between analytic, Gevrey-regular, $C^{\infty }$, or finitely smooth data; all bounds are expressed in terms of the quantity $\parallel f^{(n+1)}{\parallel }_{L^{\infty }\left({\mathbf{\Omega }}_1\right)}$ together with asymptotic properties of SG polynomials and Gamma function estimates.

It is important to note that the super-exponential decay in [4] is generally valid for $1/\alpha \in {\mathbb{Z}}^{+}$ when the derivatives of f satisfy an analytic-type growth condition of the form $\parallel f^{\left(k\right)}\parallel \le C{B}^k$, which ensures that $\parallel f^{(n+1)}{\parallel }_{L^{\infty }\left({\mathbf{\Omega }}_1\right)}$ does not dominate the polynomial–asymptotic factors. For general $C^{n+1}$ functions, however, the derivatives may grow at rates comparable to $k!$ (Gevrey), $k{!}^{\mu }$ (higher-order Gevrey), or even faster, in which case the factor $\parallel f^{(n+1)}{\parallel }_{L^{\infty }\left({\mathbf{\Omega }}_1\right)}$ may prevent super-exponential decay.

The present paper refines the convergence analysis by establishing a complete hierarchy of smoothness classes–globally analytic, Gevrey-regular ($1\le \mu <2$), general $C^{\infty }$, and finitely smooth $C^k$–and derives the corresponding decay rates for the ITE in each case. This clarifies that super-exponential convergence of the total GBFA error is achieved only when two conditions hold simultaneously: (i) f is globally analytic (ensuring super-exponential ITE decay per Theorem 1), and (ii) $1/\alpha \in {\mathbb{Z}}^{+}$, under which the classical transformation (2) renders the SG–transformed integrands (6) real-analytic (indeed, polynomial if f is polynomial). In this special regime, the $(n_q,{\lambda }_q)$–GBFA quadrature is exact for $n_q\ge n$, eliminating the quadrature truncation error entirely and yielding full super-convergence of the total error–even for functions whose derivatives grow rapidly–as confirmed in Theorem A10.

3.6. Endpoint Regularity of the Transformed Integrand for $\alpha >1$

The GBFA method relies on the transformation (2), which eliminates the algebraic kernel factor ${(t-\tau )}^{\alpha -1}$ in the RLFI. After applying (2), the RLFI integrand becomes

$$g\left(y\right):=f\left(t\left(1-y^{1/\alpha }\right)\right),y\in {\mathbf{\Omega }}_1,$$

(33)

which appears inside the SG–transformed integrals (6). Even when f is assumed smooth on ${\mathbf{\Omega }}_1$, the mapping $y\mapsto y^{1/\alpha }$ is not $C^1$ at $y=0$ when $\alpha >1$. This subsection clarifies the precise regularity of g near $y=0$ and its implications for the GBFA error structure.

Lemma 1 

(Regularity of $y^{1/\alpha }$ for $\alpha >1$). Let $\alpha >1$ and set $\beta :=1/\alpha \in {\mathbf{\Omega }}_1^{\circ }$. Then:

(i) 

$y\mapsto y^{\beta }$ is continuous on ${\mathbf{\Omega }}_1$ and analytic on ${\mathbf{\Omega }}_1^{\circ }$.

(ii) 

For every integer $m\ge 1$,

$${\displaystyle \frac{d^m}{d{y}^m}}{y}^{\beta }=c_{\beta ,m}{y}^{\beta -m},{\forall }_s{c}_{\beta ,m}\ne 0.$$

(iii) 

Since $\beta -m<0\forall m\ge 1$, all derivatives blow up at $y=0$:

$$\underset{y\to 0^{+}}{lim}\left|{\displaystyle \frac{d^m}{d{y}^m}}{y}^{\beta }\right|=+\infty .$$

In particular, since $|x^{1/\alpha }-y^{1/\alpha }{|\le |x-y|}^{1/\alpha }\forall x,y\in {\mathbf{\Omega }}_1$, the map $y\mapsto y^{\beta }$ is Hölder-continuous of order β on ${\mathbf{\Omega }}_1$, but its derivatives blow up at $y=0$; hence it belongs to $C^{0,\beta }\left({\mathbf{\Omega }}_1\right)$ but not to $C^1\left({\mathbf{\Omega }}_1\right)\forall \alpha >1$.

Theorem 6 

(Regularity of the Transformed Integrand g for $\alpha >1$). Let $\alpha >1$ and define g by (33). Assume that $f\in C^{\infty }\left({\mathbf{\Omega }}_1\right)$. Then $g\in C\left({\mathbf{\Omega }}_1\right)\setminus C^1\left({\mathbf{\Omega }}_1\right)$ and real-analytic on ${\mathbf{\Omega }}_1^{\circ }$.

Proof. 

Let $\alpha >1$ and set $\beta :=1/\alpha \in {\mathbf{\Omega }}_1^{\circ }$. Define the inner map

$$\varphi \left(y\right):=t\left(1-y^{\beta }\right),y\in {\mathbf{\Omega }}_1.$$

By Lemma 1, the function $y\mapsto y^{\beta }$ is continuous on ${\mathbf{\Omega }}_1$ and analytic on ${\mathbf{\Omega }}_1^{\circ }$, with all derivatives blowing up at $y=0$. Since $f\in C^{\infty }\left({\mathbf{\Omega }}_1\right)$, the composition $g\left(y\right)=f\left(\varphi \left(y\right)\right)$ is analytic on ${\mathbf{\Omega }}_1^{\circ }$ by the chain rule for analytic functions.

For the global regularity on ${\mathbf{\Omega }}_1$, we analyze the derivatives of g at $y=0$. By Faá di Bruno’s formula for the mth derivative of a composition $g=f\circ \varphi$ [17], we have

$$g^{\left(m\right)}\left(y\right)=\sum _{k=1}^m{f}^{\left(k\right)}\left(\varphi \left(y\right)\right)B_{m,k}\left({\varphi }^{\prime }\left(y\right),{\varphi }^{\prime \prime }\left(y\right),\dots ,{\varphi }^{(m-k+1)}\left(y\right)\right),$$

where $B_{m,k}$ are the (partial) Bell polynomials in the derivatives of $\varphi$. Each term in this sum is a finite linear combination of monomials of the form

$$f^{\left(k\right)}\left(\varphi \left(y\right)\right)\prod _{j=1}^m{\left({\varphi }^{\left(j\right)}\left(y\right)\right)}^{{\nu }_j},$$

where the nonnegative integers ${\nu }_1,\dots ,{\nu }_m$ satisfy

$$\sum _{j=1}^m{\nu }_j=k,\sum _{j=1}^{m}j{\nu }_j=m.$$

Using Lemma 1, each derivative ${\varphi }^{\left(j\right)}\left(y\right)$ behaves like

$${\varphi }^{\left(j\right)}\left(y\right)=-t{c}_{\beta ,j}{y}^{\beta -j},$$

so the product scales as

$$\prod _{j=1}^m{\left({\varphi }^{\left(j\right)}\left(y\right)\right)}^{{\nu }_j}=C_{\nu }{y}^{{\sum }_{j=1}^m(\beta -j){\nu }_j}=C_{\nu }{y}^{\beta k-m},$$

for some constant $C_{\nu }$ depending only on $t,\beta$ and the multi-index $\nu =({\nu }_1,\dots ,{\nu }_m)$. Thus, every Faà di Bruno term in $g^{\left(m\right)}\left(y\right)$ is asymptotically of the form

$$f^{\left(k\right)}\left(\varphi \left(y\right)\right)y^{\beta k-m}\mathrm{as}y\to 0^{+}.$$

To identify the most singular contribution, we note that for fixed m the exponent $\beta k-m$ is strictly increasing in k because $\beta >0$. Hence the smallest exponent (i.e., the most negative, and therefore the most singular) corresponds to the smallest possible k. Since $k\ge 1$, the dominant contribution arises from the terms with $k=1$. For such terms we have

$$\beta k-m=\beta -m,$$

and the corresponding part of $g^{\left(m\right)}\left(y\right)$ behaves like

$$g^{\left(m\right)}\left(y\right)\sim f^{\prime }\left(\varphi \left(y\right)\right)y^{\beta -m},y\to 0^{+}.$$

(34)

In particular, if $f^{\prime }\left(t\right)\ne 0$, then $f^{\prime }\left(\varphi \left(y\right)\right)\to f^{\prime }\left(t\right)\ne 0$ as $y\to 0^{+}$, and the factor $y^{\beta -m}$ controls the singular behavior of $g^{\left(m\right)}\left(y\right)$ near $y=0$. For $y^{\beta -m}$ to remain bounded as $y\to 0^{+}$, we require $\beta -m\ge 0$, i.e., $m\le \beta =1/\alpha$. Since $\beta <1$ for every $\alpha >1$, no integer $m\ge 1$ satisfies $m\le \beta$. Thus, unless $f^{\prime }\left(t\right)=0$, the dominant contribution (34) already becomes singular for $m=1$, and higher derivatives exhibit even stronger divergence. If $f^{\prime }\left(t\right)=0$, the leading singular term vanishes, but similar divergences reappear in higher-order contributions unless f possesses sufficiently high-order zeros at t. This establishes the proof.    □

Remark 5 

(Effect on GBFA Error: Quadrature vs. Interpolation). The finite endpoint regularity described in Theorem 6 affects only the quadrature error associated with approximating the SG–transformed integrals (6) using the $(n_q,{\lambda }_q)$–GBFA quadrature rule. For $\alpha >1$, the composed integrand g is merely Hölder-continuous of order $1/\alpha <1$ at $y=0$, which directly reduces the global regularity of g on ${\mathbf{\Omega }}_1$ and significantly slows the convergence of the quadrature truncation error component compared to the $0<\alpha <1$ case. This quadrature error dominates the total error for analytic f, limiting the overall convergence rate of the GBFA method and causing it to degrade as α increases (with larger α yielding slower quadrature convergence for fixed $n_q$). Importantly, the ITE in the original variable t remains unaffected and retains super-exponential decay in n. Thus, while the interpolation component preserves the full spectral rate dictated by the smoothness of f, the total GBFA accuracy for $\alpha >1$ is governed primarily by the typically algebraic quadrature decay, requiring substantially larger $n_q$ to maintain high accuracy as α grows, as confirmed numerically in Section 7.

3.7. Endpoint Regularity of the Transformed Integrand for $0<\alpha <1$, $1/\alpha \notin {\mathbb{Z}}^{+}$

We now analyze the SG–transformed integrand g defined in (33) for the regime $0<\alpha <1$. In contrast to the $\alpha >1$ case of Subsection 3.6, the inner map $\varphi \left(y\right)=t(1-y^{1/\alpha })$ exhibits higher endpoint smoothness when $1/\alpha \notin {\mathbb{Z}}^{+}$, and this improved regularity propagates to g.

Lemma 2 

(Regularity of $y^{1/\alpha }$ for $0<\alpha <1$, $1/\alpha \notin {\mathbb{Z}}^{+}$). Let $0<\alpha <1$ and $\beta :=1/\alpha >1$ with $\beta \notin {\mathbb{Z}}^{+}$. Then:

(i) 

$y\mapsto y^{\beta }$ is continuous on ${\mathbf{\Omega }}_1$ and analytic on ${\mathbf{\Omega }}_1^{\circ }$.

(ii) 

For every integer $m\ge 1$, the derivative formula from Lemma 1(ii) applies with $\beta >1$, and the exponent $\beta -m$ determines the endpoint behavior.

(iii) 

${\displaystyle \frac{d^m}{d{y}^m}}{y}^{\beta }$ remains finite at $y=0$ iff $m\le \lfloor \beta \rfloor$, and diverges for all $m\ge \lfloor \beta \rfloor +1$.

Hence $y^{\beta }\in C^{\lfloor \beta \rfloor ,\beta -\lfloor \beta \rfloor }\left({\mathbf{\Omega }}_1\right)$ and $y^{\beta }\notin C^{\lfloor \beta \rfloor +1}\left({\mathbf{\Omega }}_1\right)$.

Theorem 7 

(Regularity of the Transformed Integrand g for $0<\alpha <1$, $1/\alpha \notin {\mathbb{Z}}^{+}$). Let $0<\alpha <1$ and define g by (33). Assume $f\in C^{\infty }\left({\mathbf{\Omega }}_1\right)$. Then:

$$g\in C^{\lfloor 1/\alpha \rfloor ,1/\alpha -\lfloor 1/\alpha \rfloor }\left({\mathbf{\Omega }}_1\right)\setminus C^{\lfloor 1/\alpha \rfloor +1}\left({\mathbf{\Omega }}_1\right),$$

and g is analytic on ${\mathbf{\Omega }}_1^{\circ }$.

Proof. 

Let $\beta =1/\alpha >1$. The derivatives of $\varphi \left(y\right)=t(1-y^{\beta })$ are given by the same expressions used in the proof of Theorem 6, with the scaling ${\varphi }^{\left(j\right)}\left(y\right)\sim y^{\beta -j}$ inherited from Lemma 1. Applying Faà di Bruno’s formula [17] exactly as in Subsection 3.6 yields

$$g^{\left(m\right)}\left(y\right)\sim f^{\prime }\left(\varphi \left(y\right)\right)y^{\beta -m},y\to 0^{+},$$

since the dominant contribution again corresponds to $k=1$ in the Bell-polynomial expansion. Because $\beta >1$, the quantity $\beta -m$ is nonnegative precisely when $m\le \lfloor \beta \rfloor$, and negative for $m\ge \lfloor \beta \rfloor +1$. Thus $g^{\left(m\right)}\left(0\right)$ is finite for $m\le \lfloor \beta \rfloor$ and diverges otherwise, unless f has sufficiently high-order zeros at t. This proves the stated regularity.    □

Remark 6 

(Effect on GBFA Error for $0<\alpha <1$, $1/\alpha \notin {\mathbb{Z}}^{+}$). The endpoint smoothness in Theorem 7 improves the SG–quadrature behavior relative to the $\alpha >1$ case. Since g now possesses $C^{\lfloor 1/\alpha \rfloor }$ regularity at $y=0$, the $(n_q,{\lambda }_q)$–GBFA quadrature applied to (6) achieves a higher-order algebraic decay rate than that described in Remark 5. The ITE in t remains unchanged and retains super-exponential decay in n. Consequently, for $0<\alpha <1$, the total GBFA error is governed by a faster quadrature convergence that improves as α decreases.

4. Restoring Endpoint Regularity via a Generalized Transformation

For $\alpha >1$, the transformation (2) remains mathematically valid; however, as shown in Section 3.6, the mapping $y\mapsto y^{1/\alpha }$ is only Hölder-continuous of order $1/\alpha$ at $y=0$, and all derivatives blow up as $y\to 0^{+}$. This loss of endpoint regularity affects only the quadrature truncation error component of the GBFA method. For $0<\alpha <1$ with $1/\alpha \notin {\mathbb{Z}}^{+}$, the same mapping exhibits a complementary limitation: although $y^{1/\alpha }$ is differentiable at $y=0$, its higher-order derivatives fail to exist due to the noninteger exponent. Consequently, $\forall \alpha >0:1/\alpha \notin {\mathbb{Z}}^{+}$, the transformation (2) imposes endpoint smoothness restrictions that may reduce the global regularity of the transformed integrand. To improve the global smoothness of the transformed integrand for such values of $\alpha$, we introduce the generalized mapping

$$\tau =t\left(1-y^{\varphi \left(\alpha \right)}\right),y\in {\mathbf{\Omega }}_1,\varphi \left(\alpha \right)>0.$$

(35)

The mapping (35) is a $C^{\infty }$-diffeomorphism on ${\mathbf{\Omega }}_1^{\circ }$ for every $\alpha >0$ (see A.4), with

$${\displaystyle \frac{d\tau }{dy}}=-t\varphi \left(\alpha \right)y^{\varphi \left(\alpha \right)-1}.$$

Substituting (35) into the RLFI definition (1) yields the exact identity

(36)

Thus, for $\alpha >0$, the RLFI reduces to evaluating weighted SG-transformed integrals of the form

(37)

in place of the unweighted integrals (6). The GBFA approximation corresponding to (36) becomes

(38)

4.1. Weighted FSGIM Construction

Following the same steps used in [4, Section 3], we define

$${\mathbf{Y}}_{n_q}^{({\lambda }_q,\alpha )}:={\widehat{\mathbf{y}}}_{n_q}^{{{\lambda }_q}^{\circ p\left(\alpha \right)}}{\mathbf{1}}_{n+1}^{\top },{\tau }_{M,i,n_q}^{\left({\lambda }_q\right)}:=z_{M,i}\left({\mathbf{1}}_{n_q+1}-{\widehat{\mathbf{y}}}_{n_q}^{{{\lambda }_q}^{\circ \varphi \left(\alpha \right)}}\right),$$

and obtain the discrete GBFA approximation

(39)

for all $i\in {\mathbb{J}}_M$. This gives the weighted FSGIM defined row-wise by

(40)

which reduces to the original FSGIM in [4] when $\varphi \left(\alpha \right)=1/\alpha$.

5. Computational Complexity of the Weighted FSGIM

We analyze the computational complexity of the weighted FSGIM defined row-wise by (40). Let $n,n_q,M\in {\mathbb{Z}}^{+}$ denote the SG interpolation degree, the GBFA quadrature order, and the number of evaluation points, respectively. As in [4], the computation of the SG polynomials ${\widehat{G}}_{0:n}^{\left(\lambda \right)}$ using the three-term recurrence relation requires $\mathcal{O}\left(n\right)$ operations per degree. Since polynomial evaluations are required for all degrees up to n, this stage contributes $\mathcal{O}\left(n^2\right)$ operations per entry of ${\mathbf{z}}_M$.

For a fixed evaluation point $z_{M,i}$, the construction of requires evaluating the SG-based expressions at the $n_q+1$ transformed quadrature nodes, which incurs an additional $\mathcal{O}\left(n{n}_q\right)$ flops. The vector ${\widehat{\mathbf{y}}}_{n_q}^{{{\lambda }_q}^{\circ p\left(\alpha \right)}}$ is computed in $\mathcal{O}\left(n_q\right)$ flops, and forming the outer product ${\mathbf{Y}}_{n_q}^{({\lambda }_q,\alpha )}={\widehat{\mathbf{y}}}_{n_q}^{{{\lambda }_q}^{\circ p\left(\alpha \right)}}{\mathbf{1}}_{n+1}^{\top }$ requires $(n_q+1)(n+1)$ multiplications. The Hadamard product with ${{\mathcal{L}}_{0:n}^{\lambda }}^{\top }\left({\tau }_{M,i,n_q}^{\left({\lambda }_q\right)}\right)$ introduces another $(n_q+1)(n+1)$ multiplications. Applying $\widehat{\mathbf{P}}$ to collapse the $(n_q+1)$ rows into one row requires $(n+1)$ dot products of length $(n_q+1)$, i.e., $\mathcal{O}\left(n{n}_q\right)$ flops, and the scalar factor $\varphi \left(\alpha \right)z_{M,i}^{\alpha }/\Gamma \left(\alpha \right)$ adds $\mathcal{O}\left(n\right)$ flops.

Collecting the dominant terms, the overall per-row cost of constructing the weighted FSGIM is

$${\cos \mathrm{t}}_{\mathrm{row}}^{\mathrm{weighted}}=O\left(n^2+n{n}_q\right)=O\left(n(n+n_q)\right),$$

(41)

with constant factors at most a small multiple of the classical case. Therefore, constructing ${}^E\mathbf{Q}_n^{\alpha }\in {\mathbb{R}}^{(M+1)\times (n+1)}$ for all entries of ${\mathbf{z}}_M$ requires $O\left(Mn(n+n_q)\right)$ flops. Once ${}^E\mathbf{Q}_n^{\alpha }$ is available, applying the FSGIM to $f_{0:n}$ is a dense matrix–vector multiplication costing $O\left(Mn\right)$ flops in both the classical and weighted cases.

6. Endpoint Regularity and the Choice of $\varphi \left(\alpha \right)$

The transformed integrand in (37) is

$$g\left(y\right)=y^{p\left(\alpha \right)}f\left(t(1-y^{\varphi \left(\alpha \right)})\right),p\left(\alpha \right)=\alpha \varphi \left(\alpha \right)-1.$$

For analytic f, g is analytic on ${\mathbf{\Omega }}_1^{\circ }$, but its smoothness at $y=0$ depends on the competition between the vanishing weight $y^{p\left(\alpha \right)}$ and the singular derivatives of $y^{\varphi \left(\alpha \right)}$ as $y\to 0^{+}$. A Faá di Bruno–type exponent analysis shows that a sufficient condition for $g\in C^r\left({\mathbf{\Omega }}_1\right)$ is (A48); see Section A.5. This condition provides a practical design rule: choosing $\varphi \left(\alpha \right)$ according to (A48) guarantees that the transformed integrand is at least $C^r$ at $y=0$. Since Gauss–type quadrature for $C^r$ functions converges algebraically at a rate $\mathcal{O}\left(n_q^{-r}\right)$, the choice of r directly controls the quadrature truncation error.

Remark 7 

(Trade-offs in the Choice of $\varphi \left(\alpha \right)$). The mapping (35) enables a tunable balance between endpoint smoothness and weight concentration:

• Choosing $\varphi \left(\alpha \right)=1/\alpha$ yields perfect kernel cancellation and recovers the original GBFA identity (3), but the transformed integrand is only Hölder-continuous at $y=0$ for $\alpha >1$.
• Choosing $\varphi \left(\alpha \right)$ according to (A48) yields a $C^r$ integrand at $y=0$, improving the quadrature convergence rate to $\mathcal{O}\left(n_q^{-r}\right)$ by Theorem A12 and Corollary 1.
• Larger r increases smoothness but also increases the weight exponent $p\left(\alpha \right)=\alpha \varphi \left(\alpha \right)-1$, which may over-concentrate the integrand near $y=1$; cf. Figure 2. This concentration can degrade quadrature accuracy unless both ${\lambda }_q$ and $n_q$ are appropriately chosen. Specifically, ${\lambda }_q$ controls the node distribution, while $n_q$ must be sufficiently large to resolve the sharp peak near $y=1$${\forall }_{sl}r$.

Although the generalized $\varphi \left(\alpha \right)$-framework and its associated smoothness condition (A48) were developed primarily for $\alpha >1$, the underlying results extend naturally to all $\alpha >0$. In particular, for $0<\alpha <1$ with $1/\alpha \notin {\mathbb{Z}}^{+}$, one may enforce $C^r$ regularity via (A48), guaranteeing algebraic quadrature convergence $\mathcal{O}\left(n_q^{-r}\right)$. However, this approach is not recommended in practice. For $0<\alpha <1$, the classical transformation $\varphi \left(\alpha \right)=1/\alpha$ yields the exact, weight-free identity (3) and, for real-analytic f, induces a fractional power series expansion with exponentially decaying coefficients (cf. Remark A9). This structure leads to dramatically faster pre-asymptotic quadrature convergence–often appearing as high-order algebraic or near-exponential decay–without introducing artificial weighting or endpoint concentration near $y=1$. Consequently, the classical map is superior for practical $n_q$ (small and moderate) when $0<\alpha <1$, and the generalized $\varphi \left(\alpha \right)$-framework should be reserved ${\forall }_{sl}{n}_q$ or for $\alpha >1$.

Corollary 1 

(Effect of (A48) on the Total GBFA Error). Let f satisfy the hypotheses of Theorem 1. Let $\varphi \left(\alpha \right)$ satisfy the smoothness condition (A48) for some $r\ge 0$. Then the total GBFA error decomposes as

$$\left|{\mathfrak{I}}_t^{\alpha }f-{}^E\mathbf{Q}_n^{\alpha }{f}_{0:n}\right|\le \underset{\mathrm{ITE}}{\underbrace{\left|{\mathfrak{I}}_t^{\alpha }f-{\mathfrak{I}}_t^{\alpha }{I}_{n}f\right|}}+\underset{\mathrm{quadrature}\mathrm{error}}{\underbrace{C_{\alpha ,r}{n}_q^{-r}}},$$

where the ITE term decays super-exponentially in n by Theorem 1, and the quadrature error decays algebraically in $n_q$ with rate $\mathcal{O}\left(n_q^{-r}\right)$ due to the $C^r$ endpoint regularity enforced by (A48). In particular, the choice of $\varphi \left(\alpha \right)$ determines the dominant error component for fixed n.

Remark 8 

(Spectral Convergence Under Tunable Endpoint Regularity). Although the quadrature error in Corollary 1 decays only algebraically as $\mathcal{O}\left(n_q^{-r}\right)$, the ITE retains super-exponential decay $\mathcal{O}\left(e^{-cnlnn}\right)$ for analytic f by Theorem 1. Consequently, for any fixed $r\ge 0$, one may choose $n_q$ sufficiently large so that the quadrature error is negligible. Thereafter, $\exists N\in {\mathbb{Z}}^{+}$ such that $\forall n\ge N$, the total error is dominated by the ITE and thus decays faster than any algebraic rate $\mathcal{O}\left(n^{-k}\right)$, $\forall k\in {\mathbb{Z}}^{+}$. Hence, the method is spectrally convergent in n.

To complement the theoretical results, we now present numerical experiments illustrating the performance of the generalized GBFA method.

7. Numerical Experiments

We present comprehensive numerical experiments validating the GBFA method for arbitrary fractional orders $\alpha >0$, using both the classical transformation $\varphi \left(\alpha \right)=1/\alpha$ and the $\varphi \left(\alpha \right)$-generalized transformation developed in Section 4. All experiments use the same implementation framework of [4] with the $\varphi \left(\alpha \right)$ parameterization added for $\alpha >1$.

7.1. Experimental Setup

All simulations were performed in MATLAB R2023b on a 2.9 GHz AMD Ryzen 7 4800H processor with 16 GB RAM running Windows 11 (64-bit). Our primary test function is the cubic polynomial

$$f\left(t\right)=2t^3+8t,$$

(42)

for which the exact RLFI admits the closed form

$${\mathfrak{I}}_t^{\alpha }f={\displaystyle \frac{12t^{\alpha +3}}{\Gamma (\alpha +4)}}+{\displaystyle \frac{8t^{\alpha +1}}{\Gamma (\alpha +2)}}.$$

(43)

Unless otherwise stated, we use the following baseline parameters:

• SG interpolation degree $n=8$,
• SG parameter $\lambda =1$,
• Quadrature parameter ${\lambda }_q=0.5$,
• Evaluation point $t=0.5$,
• Quadrature order $n_q=16$,
• Fractional orders $\alpha \in \{1.5:9.5\}$.
• For the $\varphi \left(\alpha \right)$-generalized transformation, we take

  $$\varphi \left(\alpha \right)=\left\{\begin{array}{cc}{\displaystyle \frac{1}{\alpha }},\hfill & 0<\alpha <1,\hfill \\ {\displaystyle \frac{r+1}{\alpha }},\hfill & \alpha \ge 1,\hfill \end{array}\right.$$

  with smoothness parameter $r=12$, which guarantees $C^{12}\left({\mathbf{\Omega }}_1\right)$ regularity of the transformed integrand at $y=0$ for $f\in C^{12}\left({\mathbf{\Omega }}_1\right)$.

For each configuration, we compute the LARE defined as:

$$log\left(\mathrm{ARE}\right)=log\left({\displaystyle \frac{|{\mathfrak{I}}_t^{\alpha }f-{\mathfrak{I}}_t^{\alpha }{I}_{n}f|}{\left|{\mathfrak{I}}_t^{\alpha }f\right|}}\right).$$

7.2. Numerical Validation of the Classical Transformation $\tau =t(1-y^{1/\alpha })$

To support the theoretical developments established in the previous sections, we now present numerical experiments that assess the practical behavior of the classical transformation (2) for a wide range of fractional orders. These tests verify that the GBFA interpolation error in t retains its spectral character for smooth data, whereas the total approximation error for $\alpha >1$ is governed primarily by the quadrature component associated with the SG–transformed integrals (6). The results further illustrate the stability of the transformed integrand, the consistency of the computed FSGIMs, and the robustness of the method under variations in $\alpha$ and the SG parameters, while also reflecting the reduced endpoint regularity of $y\mapsto t(1-y^{1/\alpha })$ for $\alpha >1$.

Figure 3 presents a comprehensive six-panel analysis using the classical transformation $\varphi \left(\alpha \right)=1/\alpha$, systematically exploring the method’s behavior across parameter variations.

Figure 3(a) demonstrates good agreement between GBFA approximations and exact RLFI values across $\alpha \in {\mathbf{\Omega }}_{1.5,9.5}$. The GBFA values closely track the exact values over five orders of magnitude, from about $0.44$ at $\alpha =1.5$ to nearly $4.65\times {10}^{-10}$ at $\alpha =9.5$.

Figure 3(b) reveals a key phenomenon: the LARE increases monotonically from $-4.6961$ at $\alpha =1.5$ to $-2.9134$ at $\alpha =9.5$. This degradation in accuracy with increasing $\alpha$ is primarily due to the reduced endpoint regularity of the transformed integrand. As shown in Section 3.6, for $\alpha >1$, the mapping $y\mapsto y^{1/\alpha }$ is only Hölder-continuous of order $1/\alpha$ at $y=0$, with all derivatives blowing up as $y\to 0^{+}$. This endpoint singularity becomes more severe as $\alpha$ increases (since $1/\alpha$ decreases), slowing the convergence of the quadrature rule for fixed $n_q$.

Figure 3(c) investigates convergence with respect to the interpolation degree n for multiple $\alpha$ values. The data reveal a characteristic exactness plateau:

• For $n=2$: LARE ranges from approximately $-0.7$ ($\alpha =9.5$) to $-2.3$ ($\alpha =1.5$).
• For $n=4$: LARE improves dramatically to the range $[-2.9,-4.7]$.
• For $n\ge 4$: No further improvement occurs for any $\alpha$ value.

This plateau arises because f is a cubic polynomial, and SG interpolation at SGG nodes becomes exact for $n\ge 3$. The observed behavior perfectly aligns with Theorem 1: for polynomial data, the interpolation error vanishes once n exceeds the polynomial degree. The remaining error stems entirely from quadrature approximation of the SG–transformed integrals (6). The different colored curves in Figure 3(c) show that this plateau behavior is consistent across all tested $\alpha$ values.

Figure 3(d) reveals a striking result: for $\alpha >1$, the GBFA error is essentially independent of $\lambda$ across the tested range $\lambda \in {\mathbf{\Omega }}_{-0.4,2}$. For each fixed $\alpha$, all data points at different $\lambda$ values collapse onto horizontal lines, indicating constant LARE as $\lambda$ varies. This $\lambda$-insensitivity aligns with the theoretical analysis in Section 3.3. The bounding constant ${\vartheta }_{\alpha ,\lambda }$ in (13) contains the factor ${\alpha }^{-{\displaystyle \frac{1}{2}}-\alpha }$, which decays super-exponentially as $\alpha$ increases, effectively suppressing $\lambda$-dependence. This simplifies practical implementation: for $\alpha >1$, any $\lambda >-0.5$ yields comparable accuracy, though one should avoid values too close to $-0.5^{+}$ in practice since SG polynomials exhibit blow-up behavior as their indices approach this theoretical lower bound [4].

Figure 3(e) demonstrates significant sensitivity to ${\lambda }_q$, contrasting with the $\lambda$-insensitivity. The data reveal several important patterns:

• The optimal ${\lambda }_q$ is consistently ${\lambda }_q=0$ (Chebyshev quadrature) across all tested $\alpha$ values, yielding the smallest LARE. Thus, ${\lambda }_q=0$ provides robust performance with minimal sensitivity to variations in $\alpha$.
• Smaller ${\lambda }_q$ values (${\lambda }_q=-0.4,-0.2$) or larger ${\lambda }_q$ values (${\lambda }_q=1,2$) degrade accuracy substantially.

Figure 3(f) presents the most critical convergence study: varying $n_q$ while fixing $n=8$. The data reveal clear algebraic convergence patterns:

$$\begin{array}{cc}\hfill \alpha =1.5:& LARE\mathrm{improves}\mathrm{from}-2.8\mathrm{to}-5.7\mathrm{as}{n}_q\mathrm{increases}\hfill \\ \hfill & \mathrm{from}4\mathrm{to}32.\hfill \\ \hfill \alpha =4.5:& LARE\mathrm{improves}\mathrm{from}-1.9\mathrm{to}-3.9\mathrm{over}\mathrm{the}\mathrm{same}\hfill \\ \hfill & \mathrm{range}.\hfill \\ \hfill \alpha =9.5:& LARE\mathrm{improves}\mathrm{from}-1.7\mathrm{to}-3.6.\hfill \end{array}$$

The log-log plots in Figure 3(f) display generally near-linear trends, consistent with algebraic convergence $\mathcal{O}\left(n_q^{-p}\right)$. Average convergence exponents are $p\approx 3.2$ for $\alpha =1.5$ and $p\approx 2.1$ for $\alpha =9.5$, with convergence slowing as $\alpha$ increases. This algebraic convergence directly reflects the endpoint regularity limitations analyzed in Section 3.6. For $\alpha >1$, the transformed integrand g defined by (33) belongs to $C^{0,\beta }\left({\mathbf{\Omega }}_1\right)$ with $\beta =1/\alpha <1$, exhibiting only Hölder continuity at $y=0$. This reduced regularity fundamentally limits quadrature convergence to algebraic rates, consistent with Theorem 6 and Remark 5. The convergence becomes slower as $\alpha$ increases (smaller $\beta$), explaining why the $\alpha =9.5$ curve in Figure 3(f) has shallower slope than the $\alpha =1.5$ curve.

The numerical results in Figure 3 align precisely with the theoretical framework when considering the complete error composition for the classical transformation $\varphi \left(\alpha \right)=1/\alpha$: (i) Interpolation Component behaves as predicted by Theorem 1: for polynomial f, the interpolation error vanishes for $n\ge \mathrm{degree}\left(f\right)+1$, producing the observed exactness plateau in Figure 3(c); (ii) Quadrature Component dominates the total error for $\alpha >1$, exhibiting algebraic convergence $\mathcal{O}\left(n_q^{-p}\right)$ with $p>1$ (Figure 3(f)). The convergence exponent p decreases with increasing $\alpha$, from $p\approx 3.2$ at $\alpha =1.5$ to $p\approx 2.1$ at $\alpha =9.5$. This aligns perfectly with the endpoint regularity analysis in Section 3.6: the transformed integrand’s Hölder regularity $\beta =1/\alpha$ limits quadrature convergence rates.

The parameter sensitivity patterns match theoretical predictions:

• $\lambda$-insensitivity for $\alpha >1$ (Figure 3(d)) follows from ${\vartheta }_{\alpha ,\lambda }$ structure.
• The sensitivity to ${\lambda }_q$ (Figure 3(e)) indicates that quadrature accuracy can be tuned by proper parameter selection.
• $n_q$ scaling requirements increase with $\alpha$ due to worsening endpoint regularity (Figure 3(f)).

With baseline parameters ($n=8$, $n_q=16$), the method achieves $\mathrm{ARE}\sim {10}^{-5}$ to ${10}^{-3}$ across $\alpha \in {\mathbf{\Omega }}_{1.5,9.5}$ (Figure 3(b)). While not reaching machine precision, this represents useful accuracy for many applications, achievable with modest computational cost.

The high-resolution results in Figure 4 (using classical transformation $\varphi \left(\alpha \right)=1/\alpha$) provide deeper insight into the GBFA method’s behavior with extended quadrature resolution. Panel (a) demonstrates that with $n_q=1000$, the LARE values for $n\ge 4$ are substantially improved compared to those in Figure 3(c), reaching remarkably low values ranging from approximately $-6.8$ to $-10.6$ across different $\alpha$. This confirms that with sufficiently large $n_q$, the GBFA method achieves high accuracy for polynomial data, with the exactness plateau reaching LARE values below $-6.8$ across all tested $\alpha$ values.

Panel (b) reaffirms the $\lambda$-insensitivity observed in Figure 3(d), with all curves for each $\alpha$ collapsing to constant LARE values across the tested $\lambda$ range ${\mathbf{\Omega }}_{-0.4,2}$. This consistency across vastly different quadrature resolutions ($n_q=16$ versus $n_q=1000$) strengthens the conclusion that for $\alpha >1$, the SG interpolation parameter $\lambda$ has negligible effect on the total GBFA error, in agreement with the theoretical analysis in Section 3.3.

Panel (c) demonstrates that Chebyshev quadrature $({\lambda }_q=0)$ yields (near) optimal accuracy across all $\alpha$ values, achieving LARE values ranging from $-7.6$ (for $\alpha =9.5$) to $-10.8$ (for $\alpha =1.5$). Larger values $({\lambda }_q=1,2)$ progressively degrade the accuracy, with the LARE increasing more noticeably for larger $\alpha$. The sensitivity to ${\lambda }_q$ is more prominent in this high resolution setting ($n_q=1000$) than in Figure 3(e) because the quadrature error remains the dominant component of the total error, and ${\lambda }_q$ directly governs the convergence rate of this quadrature error for the transformed integrand with its specific endpoint regularity.

Panel (d) offers the most critical insight: with $n_q$ extended up to ${10}^4$, the GBFA method achieves LARE values as low as $-13.6$ (for $\alpha =1.5$) to $-9.0$ (for $\alpha =9.5$), corresponding to relative errors on the order of ${10}^{-14}$ to ${10}^{-9}$. The logarithmic x-axis clearly reveals that the convergence with respect to $n_q$ follows distinct patterns: rapid initial improvement for $n_q<100$, followed by progressively slower gains. For $\alpha =1.5$, LARE improves from $-5.7$ at $n_q=32$ to $-13.6$ at $n_q={10}^4$, demonstrating that machine-precision accuracy is attainable with sufficiently large $n_q$. However, the diminishing returns as $n_q$ increases highlight the algebraic nature of quadrature convergence for $\alpha >1$, consistent with the endpoint regularity constraints analyzed in Section 3.6.

7.3. Numerical Validation of the $\varphi \left(\alpha \right)$-Generalized Transformation

This section presents numerical experiments using the $\varphi \left(\alpha \right)$-generalized transformation (35) with smoothness parameter $r=12$. All tests use the polynomial test function $f\left(t\right)=2t^3+8t$ and its exact RLFI (43) at $t=0.5$.

Figure 5 presents a comprehensive six-panel analysis of the GBFA method with the $\varphi \left(\alpha \right)$-generalized transformation. Panel (a) demonstrates perfect agreement between GBFA approximations and exact RLFI values across $\alpha \in {\mathbf{\Omega }}_{1.5,9.5}$, with the data points being numerically identical (to near machine precision). The RLFI magnitude decreases over five orders of magnitude as $\alpha$ increases from 1.5 to 9.5.

Panel (b) reveals LARE values ranging from $-14.109$ at $\alpha =1.5$ to $-14.423$ at $\alpha =9.5$, with a minimum of $-15.104$ at $\alpha =2.5$, corresponding to relative errors of approximately ${10}^{-14.1},{10}^{-14.4}$, and ${10}^{-15.1}$, respectively. This demonstrates that with $r=12$, the $\varphi \left(\alpha \right)$-generalized transformation achieves consistent near machine-precision accuracy across all tested $\alpha$ values, with all LARE values below $-14.1$.

Panel (c) shows convergence with respect to the interpolation degree n. For $n=2$, LARE ranges from $-2.2935$ at $\alpha =1.5$ to $-0.69605$ at $\alpha =9.5$. At $n=4$, there is dramatic improvement to values ranging from $-14.102$ at $\alpha =1.5$ to $-15.954$ at $\alpha =9.5$, after which the error plateaus. This exactness plateau occurs because f is cubic, making SG interpolation exact for $n\ge 3$, leaving only quadrature error. With the $\varphi \left(\alpha \right)$-generalized transformation and $r=12$, this quadrature error is already at or near machine precision for the tested $n_q=16$, resulting in the observed LARE values near $-15$.

Panel (d) shows that under the $\varphi \left(\alpha \right)$-generalized transformation with $r=12$, the $\lambda$-sensitivity of the GBFA method exhibits a markedly different behavior from that observed under the classical mapping. By enforcing $C^{12}$ endpoint regularity, the quadrature error is reduced to near machine precision even with modest $n_q=16$, thereby exposing the ITE, which depends explicitly on $\lambda$ through ${\vartheta }_{\alpha ,\lambda }$ and the conditioning of the SG basis. This stands in contrast to the classical-transformation plots in Figure 3(d) and Figure 4(b), where the quadrature error dominates and masks the underlying $\lambda$-dependence.

In the ITE-dominated regime of the $\varphi \left(\alpha \right)$-generalized mapping, optimal $\lambda$ values (typically near $\lambda =0$ or $\lambda =0.5$) yield improving accuracy as $\alpha$ increases, consistent with the decay of the factor ${\alpha }^{-{\displaystyle \frac{1}{2}}-\alpha }$ in ${\vartheta }_{\alpha ,\lambda }$. Suboptimal choices, especially those approaching the theoretical lower bound $\lambda =-0.5^{+}$, suffer from ill-conditioning of the SG basis and produce significantly larger errors. Consequently, the observed $\lambda$-sensitivity band widens from roughly $0.11$ at $\alpha =1.5$ to nearly $1.9$ at $\alpha =8.5$. This widening is not a numerical artifact but a deterministic manifestation of the interplay between the decaying error constant and the $\lambda$-dependent stability of the spectral basis. Notably, this $\lambda$-sensitivity only emerges under the $\varphi \left(\alpha \right)$-generalized transformation; under the classical mapping $\varphi \left(\alpha \right)=1/\alpha$, the GBFA error remains practically $\lambda$-insensitive for $\alpha >1$ due to the dominance of quadrature error and the super-exponential decay of ${\vartheta }_{\alpha ,\lambda }$.

Panel (e) demonstrates that the quadrature parameter ${\lambda }_q$ plays a decisive role in the accuracy of the $\varphi \left(\alpha \right)$-generalized GBFA. The choice ${\lambda }_q=0.5$ is consistently optimal across all tested fractional orders, yielding LARE values of $-14.109$ at $\alpha =1.5$, $-15.104$ at $\alpha =2.5$, and $-14.423$ at $\alpha =9.5$. In contrast, suboptimal values such as ${\lambda }_q=-0.4$ or ${\lambda }_q=2$ degrade accuracy by 8–10 orders of magnitude. This striking disparity highlights the critical importance of selecting ${\lambda }_q=0.5$ when using the $\varphi \left(\alpha \right)$-generalized transformation, and we recommend this choice as a default setting for robust high-accuracy performance.

Panel (f) demonstrates rapid algebraic convergence with respect to quadrature points $n_q$. For $\alpha =1.5$, LARE improves from $-0.7406$ at $n_q=4$ to $-3.8205$ at $n_q=8$, $-7.5804$ at $n_q=12$, reaching near machine-precision level $-14.109$ at $n_q=16$, and stabilizing at $-14.941$ at $n_q=20$. The convergence accelerates for larger $\alpha$: for $\alpha =4.5$, LARE reaches the machine-precision at $-15.154$ by $n_q=12$, then slightly bounces back to $-14.911$ at $n_q=16$. For $\alpha =9.5$, near-optimal accuracy of $-14.463$ is achieved at $n_q=12$, stabilizing at $-14.423$ at $n_q=16$. This demonstrates that with $r=12$, modest $n_q$ values of 12–16 suffice to achieve near machine-precision accuracy, with some cases even exhibiting superconvergence (accuracy temporarily exceeding the final plateau).

Figure 6 presents a three-dimensional sensitivity analysis of the $\varphi \left(\alpha \right)$-generalized transformation with respect to the smoothness parameter r over the range $r\in \{0:2:12\}$.

Panel (a) shows LARE vs. r vs. ${\lambda }_q$ for fixed $\alpha =5.5$. The surface reveals a strong interaction: at optimal ${\lambda }_q=0.5$, LARE improves from $-3.12$ at $r=0$ to $-14.64$ at $r=12$, a gain of over $11.5$ orders of magnitude. Conversely, for suboptimal ${\lambda }_q=-0.4$ and ${\lambda }_q=2$, LARE degrades from $-2.77$ to $-9.58$ and from $-2.54$ to $-7.08$, respectively, as r increases from 0 to 12, demonstrating that inappropriate ${\lambda }_q$ combined with high r concentrates quadrature weight too aggressively near $y=1$, harming accuracy. These results indicate that the effectiveness of increasing r is highly sensitive to the quadrature parameter ${\lambda }_q$. For appropriately selected ${\lambda }_q$ (e.g., ${\lambda }_q=0.5$), larger r markedly accelerates convergence, delivering near machine-precision accuracy. However, for suboptimal ${\lambda }_q$, increasing r initially improves accuracy up to a critical value, beyond which further increases degrade performance—demonstrating that excessive smoothness without compatible quadrature design can be counterproductive.

Panel (b) shows LARE vs. r vs. $n_q$. The surface demonstrates accelerated convergence with higher r: while $r=0$ requires $n_q=32$ to reach LARE $\approx -14.6$, $r=12$ achieves LARE $<-15$ with merely $n_q=12$. This confirms Corollary 1, which predicts quadrature convergence $\mathcal{O}\left(n_q^{-r}\right)$ when the transformed integrand is $C^r$-smooth.

Panel (c) shows LARE vs. r vs. $\alpha$. The surface confirms that increasing r compensates for accuracy degradation at larger $\alpha$: for $r=0$, LARE deteriorates from $-4.70$ at $\alpha =1.5$ to $-2.91$ at $\alpha =9.5$, whereas for $r=12$, LARE remains consistently below $-14.1$ across all $\alpha \in {\mathbf{\Omega }}_{1.5,9.5}$, with the worst-case at $\alpha =8.5$ ($-14.37$) still representing near machine-precision accuracy. This demonstrates that the $\varphi \left(\alpha \right)$-generalized transformation effectively mitigates endpoint regularity limitations.

Collectively, these results demonstrate that while the GBFA method extends seamlessly to $\alpha >0$ using the classical transformation $\varphi \left(\alpha \right)=1/\alpha$, achieving high accuracy uniformly for all $\alpha >1$ with modest computational resources requires addressing the endpoint regularity limitation. The $\varphi \left(\alpha \right)$-generalized transformation introduced in Section 4 overcomes this by enforcing $C^r$ smoothness at $y=0$, achieving accelerated quadrature convergence $\mathcal{O}\left(n_q^{-r}\right)$ and enabling near machine-precision accuracy with practical $n_q$ values $\forall \alpha \in {\mathbf{\Omega }}_{1.5,9.5}$. The method exhibits strong sensitivity to the quadrature parameter ${\lambda }_q$: for $r=12$, the LARE at the optimal ${\lambda }_q=0.5$ is up to 10 orders of magnitude better than at highly suboptimal values. This dramatic difference underscores the critical importance of selecting ${\lambda }_q=0.5$—corresponding to Gauss–Legendre quadrature—when using the $\varphi \left(\alpha \right)$-generalized transformation. While $\lambda$-sensitivity exhibits a complex pattern—with optimal $\lambda$ values yielding excellent accuracy but performance degrading near $\lambda =-0.5^{+}$—the rapid $n_q$ convergence confirms the effectiveness of the $C^{12}$ regularization enforced by the $\varphi \left(\alpha \right)$ design condition (A48). The occasional superconvergence observed suggests interesting interactions between the quadrature rule and the regularized integrand structure.

7.3.1. Benchmarking GBFA Against Existing Solvers and Methods

Comparison with MATLAB’s Adaptive Quadrature

To benchmark the GBFA method with the $\varphi \left(\alpha \right)$-generalized transformation against a standard numerical approach, we compare its performance with MATLAB’s adaptive quadrature routine integral. The MATLAB function integral employs global adaptive quadrature with a relative tolerance of ${10}^{-14}$. For a fair comparison, we evaluate the RLFI ${\mathfrak{I}}_t^{\alpha }f$ for the cubic polynomial $f\left(t\right)=2t^3+8t$ at $t=0.5$ over the extended range $\alpha \in \{1.5:2:29.5\}$, using $r=14$ in the $\varphi \left(\alpha \right)$-generalized transformation to enforce $C^{14}\left({\mathbf{\Omega }}_1\right)$ regularity at $y=0$.

Figure 7 presents a comprehensive accuracy comparison between the GBFA method and MATLAB’s integral. Panel (a) shows that both methods produce RLFI values that are visually indistinguishable from the exact analytic solution across more than 41 orders of magnitude (from $\sim {10}^{-1}$ at $\alpha =1.5$ to $\sim {10}^{-42}$ at $\alpha =29.5$), confirming high accuracy for smooth data. Panel (b) displays the LARE for both methods: the GBFA method achieves LARE values in $[-15.72,-14.00]$, while MATLAB achieves LARE values in $[-15.95,-14.06]$, except at $\alpha =7.5$ where it attains machine-precision exactness (infinite LARE). Notably, both methods attain near machine-precision accuracy ($\mathrm{LARE}\le -14$) throughout the entire $\alpha$ range, with neither consistently dominating the other in accuracy.

The decisive advantage of the GBFA method lies in computational efficiency, as shown in Figure 8. Panel (a) reports computation times in microseconds ($\mu \mathrm{s}$) after FSGIM precomputation. The GBFA evaluation time ranges from $0.8$ to $104.7\mu \mathrm{s}$, with an average of $19.2\mu \mathrm{s}$. In contrast, MATLAB’s adaptive quadrature requires $192.9$–$2241.8\mu \mathrm{s}$, averaging $680.6\mu \mathrm{s}$. Panel (b) shows the speedup factor (MATLAB time / GBFA time), which ranges from $5.5\times$ at $\alpha =5.5$ to $1761.1\times$ at $\alpha =17.5$, with an average of $262.4\times$. The speedup remains substantial even for large $\alpha$, demonstrating the GBFA method’s scalability.

These results confirm that with $r=14$, the $\varphi \left(\alpha \right)$-generalized GBFA transformation achieves accuracy comparable to MATLAB’s high-tolerance adaptive quadrature while delivering orders-of-magnitude faster evaluation. This efficiency makes the method especially suitable for applications involving repeated RLFI evaluations or time-stepping schemes for fractional differential equations with arbitrary $\alpha >0$.

Comparison With SPH-Based Fractional Approximation

We compare the GBFA method with the 1D SPH-based fractional approximation of Ghorbani and Semperlotti [18]. Their scheme approximates left-handed constant-order RLFIs and derivatives via kernel-based particle quadrature on $[0,5]$, and is tested on

$$\begin{array}{cccc}\hfill f_1\left(x\right)& =sin\left(\pi x\right),\hfill & & f_2\left(x\right)=cos\left(\pi x\right),\hfill \\ \hfill f_3\left(x\right)& ={\mathrm{e}}^x,\hfill & & f_4\left(x\right)={(x-1)}^3,\hfill \end{array}$$

for the left-handed RLFI of order $\alpha =0.75$. Numerical results for ${\mathfrak{I}}_x^{0.75}{f}_j$ on ${\mathbf{\Omega }}_{0,5}$ are reported in [18, Fig. 5] for all $j\in {\mathbb{N}}_4$.

Restricting attention to $x\in {\mathbf{\Omega }}_1$ to match the present GBFA configuration, visual inspection of [18, Fig. 5] shows that the SPH approximation exhibits relatively large absolute errors:

$$\begin{array}{cc}\hfill \underset{x\in {\mathbf{\Omega }}_1}{max}|\mathrm{SPH}\left(x\right)-{\mathfrak{I}}_x^{0.75}{f}_1|& \gtrsim 2\times {10}^{-2},\hfill \\ \hfill \underset{x\in {\mathbf{\Omega }}_1}{max}|\mathrm{SPH}\left(x\right)-{\mathfrak{I}}_x^{0.75}{f}_2|& \gtrsim 5\times {10}^{-4},\hfill \\ \hfill \underset{x\in {\mathbf{\Omega }}_1}{max}|\mathrm{SPH}\left(x\right)-{\mathfrak{I}}_x^{0.75}{f}_3|& \gtrsim 5\times {10}^{-3},\hfill \\ \hfill \underset{x\in {\mathbf{\Omega }}_1}{max}|\mathrm{SPH}\left(x\right)-{\mathfrak{I}}_x^{0.75}{f}_4|& \gtrsim 2\times {10}^{-2}.\hfill \end{array}$$

Even after auxiliary-particle refinement, [18, Fig. 6] indicates

$$\underset{x\in {\mathbf{\Omega }}_1}{max}|{\mathrm{SPH}}_{\mathrm{ref}}\left(x\right)-{\mathfrak{I}}_x^{0.75}{f}_1|\gtrsim 3\times {10}^{-4}.$$

By contrast, the present GBFA method with $n=n_q=20$ and $\lambda ={\lambda }_q=0.5$, using the exact transformation (2), yields for the same four test functions and $\alpha =0.75$,

$$\underset{x\in {\mathbf{\Omega }}_1}{max}|\mathrm{GBFA}\left(x\right)-{\mathfrak{I}}_x^{0.75}{f}_j|<2\times {10}^{-7}\forall j\in {\mathbb{N}}_4,$$

with $L_2$-errors of the same order. Table 1 summarizes the MAE on ${\mathbf{\Omega }}_1$. The SPH entries are conservative visual estimates from [18, Fig. 5], while the GBFA entries follow directly from (3) and (7).

Beyond the MAE, the GBFA method achieves near-perfect goodness-of-fit. The computed $R^2$ scores exceed $0.9999999999999$ for all four functions, with specific values

$$\begin{array}{cc}\hfill R^2\left(f_1\right)& =0.999999999999879,\hfill \\ \hfill R^2\left(f_2\right)& =0.999999999999932,\hfill \\ \hfill R^2\left(f_3\right)& =0.999999999999984,\hfill \\ \hfill R^2\left(f_4\right)& =0.999999999999926.\hfill \end{array}$$

These values confirm the visual indistinguishability between the GBFA and exact RLFI curves in Figure 9.

In contrast, the SPH method reports significantly lower $R^2$ scores: $0.977319$ for $f_1$, $0.999987$ for $f_2$, $1.000000$ (rounded) for $f_3$, and $0.999994$ for $f_4$[18]. The modest $R^2$ for $f_1$ aligns with its large absolute error ($\gtrsim 2\times {10}^{-2}$), indicating a systematic deviation. Even for $f_3$, the apparently perfect SPH $R^2$ masks an absolute error above $5\times {10}^{-3}$, illustrating that $R^2$ alone may obscure significant discrepancies when the reference signal has low variance. GBFA, by contrast, achieves both near-unity $R^2$ and absolute errors below $2\times {10}^{-7}$.

Figure 9 and Figure 10 illustrate the GBFA performance. Figure 9 shows that the GBFA approximation of ${\mathfrak{I}}_x^{0.75}{f}_j$ is visually indistinguishable from the exact RLFI on ${\mathbf{\Omega }}_1$. Figure 10 displays the corresponding error distributions and statistical indicators: MAE, mean error, median error, $L_2$ accuracy, and $R^2$ score. All error curves remain below $2\times {10}^{-7}$, with $L_2$ errors $\le 1.8\times {10}^{-7}$ and $R^2>0.9999999999998$ for every test case.

The several-orders-of-magnitude accuracy gap between SPH and GBFA follows directly from their underlying approximation mechanisms. SPH employs kernel-based particle quadrature and exhibits, in practice, at most algebraic convergence with respect to particle spacing. GBFA, by contrast, exploits SG interpolation and FSGIMs within the exact transformed identity (3). Theorem 1 guarantees super-exponential decay of the ITE for analytic data, while the quadrature truncation error behaves as described in Section 3.5: its asymptotic rate is algebraic $\mathcal{O}\left(n_q^{-K-\theta }\right)$ for $0<\alpha <1$ with $1/\alpha \notin {\mathbb{Z}}^{+}$, where $K=\lfloor 1/\alpha \rfloor$ and $\theta =1/\alpha -K$. However, when f is real-analytic–as in the test functions $f_j$–the transformed integrand $g\left(y\right)=f\left(t(1-y^{1/\alpha })\right)$ admits a fractional power series with exponentially decaying coefficients, yielding dramatically faster pre-asymptotic quadrature convergence than the worst-case algebraic bound. For $\alpha =0.75$ ($1/\alpha =4/3\notin {\mathbb{Z}}^{+}$), this results in observed errors below $2\times {10}^{-7}$ with modest $n_q=20$, far surpassing SPH’s performance. Consequently, the combined GBFA error decays markedly faster in $(n,n_q)$ than the SPH error does in particle resolution, explaining the performance advantage observed in Table 1 and Figure 9 and Figure 10.

Figure 11 displays the convergence behavior of the GBFA method for the RLFI ${\mathfrak{I}}_x^{0.75}{f}_1$ at $x=0.5$, using the classical transformation $\varphi \left(\alpha \right)=1/\alpha$. The observed absolute errors decay rapidly from $2.912\times {10}^{-6}$ at $n_q=4$ to $2.166\times {10}^{-13}$ at $n_q=40$, yielding a fitted algebraic rate of $p=7.1252$ via least-squares regression in log–log space. This empirical rate significantly exceeds the theoretical asymptotic prediction $\mathcal{O}\left(n_q^{-K-\theta }\right)=\mathcal{O}\left(n_q^{-4/3}\right)$ given by Corollary A3(ii), where $K=\lfloor 1/\alpha \rfloor =1$ and $\theta =1/\alpha -K=1/3$. As explained in Section 3.5 and Remark A9, this discrepancy arises because $f_1$ is real-analytic on a neighborhood of ${\mathbf{\Omega }}_x$, inducing a fractional power series expansion

$$g\left(y\right)=f_1\left(x(1-y^{1/\alpha })\right)=\sum _{m=0}^{\infty }{a}_m{y}^{m/\alpha },\left|a_m\right|\le \tilde{C}{\tilde{B}}^m,$$

with exponentially decaying coefficients. Consequently, the tail of the series beyond mode $m\gtrsim n_q$ is extremely small, causing Gauss-type quadrature to exhibit dramatically faster pre-asymptotic convergence than the worst-case $\mathcal{O}\left(n_q^{-4/3}\right)$ rate–which governs only the limit $n_q\to \infty$ under minimal $C^{1,\theta }$ regularity. Thus, the near-exponential decay observed here is fully consistent with the refined smoothness analysis in Section A.3.

Generalized Transformation for $0<\alpha <1$ with $1/\alpha \notin {\mathbb{Z}}^{+}$

Although the generalized $\varphi \left(\alpha \right)$-transformation was primarily designed for $\alpha >1$, it can also be applied to $0<\alpha <1$ when $1/\alpha \notin {\mathbb{Z}}^{+}$ to enforce algebraic convergence $\mathcal{O}\left(n_q^{-r}\right)$ for any prescribed $r>0$. This capability follows directly from Theorem A12 and Corollary 1, which guarantee that choosing $\varphi \left(\alpha \right)=(r+1)/\alpha$ yields a transformed integrand in $C^r\left({\mathbf{\Omega }}_1\right)$, thereby ensuring quadrature error decay $\mathcal{O}\left(n_q^{-r}\right)$. To illustrate, we revisit the RLFI ${\mathfrak{I}}_x^{0.75}sin\left(\pi x\right)$ at $x=0.5$ with $r=15$. As shown in Figure 12, the absolute error decays slowly for small $n_q$ due to endpoint weight concentration near $y=1$, but transitions into the asymptotic regime $\mathcal{O}\left(n_q^{-r}\right)$ for $n_q\gtrsim 30$, achieving near machine precision by $n_q=40$. The fitted convergence rate $p\approx 14.33$ closely matches the theoretical expectation $p=r=15$. However, this accelerated convergence emerges only beyond a threshold quadrature order ($n_q\approx 30$); below this, the classical map $\varphi \left(\alpha \right)=1/\alpha$ yields superior accuracy due to its exact kernel cancellation and absence of artificial weighting. Moreover, for real-analytic f, the classical transformation induces a fractional power series with exponentially decaying coefficients, as we demonstrated earlier, leading to dramatically faster pre-asymptotic convergence than the worst-case algebraic bound. Consequently, the generalized transformation is not recommended for $0<\alpha <1$ with $1/\alpha \notin {\mathbb{Z}}^{+}$ unless (near) machine precision is required; in such cases, it may be employed to enforce a prescribed algebraic rate $\mathcal{O}\left(n_q^{-r}\right)$${\forall }_{sl}r,n_q$. Otherwise, the classical GBFA framework remains optimal both theoretically and computationally.

Comparison With Hybrid Function (HF) Operational Matrices

We now compare the GBFA method with the orthogonal HF approach of Damarla and Kundu [19], which constructs generalized one-shot operational matrices for the RLFI. Their method approximates the RLFI of $f\left(t\right)=t$ on $t\in {\mathbf{\Omega }}_1$ using a step size $h=0.125$ (i.e., $m=8$ subintervals) and reports the ∞-norm error against the exact RLFI

$${\mathfrak{I}}_t^{\alpha }f={\displaystyle \frac{t^{\alpha +1}}{\Gamma (\alpha +2)}}.$$

Following the experimental protocol of [19, Table 1], we evaluate the GBFA of the same RLFI at the HF grid points $t_{\mathrm{HF}}=\{0,h,2h,\dots ,1\}$ and compute the ∞-norm error over this discrete set. The GBFA parameters are fixed at $n=2$, $n_q=12$, $\lambda =1$, ${\lambda }_q=0.5$, and $r=15$.

Table 2 summarizes the comparison between both methods. The results demonstrate that both methods achieve machine-precision accuracy $\forall \alpha \in {\mathbf{\Omega }}_{0.5,5}$. Crucially, the computational cost of the GBFA evaluation is up to nearly three orders of magnitude lower. As shown in Table 2, the GBFA CPU time is consistently under 1 millisecond (ms), with a maximum of $0.76$ ms at $\alpha =0.5$. In contrast, the HF method requires $\sim 21.7$ ms on average for the same task, as reported in [19, Table 1]. This represents a speedup factor of $\sim 28\times$ to $\sim 2892\times$, computed as the ratio $\mathrm{HF}\mathrm{CPU}/\mathrm{GBFA}\mathrm{CPU}$.

This superior efficiency stems from the spectral nature of GBFA: it requires only $n+1=3$ function evaluations to achieve machine precision, whereas the HF method relies on a piecewise representation over $m=8$ intervals (16 coefficients). Thus, for repeated evaluations of the RLFI with a fixed $\alpha$, the GBFA method offers a dramatically faster online phase. This advantage is critical for applications involving fractional differential equations, where the RLFI operator is evaluated repeatedly in a time-stepping loop.

Comparison With Bernstein Approximation Method

We compare the GBFA method with the Bernstein approximation method of Usta [20] applied to solve the following second-kind FVIE:

$$u\left(x\right)=\sqrt{\pi }{(1+x)}^{-3/2}-0.02{\displaystyle \frac{x^3}{1+x}}+0.01x^{5/2}{\mathfrak{I}}_x^{1/2}u,x\in {\mathbf{\Omega }}_1,$$

(44)

with exact solution

$$u\left(x\right)=\sqrt{\pi }{(1+x)}^{-3/2},$$

(45)

which is real-analytic on ${\mathbf{\Omega }}_1$, with its nearest singularity located at $x=-1$. Although analytic, its derivatives grow factorially, placing it in the Gevrey-1 class (see A.6). The GBFA solver uses SGG nodes with $n=12$, $\lambda =0$, $n_q=12$, and ${\lambda }_q=0.5$, yielding a MAE of $6.080\times {10}^{-13}$, a CPU time of $3\times {10}^{-3}\mathrm{s}$, and a RmsE of $3.601\times {10}^{-13}$, which is defined as:

$$\mathrm{RmsE}=\sqrt{{\displaystyle \frac{1}{N}}\sum _{i=1}^N{\left|u\left(x_i\right)-u_{\mathrm{exact}}\left(x_i\right)\right|}^2},$$

with $x_{1:N}$ being the evaluation nodes and N is their total number.

As shown in Table 3, the GBFA method achieves a MAE of $6.08\times {10}^{-13}$ at $n=12$, which is two orders of magnitude smaller than the Bernstein result ($1.05\times {10}^{-11}$) at the same resolution, while requiring over 60 times less CPU time (3 ms vs. 186 ms). Although the Bernstein method is analyzed using a theoretical convergence bound of $\mathcal{O}\left(n^{-1}\right)$ [20, Eq. (1.7)], its empirical performance for smooth solutions is significantly better–reaching near machine precision at modest n. Nevertheless, the GBFA method maintains a clear advantage in both accuracy per degree of freedom and computational efficiency, underpinned by its provable super-exponential convergence (Theorem 1).

Although the exact solution u is analytic on ${\mathbf{\Omega }}_1$, its derivatives satisfy (A49), placing it in the Gevrey-1 class rather than the stricter analytic class assumed in Theorem 1; cf. Section A.6. This distinction does not affect the observed convergence because the radius of analyticity is large (distance 1 to the singularity at $x=-1$). Moreover, the Gevrey-$\mu$ convergence rate in Theorem 2 continuously recovers the analytic rate of Theorem 1 as $\mu \to 1^{+}$: the log-stretched-exponential decay $e^{-c{n}^{1/\mu }lnn}$ tends to the super-exponential form $e^{-cnlnn}$. Hence, Gevrey-1 functions exhibit asymptotically the same GBFA convergence as analytic ones. Crucially, the fractional order satisfies $1/\alpha =2\in {\mathbb{Z}}^{+}$, so by Theorem A10, the SG–transformed integrands are real-analytic (polynomial if f is polynomial), and the $(n_q,{\lambda }_q)$–GBFA quadrature is exact for $n_q\ge n$, eliminating quadrature truncation error. Consequently, the total error inherits the super-exponential ITE decay, yielding near machine-precision accuracy at modest resolutions–e.g., MAE $\approx 6.08\times {10}^{-13}$ at $n=12$.

Figure 13 visually confirms the near-indistinguishability of the GBFA solution from the exact one on the SGG grid.

Figure 14 displays the convergence rate of the GBFA method, where the MAE is plotted against the polynomial degree n on a semi-log scale. The data in Table 4 show a rapid decay of the MAE from $\approx 2.8\times {10}^{-8}$ at $n=6$ to $\approx 7.8\times {10}^{-16}$ at $n=16$, consistent with near machine-precision accuracy. A linear regression of $ln\left(\mathrm{MAE}\right)$ against $nlnn$ yields a slope of approximately $-0.54$, confirming the super-exponential convergence rate $\mathcal{O}\left(e^{-cnlnn}\right)$ predicted by Theorems 1 and 2 in the limit $\mu \to 1^{+}$. Although the exact solution u belongs to the Gevrey-1 class, its analyticity on ${\mathbf{\Omega }}_1$ and large radius of convergence ensure that its numerical behavior is indistinguishable from the analytic case, producing the observed super-exponential decay. This is fully consistent with the error hierarchy in Remark 3 and demonstrates the GBFA method’s exceptional efficiency for smooth problems.

Comparison With Quadratic and Cubic Spline-Based Schemes

We now compare the GBFA method with the Quadratic (S1) approximation scheme of Kumar et al. [23] for the RLFI. Their approach constructs piecewise polynomial interpolants over uniform grids with step size h, achieving theoretical convergence order $\mathcal{O}\left(h^3\right)$. The test function $f\left(t\right)=sin\left(t\right)$ at $t=1$ with $\alpha =1.5$ provides a direct benchmark against their reported results in [23, Table 3].

Using the $\varphi \left(\alpha \right)$-generalized GBFA transformation with $r=10$, $n=10$, $n_q=14$, $\lambda =0.1$, and ${\lambda }_q=0.5$, we compute ${\mathfrak{I}}_1^{1.5}sin$ with absolute error $3.045\times {10}^{-11}$ in $2.264$ ms. In contrast, the S1 scheme requires $h=6.25\times {10}^{-3}$ ($n=80$ i.e. 160 subintervals) to reach error $3.05977\times {10}^{-11}$. Thus, the GBFA method attains slightly higher accuracy than the S1 scheme, despite using far fewer degrees of freedom ($n+1=11$ vs. 161 function evaluations) and significantly lower computational effort.

This dramatic advantage stems from the spectral convergence of GBFA, as established in Theorem 1 for analytic f, versus the algebraic convergence inherent to spline-based methods. The classical GBFA transformation already ensures exact kernel cancellation via (3), while the generalized transformation (35) enforces $C^{10}$ endpoint regularity through (A48), thereby accelerating quadrature convergence per Corollary 1. Consequently, GBFA achieves near machine-precision accuracy where piecewise polynomial methods remain limited by their fixed convergence rates, even at fine mesh resolutions.

Having established both theoretical and numerical evidence, we conclude by summarizing the broader implications of the generalized GBFA framework.

8. Conclusions and Discussion

This work establishes that the GBFA method, originally formulated for $\alpha \in {\mathbf{\Omega }}_1^{\circ }$, extends naturally to all $\alpha \in {\mathbb{R}}^{+}$ through two complementary approaches: (i) the classical transformation (2) with $\varphi \left(\alpha \right)=1/\alpha$, yielding the exact identity (3) but suffering from reduced endpoint regularity for $\alpha >1$; and (ii) the $\varphi \left(\alpha \right)$-generalized transformation (35), producing the weighted exact form (36) with a tunable $\varphi \left(\alpha \right)$ that restores smoothness at $y=0$. The total GBFA error decomposes per Corollary 1 into an ITE governed by Theorems 1–5 and a quadrature component scaling as $\mathcal{O}\left(n_q^{-r}\right)$ when $\varphi \left(\alpha \right)$ satisfies (A48). Numerical experiments confirm that for $r=12$, the generalized transformation achieves near machine-precision accuracy (LARE $<-14$) across $\alpha \in {\mathbf{\Omega }}_{1.5,9.5}$ using only $n_q=16$ (Figure 5f), whereas the classical transformation requires $n_q\gtrsim {10}^4$ to attain comparable accuracy (Figure 4d).

A refined smoothness analysis in Appendix A.3 reveals that full super-exponential convergence–including exact quadrature for polynomial data–is achieved by the classical GBFA transformation iff $1/\alpha \in {\mathbb{Z}}^{+}$. In this case, the transformed integrand is real-analytic on ${\mathbf{\Omega }}_1$, and Gauss-type quadrature becomes exact for sufficiently large $n_q$. Otherwise, the integrand belongs only to the Hölder space $C^{K,\theta }\left({\mathbf{\Omega }}_1\right)$ with $K=\lfloor 1/\alpha \rfloor$ and $\theta =1/\alpha -K$, limiting the asymptotic quadrature rate to $\mathcal{O}\left(n_q^{-K-\theta }\right)$. Crucially, as clarified in Remark A9, when f is real-analytic, the transformed integrand admits a fractional power series with exponentially decaying coefficients, causing the pre-asymptotic quadrature error to decay dramatically faster than the worst-case algebraic bound, thereby reconciling the near-exponential convergence often observed numerically even when $1/\alpha \notin {\mathbb{Z}}^{+}$.

Importantly, while the generalized $\varphi \left(\alpha \right)$-framework was developed primarily for $\alpha >1$, our numerical study in Section 7 confirms its applicability to the regime $0<\alpha <1$ with $1/\alpha \notin {\mathbb{Z}}^{+}$. In such cases, the classical map already yields rapidly decaying pre-asymptotic errors due to exponential coefficient decay in the fractional power series; consequently, the generalized transformation is not recommended unless (near) machine precision is explicitly required, as it incurs degraded accuracy for moderate $n_q$ and introduces sensitivity to ${\lambda }_q$. When such extreme accuracy is needed, however, the generalized approach can enforce a prescribed algebraic rate $\mathcal{O}\left(n_q^{-r}\right)$ for all $r>0$, achieving the target precision beyond a threshold $n_q\approx 30$ (see Figure 12).

The extension to arbitrary $\alpha >0$ enables direct application of GBFA to a broad class of fractional models, grounded in the exact identities (1), (3), (36) and the error hierarchy of Theorems 1–5. In particular, the GBFA can be systematically applicable to:

1.

High-order fractional viscoelastic models ($\alpha >1$) [2], where the $\varphi \left(\alpha \right)$-transformation maintains accuracy with modest $n_q$.

2.

Super-diffusive transport ($\alpha >1$) [3], efficiently evaluated via FSGIMs (40).

3.

Integer-order repeated integrals ($\alpha =m\in {\mathbb{Z}}^{+}$), for which (36) provides a spectrally accurate alternative to classical quadrature.

4.

Variable-order operators $\alpha =\alpha \left(t\right)>0$, where precomputed FSGIMs ${}^E\mathbf{Q}_n^{\alpha }$ (classical) or (40) (generalized) can be tabulated and interpolated over $\alpha$-nodes.

Key limitations and directions for further work include: (i) For $f\in C^k\left({\mathbf{\Omega }}_1\right)$ (Theorem 5), a rigorous analysis of the interplay among $\varphi \left(\alpha \right)$, ITE, and quadrature error–particularly in the presence of endpoint singularities–is still lacking; (ii) For variable-order $\alpha \left(t\right)$, pointwise-adaptive $\varphi \left(\alpha \right(t\left)\right)$ satisfying (A48) ${\forall }_{a}t\in {\mathbf{\Omega }}_1$ may significantly improve time-integration efficiency.

The $\varphi \left(\alpha \right)$-generalized transformation (35) thus elevates GBFA from a method restricted to $0<\alpha <1$ to a unified, high-accuracy framework for arbitrary $\alpha \in {\mathbb{R}}^{+}$, with tunable smoothness and efficiency controlled by the phi-design strategy via (A48).

8.1. Practical Implementation Guidelines

While the analytical structure of GBFA is unchanged for $\alpha >0$, numerical behavior hinges critically on the choice between transformations. The following guidelines emerge from [4], Theorems 1–5, and the numerical experiments in Section 7.2 and Section 7.3.

• For $0<\alpha <1$, the classical transformation (2) with $\varphi \left(\alpha \right)=1/\alpha$ yields the exact identity (3). The smoothness of the transformed integrand depends critically on whether $1/\alpha \in {\mathbb{Z}}^{+}$:

  -

  If $1/\alpha \in {\mathbb{Z}}^{+}$, the integrand is real-analytic on ${\mathbf{\Omega }}_1$, and GBFA achieves full super-exponential convergence (including exact quadrature for polynomial data). Parameter recommendations from [4] then apply directly.

  -

  If $1/\alpha \notin {\mathbb{Z}}^{+}$, the integrand lies in $C^{K,\theta }\left({\mathbf{\Omega }}_1\right)$ with $K=\lfloor 1/\alpha \rfloor$, $\theta =1/\alpha -K$, limiting the asymptotic quadrature rate to $\mathcal{O}\left(n_q^{-K-\theta }\right)$. However, for real-analytic f, the fractional power series of the integrand has exponentially decaying coefficients (Remark A9), yielding dramatically faster pre-asymptotic convergence—often appearing as high-order algebraic or near-exponential decay for practical $n_q$. In this regime, the generalized $\varphi \left(\alpha \right)$-transformation is not recommended unless (near) machine precision is required (see Section 7).

  When $1/\alpha \in {\mathbb{Z}}^{+}$, parameter recommendations from [4] depend on the scale of n and $n_q$:

  -

  For relatively small n and $n_q$: any $\lambda \in {\mathcal{T}}_{c,r}$ is feasible:

  $${\mathcal{T}}_{c,r}=\left\{\gamma \mid -{\displaystyle \frac{1}{2}}+\epsilon \le \gamma \le r,0<\epsilon \ll 1,r\in {\mathbf{\Omega }}_{1,2}\right\}.$$

  The set ${\mathcal{T}}_{c,r}$ is often the effective operational range for the SG parameters used in SG-based collocation regimes; the subscript `c’ in ${\mathcal{T}}_{c,r}$ denotes “collocation.” Choosing smaller ${\lambda }_q$ from this range generally improves quadrature accuracy, with the notable exception that ${\lambda }_q=0.5$ often minimizes the quadrature error.

  -

  For relatively large n and $n_q$:

  *

  Precision computations: select $\lambda \in {\mathbf{\Omega }}_{-0.5+\epsilon ,0}\setminus N_{\delta }\left({\lambda }^{*}\right)$ with ${\lambda }^{*}\approx -0.1351$, and choose ${\lambda }_q\in {\mathcal{T}}_{c,r}$ such that

  $${\lambda }_q<\left\{\begin{array}{cc}\lambda ,\hfill & \mathrm{if}n\sim n_q,\hfill \\ 3/2,\hfill & \mathrm{if}n\gg n_q,\hfill \end{array}\right.$$

  where $N_{\delta }\left({\lambda }^{*}\right)$ excludes a $\delta$-neighborhood of ${\lambda }^{*}$ to avoid amplification of interpolation error.

  *

  Standard computational scenarios: utilize $\lambda ={\lambda }_q=0$, corresponding to shifted Chebyshev approximation, which offers a robust default balancing accuracy and efficiency for smooth functions.

• For $\alpha >1$, two strategies apply:

  (i)

  Classical transformation: use (2) with $\varphi \left(\alpha \right)=1/\alpha$, but anticipate algebraic quadrature convergence $\mathcal{O}\left(n_q^{-p}\right)$ with p decreasing as $\alpha$ increases (Figure 3f). Set ${\lambda }_q=0$ for optimal performance and use $n_q\sim {10}^3$–${10}^4$ to achieve high accuracy (Figure 4(c)).

  (ii)

  Generalized transformation: use (35) with $\varphi \left(\alpha \right)$ chosen via (A48) to enforce $C^r$ smoothness. As shown in Figure 5, with $r=12$ and $n_q=16$, this achieves LARE $<-14$ (ARE $<{10}^{-14}$) across $\alpha \in {\mathbf{\Omega }}_{1.5,9.5}$. Use ${\lambda }_q=0.5$ for best results (Figure 5e).

• Parameter sensitivity (Figure 3d,e, Figure 4b,c and Figure 5d,e):

  -

  $\lambda$: weakly influences the total error for $\alpha >1$ under the classical transformation (quadrature-dominated regime). For the $\varphi \left(\alpha \right)$-generalized transformation with sufficiently large r, the total error becomes ITE-dominated, and the choice of $\lambda$ is therefore more critical than in the classical, quadrature-dominated regime. Choosing $\lambda \in {\mathbf{\Omega }}_1$ ensures spectral stability.

  -

  ${\lambda }_q$: has a dominant influence on quadrature accuracy. For $\alpha >1$, the classical transformation performs best with ${\lambda }_q=0$. For the $\varphi \left(\alpha \right)$-generalized transformation, ${\lambda }_q=0.5$ is consistently optimal, and can improve LARE by 8–10 orders of magnitude compared with suboptimal choices.

  -

  n: once n exceeds the polynomial degree of a polynomial function f (or achieves the target ITE per Theorem 1), further refinement yields diminishing returns; focus instead on $n_q$.

  -

  $n_q$: governs quadrature error. For the generalized transformation with $\alpha >1$, $r=12$, $n_q=16$ often suffices for near machine precision (Figure 5f).

• For repeated evaluations with fixed $\alpha$, precompute the appropriate FSGIM: ${}^E\mathbf{Q}_n^{\alpha }$ for the classical transformation or (40) for the generalized transformation. For the classical FSGIM of [4], the cost analysis following [4] shows that the construction of ${}^E\mathbf{Q}_n^{\alpha }$ requires $O\left(n(n+n_q)\right)$ operations per single-point evaluation. After this offline stage, evaluating ${\mathfrak{I}}_t^{\alpha }f$ via (7) reduces to a single matrix–vector multiplication, giving an online cost $\mathcal{O}\left(n\right)$. For $M>1$ evaluation points, the classical FSGIM construction incurs $O\left(Mn(n+n_q)\right)$ operations. Once ${}^E\mathbf{Q}_n^{\alpha }$ is formed, the online evaluation of ${\mathfrak{I}}_{{\mathbf{z}}_M}^{\alpha }f$ requires only $\mathcal{O}\left(Mn\right)$ operations.
  For the weighted FSGIM defined row-wise by (40), Section 5 establishes that the per-row construction cost is $O\left(n(n+n_q)\right)$, identical in asymptotic order to the classical case. Hence, for $M+1$ evaluation points, the total offline cost to assemble the full weighted FSGIM is $O\left(Mn(n+n_q)\right)$. The subsequent online application to $f_{0:n}$ remains a dense matrix–vector product, incurring $\mathcal{O}\left(Mn\right)$ operations. Thus, both the classical and weighted FSGIMs share the same asymptotic offline/online computational structure, differing only in constant factors arising from the evaluation of the weight $y^{p\left(\alpha \right)}$ and the modified parameter $\varphi \left(\alpha \right)$.

8.2. Performance Comparison and Trade-Offs

• Classical transformation ($\varphi \left(\alpha \right)=1/\alpha$):

  -

  Advantage: Exact kernel cancellation in (3); optimal $\forall \alpha <1$.

  -

  Limitation: Algebraic quadrature convergence for $\alpha >1$ (Figure 3f and Figure 4d).

  -

  Best for: $\alpha <1$: either $1/\alpha \in {\mathbb{Z}}^{+}$ or moderate-accuracy regimes. In the latter case, although $1/\alpha \notin {\mathbb{Z}}^{+}$ implies only algebraic asymptotic quadrature convergence $\mathcal{O}\left(n_q^{-K-\theta }\right)$ with $K=\lfloor 1/\alpha \rfloor$, $\theta =1/\alpha -K$, the transformed integrand admits a fractional power series with exponentially decaying coefficients for real-analytic f, yielding dramatically faster pre-asymptotic decay that often appears near-exponential in practice (see Remark A9).

• $\varphi \left(\alpha \right)$-generalized transformation:

  -

  Advantage: Enforces $C^r$ smoothness via (A48), yielding near machine precision with $n_q=\mathcal{O}\left(10\right)$ (Figure 5f).

  -

  Limitation: Introduces weight $y^{p\left(\alpha \right)}$ with $p\left(\alpha \right)=\alpha \varphi \left(\alpha \right)-1$, necessitating tuned ${\lambda }_q$; moreover, for $0<\alpha <1$ with $1/\alpha \notin {\mathbb{Z}}^{+}$, achieving the asymptotic regime $\mathcal{O}\left(n_q^{-r}\right)$ requires $n_q$ to be sufficiently large (cf. Figure 12).

  -

  Best for: $\alpha >1$ or $0<\alpha <1:1/\alpha \notin {\mathbb{Z}}^{+}$ under high-accuracy demands (ARE $\sim {10}^{-15}$).

Appendix A.1. Orthogonality of SG Polynomials

The SG polynomials inherit the following well-known property.

Theorem A8 

(Orthogonality of SG Polynomials). The SG polynomials ${\widehat{G}}_n^{\left(\lambda \right)}$, defined for $t\in {\mathbf{\Omega }}_1$ and $\lambda >-1/2$, form an orthogonal basis on ${\mathbf{\Omega }}_1$ with respect to the weight function ${(t-t^2)}^{\lambda -1/2}$. Specifically, they satisfy the orthogonality relation:

where ${\delta }_{mn}$ is the Kronecker delta and ${\widehat{h}}_n^{\left(\lambda \right)}$ is the normalization constant given by:

$${\widehat{h}}_n^{\left(\lambda \right)}={\displaystyle \frac{\pi 2^{1-4\lambda }\Gamma (n+2\lambda )}{n!(n+\lambda ){\Gamma }^2\left(\lambda \right)}}.$$

Proof. 

The SG polynomials are defined in terms of standard Gegenbauer polynomials $G_n^{\left(\lambda \right)}$ through the transformation:

$${\widehat{G}}_n^{\left(\lambda \right)}=G_n^{\left(\lambda \right)}(2t-1),t\in {\mathbf{\Omega }}_1.$$

Standard Gegenbauer polynomials are orthogonal on $[-1,1]$ with respect to the weight function ${(1-x^2)}^{\lambda -1/2}$:

$${\mathcal{I}}_{-1,1}^{\left(t\right)}\left(G_m^{\left(\lambda \right)}{G}_n^{\left(\lambda \right)}{(1-x^2)}^{\lambda -1/2}\right)={\delta }_{mn}{h}_n^{\left(\lambda \right)},$$

where

$$h_n^{\left(\lambda \right)}={\displaystyle \frac{\pi 2^{1-2\lambda }\Gamma (n+2\lambda )}{n!(n+\lambda ){\Gamma }^2\left(\lambda \right)}}.$$

Substituting $x=2t-1$ into the standard orthogonality relation:

$$\begin{array}{cc}\hfill & {\mathcal{I}}_{-1,1}^{\left(t\right)}\left(G_m^{\left(\lambda \right)}{G}_n^{\left(\lambda \right)}{(1-x^2)}^{\lambda -1/2}\right)\hfill \\ \hfill & ={\mathcal{I}}_1^{\left(t\right)}\left(G_m^{\left(\lambda \right)}2t-1G_n^{\left(\lambda \right)}4(t-t^2){]}^{\lambda -1/2}·2\right)\hfill \\ \hfill & =2·4^{\lambda -1/2}{\mathcal{I}}_1^{\left(t\right)}\left({\widehat{G}}_m^{\left(\lambda \right)}{\widehat{G}}_n^{\left(\lambda \right)}{(t-t^2)}^{\lambda -1/2}\right).\hfill \end{array}$$

Therefore:

$${\mathcal{I}}_1^{\left(t\right)}\left({\widehat{G}}_m^{\left(\lambda \right)}{\widehat{G}}_n^{\left(\lambda \right)}{(t-t^2)}^{\lambda -1/2}\right)={\displaystyle \frac{{\delta }_{mn}{h}_n^{\left(\lambda \right)}}{2^{2\lambda }}}={\delta }_{mn}{\widehat{h}}_n^{\left(\lambda \right)},$$

where ${\widehat{h}}_n^{\left(\lambda \right)}=h_n^{\left(\lambda \right)}/2^{2\lambda }$. This completes the proof of orthogonality. The completeness of the polynomial basis follows from standard results in approximation theory. □

Appendix A.2. Bijectivity of the RLFI Transformation

Theorem A9 

(Bijectivity of the RLFI Transformation). For any fixed $t>0$ and $\alpha >0$, the mapping $\tau :{\mathbf{\Omega }}_1^{\circ }\to {\mathbf{\Omega }}_t^{\circ }$ defined by $\tau \left(y\right)=t(1-y^{1/\alpha })$ is bijective. That is, it is both injective and surjective, establishing a perfect correspondence between the intervals ${\mathbf{\Omega }}_1^{\circ }$ and ${\mathbf{\Omega }}_t^{\circ }$.

Proof. 

A function is bijective if and only if it is both injective and surjective.

1.

Injectivity (One-to-One): Assume $\tau \left(y_1\right)=\tau \left(y_2\right)$ for $y_1,y_2\in {\mathbf{\Omega }}_1^{\circ }$. Then:

$$t(1-y_1^{1/\alpha })=t(1-y_2^{1/\alpha }).$$

Since $t>0$, we can divide both sides by t:

$$1-y_1^{1/\alpha }=1-y_2^{1/\alpha }\Rightarrow y_1^{1/\alpha }=y_2^{1/\alpha }\Rightarrow y_1=y_2.$$

Thus, $\tau$ is injective.

2.

Surjectivity (Onto): ${\forall }_a$${\tau }_0\in {\mathbf{\Omega }}_t^{\circ }$, we need to find $y_0\in {\mathbf{\Omega }}_1^{\circ }$ such that $\tau \left(y_0\right)={\tau }_0$. Solving:

$${\tau }_0=t(1-y_0^{1/\alpha })\Rightarrow 1-y_0^{1/\alpha }={\displaystyle \frac{{\tau }_0}{t}}\Rightarrow y_0^{1/\alpha }=1-{\displaystyle \frac{{\tau }_0}{t}}.$$

Since ${\tau }_0\in {\mathbf{\Omega }}_t^{\circ }$, we have ${\displaystyle \frac{{\tau }_0}{t}}\in {\mathbf{\Omega }}_1^{\circ }$, so $1-{\displaystyle \frac{{\tau }_0}{t}}\in {\mathbf{\Omega }}_1^{\circ }$. Raising both sides to the power $\alpha$:

$$y_0={\left(1-{\displaystyle \frac{{\tau }_0}{t}}\right)}^{\alpha }\in {\mathbf{\Omega }}_1^{\circ }.$$

Thus, $\forall {\tau }_0\in {\mathbf{\Omega }}_t^{\circ },\exists y_0={\left(1-{\displaystyle \frac{{\tau }_0}{t}}\right)}^{\alpha }\in {\mathbf{\Omega }}_1^{\circ }:\tau \left(y_0\right)={\tau }_0$. Therefore, $\tau$ is surjective.

3.

Explicit Inverse: The inverse function $y:{\mathbf{\Omega }}_t^{\circ }\to {\mathbf{\Omega }}_1^{\circ }$ is explicitly given by:

$$y\left(\tau \right)={\left(1-{\displaystyle \frac{\tau }{t}}\right)}^{\alpha }.$$

This confirms the bijection since a function is bijective if and only if it has an inverse.

4.

Boundary Behavior:

• By definition: As $y\to 0^{+}$, $\tau \to t^{-}$; as $y\to 1^{-}$, $\tau \to 0^{+}$.
• By injectivity and the inverse function: As $\tau \to 0^{+}$, $y\to 1^{-}$; as $\tau \to t^{-}$, $y\to 0^{+}$.

Thus, the mapping preserves the open intervals: ${\mathbf{\Omega }}_1^{\circ }\leftrightarrow {\mathbf{\Omega }}_t^{\circ }$.

□

Corollary A2 

(Diffeomorphism Property). The transformation (2) constitutes a $C^{\infty }$-diffeomorphism between ${\mathbf{\Omega }}_1^{\circ }$ and ${\mathbf{\Omega }}_t^{\circ }$. Specifically, it has non-vanishing Jacobian

$${\displaystyle \frac{d\tau }{dy}}=-{\displaystyle \frac{t}{\alpha }}{y}^{1/\alpha -1},$$

(A46)

and both $\tau \left(y\right)$ and its inverse

$$y\left(\tau \right)={\left(1-{\displaystyle \frac{\tau }{t}}\right)}^{\alpha },$$

(A47)

are $C^{\infty }$ on their respective domains.

Proof. 

For fixed $t>0$ and $\alpha >0$, the map $y\mapsto y^{1/\alpha }$ is well-defined and real-analytic on ${\mathbf{\Omega }}_1^{\circ }$. Hence (2) is real-analytic, in particular $\tau \in C^{\infty }\left({\mathbf{\Omega }}_1^{\circ }\right)$. Differentiating (2) with respect to y gives (A46). Since $t>0$, $\alpha >0$, and $y\in {\mathbf{\Omega }}_1^{\circ }$, we have $y^{1/\alpha -1}>0$, so

$${\displaystyle \frac{d\tau }{dy}}\left(y\right)\ne 0,\forall y\in {\mathbf{\Omega }}_1^{\circ }.$$

By Theorem A9, the mapping $\tau :{\mathbf{\Omega }}_1^{\circ }\to {\mathbf{\Omega }}_t^{\circ }$ is bijective. Because $\tau$ is $C^{\infty }$ with non-vanishing derivative on ${\mathbf{\Omega }}_1^{\circ }$, the inverse function theorem implies that the inverse map (A47) is also $C^{\infty }$ on ${\mathbf{\Omega }}_t^{\circ }$. Thus $\tau$ is a $C^{\infty }$-diffeomorphism between ${\mathbf{\Omega }}_1^{\circ }$ and ${\mathbf{\Omega }}_t^{\circ }$. □

Appendix A.3. Smoothness Characterization of GBFA-Transformed Integrands and Its Implications for Convergence Rates

Theorem A10 

(Exact Smoothness of the Transformed Integrand under the Classical GBFA Map). Let $0<\alpha <1$, $\beta =1/\alpha >1$, $t\in {\mathbf{\Omega }}_1$, and assume f is real-analytic on a neighborhood of ${\mathbf{\Omega }}_t$. Define the transformed integrand

$$g\left(y\right)=f\left(t(1-y^{\beta })\right),y\in {\mathbf{\Omega }}_1.$$

Set $K=\lfloor \beta \rfloor$ and $\theta =\beta -K\in [0,1)$. Then:

(i) 

If $\theta >0$ (i.e., $\beta \notin {\mathbb{Z}}^{+}$), then $g\in C^{K,\theta }\left({\mathbf{\Omega }}_1\right)$ and, generically, $g\notin C^{K+1}\left({\mathbf{\Omega }}_1\right)$.

(ii) 

If $\theta =0$ (i.e., $\beta \in {\mathbb{Z}}^{+}$), then g is real-analytic on ${\mathbf{\Omega }}_1$.

Proof. 

Define $\varphi \left(y\right)=t(1-y^{\beta })$. Since f is real-analytic on a neighborhood of ${\mathbf{\Omega }}_t$ and $\varphi \left(y\right)\in {\mathbf{\Omega }}_t$ for all $y\in {\mathbf{\Omega }}_1$, the composition $g=f\circ \varphi$ is real-analytic on ${\mathbf{\Omega }}_1^{\circ }$ by the chain rule.

To analyze regularity at $y=0$, first note that for any $m\in {\mathbb{Z}}^{+}$ and any $\beta >1$,

$${\displaystyle \frac{d^m}{d{y}^m}}{y}^{\beta }=c_{\beta ,m}{y}^{\beta -m},$$

for some constant $c_{\beta ,m}$ (and $c_{\beta ,m}\ne 0$ whenever $\beta \notin {\mathbb{Z}}^{+}$ or $m\le \beta$). Thus

$${\varphi }^{\left(m\right)}\left(y\right)=-t{c}_{\beta ,m}{y}^{\beta -m}.$$

By Faà di Bruno’s formula [17],

$$g^{\left(m\right)}\left(y\right)=\sum _{k=1}^m{f}^{\left(k\right)}\left(\varphi \left(y\right)\right)B_{m,k}\left({\varphi }^{\prime }\left(y\right),\dots ,{\varphi }^{(m-k+1)}\left(y\right)\right),$$

where each monomial in $B_{m,k}$ is a product of powers of ${\varphi }^{\left(j\right)}\left(y\right)$ and scales as $y^{\beta k-m}$. The most singular contribution arises when $k=1$, yielding asymptotically

$$g^{\left(m\right)}\left(y\right)\sim f^{\prime }\left(\varphi \left(y\right)\right)y^{\beta -m},y\to 0^{+}.$$

Since $f^{\prime }\left(\varphi \left(y\right)\right)\to f^{\prime }\left(t\right)$ as $y\to 0^{+}$ and $f^{\prime }\left(t\right)\ne 0$ generically, $g^{\left(m\right)}$ remains bounded near $y=0$ if and only if $\beta -m\ge 0$, i.e., $m\le \beta$.

We now distinguish two cases.

Case (i): $\beta \notin {\mathbb{Z}}^{+}$. Then $K=\lfloor \beta \rfloor <\beta <K+1$, so $g^{\left(m\right)}$ is continuous on ${\mathbf{\Omega }}_1\forall m\le K$, while for $m=K+1$ we have

$$g^{(K+1)}\left(y\right)\sim C{y}^{\theta -1},\theta =\beta -K\in {\mathbf{\Omega }}_1^{\circ },$$

which diverges as $y\to 0^{+}$. Moreover,

$$g^{\left(K\right)}\left(y\right)\sim C{y}^{\theta },y\to 0^{+},$$

so $g^{\left(K\right)}$ is Hölder-continuous of order $\theta$ at $y=0$. Hence $g\in C^{K,\theta }\left({\mathbf{\Omega }}_1\right)$ and, generically, $g\notin C^{K+1}\left({\mathbf{\Omega }}_1\right)$.

Case (ii): $\beta =m\in {\mathbb{Z}}^{+}$. In this case $y^{\beta }=y^m$ is a polynomial, so $\varphi \left(y\right)=t(1-y^m)$ is also a polynomial and therefore real-analytic on ${\mathbf{\Omega }}_1$. Since f is real-analytic on a neighborhood of ${\mathbf{\Omega }}_t$, the composition $g=f\circ \varphi$ is real-analytic on ${\mathbf{\Omega }}_1$. In particular, $g\in C^{\infty }\left({\mathbf{\Omega }}_1\right)$.

This completes the proof. □

The preceding analysis leads directly to the following corollary, formalizing the resulting regularity and quadrature behavior.

Corollary A3. 

Let $f\in C^{\infty }\left({\mathbf{\Omega }}_1\right)$ and define $g\left(y\right)=f\left(t(1-y^{1/\alpha })\right)$. Then the following hold:

(i)

If $1/\alpha =m\in {\mathbb{Z}}^{+}$, then g is real-analytic on ${\mathbf{\Omega }}_1$, and the $(n_q,{\lambda }_q)$-GBFA quadrature applied to the SG–transformed integrals in (6) converges super-exponentially in $n_q$. If, in addition, f is a polynomial of degree d, then g is a polynomial of degree at most $dm$, and the quadrature is exact for all $n_q\ge dm$. Consequently, the total GBFA error reduces to the ITE, which vanishes for all $n\ge d$, yielding exact evaluation of ${\mathfrak{I}}_t^{\alpha }f$ for every $n\ge d$ and $n_q\ge d/\alpha$.

(ii)

If $1/\alpha \notin {\mathbb{Z}}^{+}$, then $g\in C^{K,\theta }\left({\mathbf{\Omega }}_1\right)$ with $K=\lfloor 1/\alpha \rfloor$ and $\theta =1/\alpha -K\in {\mathbf{\Omega }}_1^{\circ }$, and generically $g\notin C^{K+1}\left({\mathbf{\Omega }}_1\right)$; cf. Theorem A10(i). In this case, the quadrature error decays algebraically as $\mathcal{O}\left(n_q^{-K-\theta }\right)$ under Gauss-type rules.

Remark A9. 

The algebraic rate $\mathcal{O}\left(n_q^{-K-\theta }\right)$ in Corollary A3(ii) represents the worst-case convergence for arbitrary $f\in C^{\infty }\left({\mathbf{\Omega }}_1\right)$. However, if f is real-analytic on a neighborhood of ${\mathbf{\Omega }}_t$, then the transformed integrand

$$g\left(y\right)=f\left(t(1-y^{1/\alpha })\right),y\in {\mathbf{\Omega }}_1,$$

admits a convergent fractional power series expansion near $y=0$ of the form

$$g\left(y\right)=\sum _{m=0}^{\infty }{a}_m{y}^{m/\alpha },a_m=c_m{(-t)}^m,$$

where $f\left(z\right)={\sum }_{m=0}^{\infty }{c}_m{(z-t)}^m$ with $|c_m|\le C{B}^m$ for some $C,B>0$. Consequently, the coefficients satisfy $|a_m|\le \tilde{C}{\tilde{B}}^m$ for $\tilde{C}=C$ and $\tilde{B}=\left|t\right|B$, exhibiting exponential-type decay. This rapid decay implies that the tail of the series beyond indices $m\gtrsim n_q$ is extremely small. As a result, Gauss-type quadrature rules–though asymptotically limited to an algebraic rate by the endpoint Hölder regularity–often exhibit dramatically faster pre-asymptotic convergence in practice. The observed convergence can thus appear nearly exponential or high-order algebraic over typical computational ranges ($n_q\lesssim {10}^2$–${10}^3$), even though the strict asymptotic rate $\mathcal{O}\left(n_q^{-K-\theta }\right)$ governs only the limit $n_q\to \infty$.

Appendix A.4. Diffeomorphism Property of the ϕ(α)-Mapping

We now formalize the smoothness and invertibility of the generalized RLFI transformation (35) introduced in Section 4. The result is the exact analogue of Corollary A2 for the mapping (2).

Theorem A11 

($C^{\infty }$-diffeomorphism of the $\varphi \left(\alpha \right)$-mapping). Let $t\in {\mathbb{R}}^{+}$, $\alpha >0$, and let $\varphi \left(\alpha \right)\in {\mathbb{R}}^{+}$ be fixed. Consider the mapping

$$\tau :{\mathbf{\Omega }}_1^{\circ }\to {\mathbf{\Omega }}_t^{\circ },\tau \left(y\right)=t\left(1-y^{\varphi \left(\alpha \right)}\right),$$

as defined in (35). Then τ is a $C^{\infty }$-diffeomorphism between ${\mathbf{\Omega }}_1^{\circ }$ and ${\mathbf{\Omega }}_t^{\circ }$; that is,

(i) 

τ is bijective,

(ii) 

$\tau \in C^{\infty }\left({\mathbf{\Omega }}_1^{\circ }\right)$ with

$${\displaystyle \frac{d\tau }{dy}}\left(y\right)=-t\varphi \left(\alpha \right)y^{\varphi \left(\alpha \right)-1}\ne 0,\forall y\in {\mathbf{\Omega }}_1^{\circ },$$

(iii) 

its inverse $y:{\mathbf{\Omega }}_t^{\circ }\to {\mathbf{\Omega }}_1^{\circ }$,

$$y\left(\tau \right)={\left(1-{\displaystyle \frac{\tau }{t}}\right)}^{1/\varphi \left(\alpha \right)},$$

belongs to $C^{\infty }\left({\mathbf{\Omega }}_t^{\circ }\right)$.

Proof. 

To prove the smoothness and non-vanishing Jacobian, notice that ${\forall }_a\alpha >0$ and $\varphi \left(\alpha \right)>0$, the map $y\mapsto y^{\varphi \left(\alpha \right)}$ is real-analytic on ${\mathbf{\Omega }}_1^{\circ }$; thus the map (35) is real-analytic on ${\mathbf{\Omega }}_1^{\circ }$, in particular $\tau \in C^{\infty }\left({\mathbf{\Omega }}_1^{\circ }\right)$. Differentiating (35) with respect to y yields

$${\displaystyle \frac{d\tau }{dy}}\left(y\right)=-t\varphi \left(\alpha \right)y^{\varphi \left(\alpha \right)-1},y\in {\mathbf{\Omega }}_1^{\circ }.$$

Since $t>0$, $\varphi \left(\alpha \right)>0$, and $y^{\varphi \left(\alpha \right)-1}>0\forall y\in {\mathbf{\Omega }}_1^{\circ }$, it follows that

$${\displaystyle \frac{d\tau }{dy}}\left(y\right)\ne 0,\forall y\in {\mathbf{\Omega }}_1^{\circ }.$$

The bijectivity proof argument parallels that of Theorem A9, but now for (35). In particular, to verify the map’s injectivity, assume $\tau \left(y_1\right)=\tau \left(y_2\right)$ for some $y_1,y_2\in {\mathbf{\Omega }}_1^{\circ }$. Then

$$t\left(1-y_1^{\varphi \left(\alpha \right)}\right)=t\left(1-y_2^{\varphi \left(\alpha \right)}\right).$$

Dividing by $t>0$ gives

$$1-y_1^{\varphi \left(\alpha \right)}=1-y_2^{\varphi \left(\alpha \right)}\Rightarrow y_1^{\varphi \left(\alpha \right)}=y_2^{\varphi \left(\alpha \right)}\Rightarrow y_1=y_2,$$

since the map $y\mapsto y^{\varphi \left(\alpha \right)}$ is strictly increasing on ${\mathbf{\Omega }}_1^{\circ }$. Hence $\tau$ is injective. To prove surjectivity, let ${\tau }_0\in {\mathbf{\Omega }}_t^{\circ }$ be arbitrary. We seek $y_0\in {\mathbf{\Omega }}_1^{\circ }$ such that $\tau \left(y_0\right)={\tau }_0$. Solving

$${\tau }_0=t\left(1-y_0^{\varphi \left(\alpha \right)}\right)$$

for $y_0$ gives

$$y_0={\left(1-{\displaystyle \frac{{\tau }_0}{t}}\right)}^{1/\varphi \left(\alpha \right)}.$$

Since ${\tau }_0\in {\mathbf{\Omega }}_t^{\circ }$, we have ${\displaystyle \frac{{\tau }_0}{t}}\in {\mathbf{\Omega }}_1^{\circ }$, hence

$$1-{\displaystyle \frac{{\tau }_0}{t}}\in {\mathbf{\Omega }}_1^{\circ }\Rightarrow y_0\in {\mathbf{\Omega }}_1^{\circ }.$$

Thus, $\forall {\tau }_0\in {\mathbf{\Omega }}_t^{\circ },\exists y_0\in {\mathbf{\Omega }}_1^{\circ }$ with $\tau \left(y_0\right)={\tau }_0$, and $\tau$ is surjective. Combining injectivity and surjectivity, $\tau$ is bijective with inverse

$$y\left(\tau \right)={\left(1-{\displaystyle \frac{\tau }{t}}\right)}^{1/\varphi \left(\alpha \right)},\tau \in {\mathbf{\Omega }}_t^{\circ }.$$

The smoothness of the inverse follows by realizing that the map

$$\tau \mapsto 1-{\displaystyle \frac{\tau }{t}}$$

is real-analytic on ${\mathbf{\Omega }}_t^{\circ }$, and its image lies in ${\mathbf{\Omega }}_1^{\circ }$. The map

$$y\mapsto y^{1/\varphi \left(\alpha \right)},$$

with $\varphi \left(\alpha \right)>0$ fixed, is real-analytic on ${\mathbf{\Omega }}_1^{\circ }$. Therefore, their composition

$$y\left(\tau \right)={\left(1-{\displaystyle \frac{\tau }{t}}\right)}^{1/\varphi \left(\alpha \right)}$$

is real-analytic on ${\mathbf{\Omega }}_t^{\circ }$, in particular $y\in C^{\infty }\left({\mathbf{\Omega }}_t^{\circ }\right)$. In summary, we have shown that $\tau :{\mathbf{\Omega }}_1^{\circ }\to {\mathbf{\Omega }}_t^{\circ }$ is bijective, $\tau \in C^{\infty }\left({\mathbf{\Omega }}_1^{\circ }\right)$ with $\frac{d\tau }{dy}\ne 0$ on ${\mathbf{\Omega }}_1^{\circ }$, and its inverse y is $C^{\infty }$ on ${\mathbf{\Omega }}_t^{\circ }$. Hence, by definition, $\tau$ is a $C^{\infty }$-diffeomorphism between ${\mathbf{\Omega }}_1^{\circ }$ and ${\mathbf{\Omega }}_t^{\circ }$. □

Appendix A.5. Smoothness of g via Faà di Bruno Exponent Analysis

Theorem A12 

(Sufficient Smoothness Condition for $g\left(y\right)=y^{\alpha \varphi \left(\alpha \right)-1}f\left(t(1-y^{\varphi \left(\alpha \right)})\right)$). Let $\alpha >0$, $t\in {\mathbf{\Omega }}_1$, and $f\in C^{\infty }\left({\mathbf{\Omega }}_1\right)$. Define

$$g\left(y\right)=y^{p\left(\alpha \right)}f\left(t(1-y^{\varphi \left(\alpha \right)})\right),p\left(\alpha \right)=\alpha \varphi \left(\alpha \right)-1,$$

for $y\in {\mathbf{\Omega }}_1$. Then, for any $r\in \mathbb{N}$, a sufficient condition for $g\in C^r\left({\mathbf{\Omega }}_1\right)$ is

$$\varphi \left(\alpha \right)\ge {\displaystyle \frac{r+1}{\alpha }}.$$

(A48)

Proof. 

Set $\varphi =\varphi \left(\alpha \right)$ and $p=p\left(\alpha \right)=\alpha \varphi -1$ for brevity, and define $\psi \left(y\right)=t(1-y^{\varphi })$. Then $g\left(y\right)=y^{p}f\left(\psi \left(y\right)\right)$. Since $f\in C^{\infty }\left({\mathbf{\Omega }}_1\right)$ and $\psi \in C^{\infty }\left({\mathbf{\Omega }}_1^{\circ }\right)$, g is smooth on ${\mathbf{\Omega }}_1^{\circ }$. We analyze the behavior of $g^{\left(r\right)}\left(y\right)$ as $y\to 0^{+}$. By the Leibniz rule,

$$g^{\left(r\right)}\left(y\right)=\sum _{\ell =0}^r\left({\displaystyle \genfrac{}{}{0pt}{}{r}{\ell }}\right){\displaystyle \frac{d^{\ell }}{d{y}^{\ell }}}{y}^p·{\displaystyle \frac{d^{r-\ell }}{d{y}^{r-\ell }}}f\left(\psi \left(y\right)\right).$$

Step 1: The $\ell =r$ term. Repeated differentiation of $y^p$ shows that

$${\displaystyle \frac{d^r}{d{y}^r}}{y}^p=c_{p,r}{y}^{p-r},$$

for some finite constant $c_{p,r}$ (possibly zero only if $p\in {\mathbb{Z}}_0^{+}$ and $r>p$). Since $\psi \left(y\right)\to t$ as $y\to 0^{+}$ and f is smooth on ${\mathbf{\Omega }}_1$, we have

$$f\left(\psi \left(y\right)\right)=f\left(t\right)+o\left(1\right),y\to 0^{+}.$$

Thus, unless $c_{p,r}f\left(t\right)=0$, the $\ell =r$ term behaves asymptotically as

$${\displaystyle \frac{d^r}{d{y}^r}}{y}^p·f\left(\psi \left(y\right)\right)\sim C{y}^{p-r},y\to 0^{+},$$

for some constant $C\ne 0$. In particular, its exponent is $p-r=\alpha \varphi -1-r$. If $c_{p,r}f\left(t\right)=0$, then this contribution either vanishes identically or is replaced by a term with a strictly larger exponent, which can only improve regularity at $y=0$.

Step 2: The terms with $\ell <r$. Each such term contains at least one derivative of the composition $f\circ \psi$. By Faà di Bruno’s formula, we can write

$${\displaystyle \frac{d^{r-\ell }}{d{y}^{r-\ell }}}f\left(\psi \left(y\right)\right)=\sum _{k=1}^{r-\ell }{f}^{\left(k\right)}\left(\psi \left(y\right)\right)B_{r-\ell ,k}\left({\psi }^{\prime }\left(y\right),\dots ,{\psi }^{(r-\ell -k+1)}\left(y\right)\right),$$

where each monomial in the Bell polynomial $B_{r-\ell ,k}$ is of the form

$$f^{\left(k\right)}\left(\psi \left(y\right)\right)\prod _{j=1}^{r-\ell }{\left({\psi }^{\left(j\right)}\left(y\right)\right)}^{u_j},$$

with nonnegative integers $u_1,\dots ,u_{r-\ell }$ satisfying

$$\sum _{j=1}^{r-\ell }j{u}_j=r-\ell ,\sum _{j=1}^{r-\ell }{u}_j=k.$$

Since ${\psi }^{\left(j\right)}\left(y\right)=-t\varphi (\varphi -1)\dots (\varphi -j+1)y^{\varphi -j}\sim C_j{y}^{\varphi -j}$ as $y\to 0^{+}$, the product scales as

$$\prod _{j=1}^{r-\ell }{\left({\psi }^{\left(j\right)}\left(y\right)\right)}^{u_j}\sim \tilde{C}{y}^{{\sum }_{j=1}^{r-\ell }(\varphi -j)u_j}=\tilde{C}{y}^{\varphi k-(r-\ell )}.$$

Multiplying by $\frac{d^{\ell }}{d{y}^{\ell }}{y}^p\sim c_{p,\ell }{y}^{p-\ell }$ gives a total factor of

$$y^{p-\ell }·y^{\varphi k-(r-\ell )}=y^{p+\varphi k-r}.$$

Since $k\ge 1$ and $\varphi >0$, we have

$$p+\varphi k-r\ge p+\varphi -r>p-r.$$

Thus every term with $\ell <r$ carries an exponent strictly larger than $p-r$, and therefore decays faster as $y\to 0^{+}$ than the generic $\ell =r$ contribution.

Step 3: Boundedness of $g^{\left(r\right)}$ at $y=0$. Combining the two steps, we see that, in the generic case where $c_{p,r}f\left(t\right)\ne 0$, the smallest exponent among all contributions to $g^{\left(r\right)}\left(y\right)$ is $p-r$. Hence the worst-case behaviour near $y=0$ is governed by

$$g^{\left(r\right)}\left(y\right)=\mathcal{O}\left(y^{p-r}\right)=\mathcal{O}\left(y^{\alpha \varphi -1-r}\right),y\to 0^{+}.$$

Therefore $g^{\left(r\right)}$ remains bounded at $y=0$ whenever

$$\alpha \varphi -1-r\ge 0,$$

which gives the smoothness condition (A48). If cancellations occur (e.g., $c_{p,r}f\left(t\right)=0$), the first nonzero term in $g^{\left(r\right)}$ has exponent strictly larger than $p-r$, so the resulting regularity at $y=0$ is even higher. Thus (A48) is a sufficient condition for $g\in C^r\left({\mathbf{\Omega }}_1\right)$. □

Appendix A.6. Gevrey Regularity of the Exact Solution (45)

We establish the precise smoothness class of the exact solution (45) of (44). Although u is real-analytic on ${\mathbf{\Omega }}_1$, its derivatives grow factorially. We prove that u belongs to the Gevrey-1 class, which is equivalent to the class of real-analytic functions, but does not satisfy the derivative bound (9) required by Theorem 1.

Theorem A13 

(Gevrey-1 regularity of $u\left(x\right)=\sqrt{\pi }{(1+x)}^{-3/2}$). The function u defined by (45) is real-analytic on ${\mathbf{\Omega }}_1$ and belongs to the Gevrey-1 class. Specifically, its derivatives satisfy

$$\parallel u^{\left(k\right)}{\parallel }_{L^{\infty }\left({\mathbf{\Omega }}_1\right)}\le C{B}^{k}k!,\forall k\in {\mathbb{Z}}_0^{+},$$

(A49)

for some constants $C,B>0$. Consequently, u does not satisfy the derivative growth condition (9).

Proof. 

The function $u\left(x\right)=\sqrt{\pi }{(1+x)}^{-3/2}$ is the restriction to ${\mathbf{\Omega }}_1$ of the complex function $z\mapsto \sqrt{\pi }{(1+z)}^{-3/2}$, which is analytic on $\mathbb{C}\setminus \{z\in \mathbb{R}:z\le -1\}$. Since the singularity at $z=-1$ is at a distance of 1 from ${\mathbf{\Omega }}_1$, u is real-analytic on ${\mathbf{\Omega }}_1$.

We now compute the kth derivative. For $\nu =3/2$, the kth derivative of ${(1+x)}^{-\nu }$ is

$${\displaystyle \frac{d^k}{d{x}^k}}{(1+x)}^{-\nu }={(-1)}^k{\displaystyle \frac{\Gamma (k+\nu )}{\Gamma \left(\nu \right)}}{(1+x)}^{-(k+\nu )}.$$

(A50)

Using $\Gamma (3/2)=\sqrt{\pi }/2$ and $u\left(x\right)=\sqrt{\pi }{(1+x)}^{-3/2}$, we obtain

$$|u^{\left(k\right)}\left(x\right)|\le 2\Gamma \left(k+{\displaystyle \frac{3}{2}}\right),\forall x\in {\mathbf{\Omega }}_1.$$

Applying the property of the Gamma function [22, Eq. (6.1.12)],

$$\Gamma (k+3/2)={\displaystyle \frac{(2k+1)!!}{2^{k+1}}}\sqrt{\pi },$$

and the bound for the double factorial $(2k+1)!!\le (2k+1)!$, we can use Stirling’s approximation to find that there exist constants $C_1,B_1>0$ such that

$$\Gamma \left(k+{\displaystyle \frac{3}{2}}\right)\le C_1{B}_1^{k}k!,\forall k\in {\mathbb{Z}}_0^{+}.$$

Thus, for $C=2C_1$ and $B=B_1$, the bound (A49) holds. Consequently, u fails to meet the hypothesis of Theorem 1. □