## 4 Error Analysis
 \tag{47}\]  
_with \(D^{\lambda}>1\) being a \(\lambda\)-dependent constant, and \(\theta_{\alpha,\lambda},\mathrm{Y}_{\sigma^{\lambda}}(n),\gamma_{n_{q}}^{ \lambda}\), and \(\Theta_{\lambda_{q}}\) are constants with the definitions and properties outlined in Theorems 2 and 4, and Lemmas A1 and A3, and \(\omega^{\text{upper},\pm}\) are as defined by [38] (Formulas (B.4) and (B.15))._  
Proof. : The total truncation error combines the interpolation error from Theorem 2 and the accumulated quadrature errors from Theorem 4 for \(j\in\mathbb{J}_{n}^{+}\):  
\[\begin{array}{l}{}^{\alpha}\xi_{n,n_{q}}^{\lambda,\lambda_{q}}(t,\xi,\eta) ={}^{\alpha}\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! The proof is completed by applying Theorems 2 and 4 on Formula (48), noting that \(\max_{j\in N_{n_{q}+1n}}\left|\Delta_{j,n_{q}}^{\lambda_{q}}(\eta)\right|\) occurs at \(j=n\). The total truncation error bound presented in Theorem 5 reveals several important insights about the convergence behavior of the RLFI approximation: (i) The total error consists of the interpolation error term \({}^{\alpha}\mathcal{T}_{n}^{\lambda}(t,\xi)\) and the accumulated quadrature error term. The interpolation error term decays at a super-exponential rate \(\,\forall t\,n\) due to the factor \(\left(\frac{\epsilon}{4n}\right)^{n}\). The quadrature error either vanishes when \(n_{q}\geq n\), decays exponentially when \(n_{q}\sim n\) under Condition (41), or typically diverges when \(n\gg n_{q}\) unless Condition (43) is fulfilled. Therefore, the quadrature nodes should scale appropriately with the interpolation mesh size in practice. (ii) The interpolation error bound tightens as \(\alpha\) increases, due to the \(\theta_{\alpha,\lambda}\) factor's decay. The \(\lambda\) parameter shows a unimodal influence on the interpolation error, with potentially maximum error size at about \(\lambda^{*}\approx-0.1351\). The \(\lambda_{q}\) parameter should generally be chosen smaller than \(\lambda\) to accelerate the convergence of the quadrature error when \(n_{q}\sim n\). This analysis suggests that the proposed method achieves spectral convergence when the parameters are chosen appropriately. The quadrature precision should be selected to maintain balance between the two error components based on the desired accuracy and computational constraints. **Remark 2**: _Theorem 5 provides a theoretical guarantee for the accuracy of the GBFA method when computing the RLFI. In simple terms, it shows that the error in our method's approximations decreases extremely rapidly as we increase the number of computational points, particularly for smooth functions. This rapid error reduction, often referred to as super-exponential convergence, means that with just a few additional points, the GBFA method can achieve results that are nearly as accurate as the computer's maximum precision allows. The theorem considers key parameters of the GBFA method: the polynomial degree \(n\), the quadrature degree \(n_{q}\), and the Gegenbauer parameters \(\lambda\) and \(\lambda_{q}\). It predicts that when these parameters are chosen wisely, the error shrinks faster than exponentially, as we demonstrate later in our numerical experiments. This is especially powerful for modeling complex systems with memory, such as those in viscoelasticity or anomalous diffusion, where high precision is critical. While the mathematical details of the theorem are complex, its core message is that the GBFA method is highly reliable and efficient, producing accurate results with minimal computational effort for a wide range of problems._