An Exact-Moment Local Legendre Frame Method with Block Convolution for Caputo Fractional Differentiation

Abstract: We propose a local Legendre frame method for the accurate computation of Caputo fractional derivatives of order (0<α<1). On each local subinterval, the function is represented by a restricted Legendre frame obtained from scaled Legendre polynomials on an extended interval. The local coefficients are computed from equispaced samples by an exponentially weighted GTSVD regularization. The Caputo derivative is then evaluated by applying the weakly singular fractional integral to the derivatives of the local frame basis functions. Since these derivatives are polynomials, the corresponding Caputo weights can be written in terms of finite weighted moments, so that the singular kernel is treated analytically rather than by a low-order quadrature rule. For uniform partitions, the history weights have a block-dependent structure and can be reused efficiently. The error analysis separates the local frame reconstruction from the Caputo integration. In particular, the Caputo error is bounded by the derivative reconstruction error, while the latter is obtained from the (L^2) reconstruction error and a weighted smoothness bound of the GTSVD approximation through an interpolation argument. For analytic local functions with exponential coefficient decay, this leads to exponential-type convergence of the derivative and hence of the Caputo approximation. Numerical experiments confirm the accuracy of the exact moment weights, the effectiveness of the local weighted reconstruction, and the efficiency of the block implementation.

• The paper develops a regularized local Legendre frame method that computes Caputo derivatives from equispaced samples using exact weak-singularity moments and stable GTSVD coefficient recovery.
• The method achieves roughly 10^-13 to 10^-12 errors on analytic test functions, outperforming L1 discretization by 8–10 orders of magnitude and reaching 10^-12 accuracy in a fractional differential-equation solver.
• Uniform blocks create lag-dependent history weights, reducing storage to O(MN_b) and enabling FFT-based application in O(N_bM log M), while limitations remain for noisy, nonsmooth, and initially singular solutions.

Overview

The paper develops a numerical method for computing the left-sided Caputo fractional derivative of order $0<\alpha<1$ from equispaced samples, combining three ingredients: a restricted scaled Legendre frame representation on each subinterval, an exponentially weighted truncated generalized singular value decomposition (GTSVD) for stable coefficient recovery, and exact analytic moment formulas for the weakly singular Caputo kernel. The method is designed for uniform partitions, where the history weights acquire a lag-dependent block-convolution structure that permits FFT-based accumulation with nearly linear complexity in the number of samples.

Method construction

On each subinterval $I_k$ of a partition of $[0,b]$, the local function is represented by the restricted frame $\{P_j^{(T)}|_{[-1,1]}\}_{j=0}^{N_b-1}$, where $P_j^{(T)}(t)=T^{-1/2}p_j(t/T)$ are Legendre polynomials orthonormal on the extended interval $[-T,T]$ with $T>1$. The restriction is essential: the system is a redundant polynomial frame on the reference interval rather than an orthogonal basis, so the least-squares systems are ill-conditioned and require regularization. Coefficients are recovered via a GTSVD with exponential weight matrix $R_\beta=\operatorname{diag}(e^{\beta j})$, which selects smooth coefficient vectors by damping high-index modes. Because all subintervals share the same reference nodes and frame, the factorization is computed once and reused.

The Caputo derivative is approximated by applying the fractional integral to the derivative of the local reconstruction. Since derivatives of the frame functions are polynomials, each local weight reduces to finite moments

$$J_r^{(\alpha)}=\int_{t_L}^{t_U}(\lambda-t)^{-\alpha}t^r\,dt,$$

which admit closed-form expressions via the auxiliary integrals $H_s(\delta,L)$, including the singular case $I_k$0. The weak singularity is therefore handled analytically rather than by quadrature. For uniform partitions evaluated at identical relative positions within each block, the current-block weights depend only on $I_k$1 and the history weights only on the lag $I_k$2, converting the nonlocal Caputo history into one-sided discrete convolutions over the block index. Storage for the weights is $I_k$3 for $I_k$4 total samples, and FFT-based history accumulation costs $I_k$5, compared with $I_k$6 storage for a dense global differentiation matrix.

Error analysis

The analysis separates cleanly into two parts. First, the Caputo error admits the exact identity as a fractional integral of the derivative reconstruction error, yielding the stability bound

$I_k$7

so the singular kernel introduces no additional discretization error once the moment weights are exact; convergence of the Caputo approximation follows entirely from convergence of the derivative reconstruction. Second, the derivative error is bounded through an interpolation argument: given an $I_k$8 reconstruction error $I_k$9 and a weighted smoothness bound $[0,b]$0 on the error, the derivative error satisfies $[0,b]$1. Under an exponential coefficient-weight bound on the error and exponential decay of the $[0,b]$2 reconstruction error—appropriate for analytic local functions—the derivative error retains exponential-type decay, which the Caputo estimate then inherits. A notable feature is that the analysis is formulated for the actual regularized GTSVD reconstruction, not identified with an ordinary Legendre projection; the authors are explicit that this identification does not hold. Localization also improves approximability: under the affine map, an oscillatory component $[0,b]$3 acquires effective frequency $[0,b]$4, so local coefficient decay accelerates as blocks shrink. A consistency proposition confirms that if the local derivative is represented exactly, the moment weights reproduce the Caputo derivative exactly.

Numerical results

The experiments use fixed parameters throughout ($[0,b]$5, sampling ratio $[0,b]$6, $[0,b]$7, truncation threshold $[0,b]$8, $[0,b]$9), inherited from a prior parameter study rather than tuned per example.

Moment verification: For polynomials $\{P_j^{(T)}|_{[-1,1]}\}_{j=0}^{N_b-1}$0 with exactly representable frame coefficients, errors against closed-form Caputo derivatives remain at roundoff level (below $\{P_j^{(T)}|_{[-1,1]}\}_{j=0}^{N_b-1}$1 relative even at $\{P_j^{(T)}|_{[-1,1]}\}_{j=0}^{N_b-1}$2), confirming that observed errors in subsequent tests stem solely from the reconstruction stage.

Comparison with L1: With comparable sample counts ($\{P_j^{(T)}|_{[-1,1]}\}_{j=0}^{N_b-1}$3 vs. 241), the LLF-Caputo method achieves errors between $\{P_j^{(T)}|_{[-1,1]}\}_{j=0}^{N_b-1}$4 and $\{P_j^{(T)}|_{[-1,1]}\}_{j=0}^{N_b-1}$5 across $\{P_j^{(T)}|_{[-1,1]}\}_{j=0}^{N_b-1}$6, $\{P_j^{(T)}|_{[-1,1]}\}_{j=0}^{N_b-1}$7, and a composite analytic function, while L1 errors range from $\{P_j^{(T)}|_{[-1,1]}\}_{j=0}^{N_b-1}$8 to $\{P_j^{(T)}|_{[-1,1]}\}_{j=0}^{N_b-1}$9 — roughly eight to ten orders of magnitude improvement, reflecting algebraic versus roundoff-limited convergence.

Comparison with global Chebyshev spectral benchmark: The Chebyshev method on global nonuniform nodes reaches $P_j^{(T)}(t)=T^{-1/2}p_j(t/T)$0-level accuracy with fewer nodes for analytic functions, and the paper states plainly that global spectral methods remain more accurate when such nodes are available. However, the Chebyshev benchmark degrades sharply for $P_j^{(T)}(t)=T^{-1/2}p_j(t/T)$1 when nodes are reduced (error rising to $P_j^{(T)}(t)=T^{-1/2}p_j(t/T)$2 at 29 nodes), whereas the LLF method maintains $P_j^{(T)}(t)=T^{-1/2}p_j(t/T)$3-level accuracy from local equispaced data. This positions the method as targeting settings where equispaced data, localized structure, or blockwise marching are required.

Cost: Weight storage grows linearly with sample count while a dense-matrix proxy grows quadratically; direct summation timing favors the LLF implementation beyond roughly 450 samples ($P_j^{(T)}(t)=T^{-1/2}p_j(t/T)$4).

Fractional differential equation: Embedding the construction into a block solver for $P_j^{(T)}(t)=T^{-1/2}p_j(t/T)$5 with reduced rank $P_j^{(T)}(t)=T^{-1/2}p_j(t/T)$6 yields $P_j^{(T)}(t)=T^{-1/2}p_j(t/T)$7-level solution errors for manufactured solutions, several orders better than an L1 solver at equal unknown count, demonstrating feasibility of block time-marching.

Limitations and open questions

Several restrictions are acknowledged directly. The exponential convergence theory assumes analytic local functions and relies on the exponential weighted bound holding for the GTSVD reconstruction error; the paper does not prove this bound from first principles but takes it as the characteristic output of the regularization. The complexity advantage depends on fixed local parameters, and the reported timings are implementation-dependent. Noisy data, solutions with initial weak singularities, and nonsmooth functions are not treated: the authors note that handling piecewise smooth functions would require singularity detection, interface alignment, and one-sided reconstruction, and that truncation levels, weight parameters, singularity subtraction, and graded partitions for such cases await systematic study. Nonlinear fractional ODEs and fractional PDEs are likewise outside the present scope.

Conclusion

The paper provides a complete pipeline—from equispaced local data to Caputo derivatives—with analytic treatment of the singular kernel, a stability-oriented regularized reconstruction whose error propagates to the Caputo output without additional discretization loss, and a block-convolution implementation with linear storage and nearly linear application cost. The numerical evidence supports both the exactness of the moment weights and the practical accuracy claims, while honestly situating the method below global spectral benchmarks in the regime where nonuniform global nodes are available. The main open questions concern extension to nonsmooth and noisy data and to nonlinear and multidimensional fractional problems.