GLI – Improved Grünwald-Letnikov Fractional Derivative

Overview

GLI computes the improved Grünwald-Letnikov (GL) fractional derivative of a function over a specified domain, using quadratic (3-point Lagrange) interpolation for higher accuracy. It implements the method described in:

Oldham, K. & Spanier, J. (1974). The Fractional Calculus: Theory and Applications of Differentiation and Integration to Arbitrary Order. Academic Press.

This approach provides a higher-order, more accurate approximation of the fractional derivative than the standard GL method, particularly when working with smooth functions or limited discretization points.

---

Mathematical Background

The improved GL algorithm utilizes a 3-point Lagrange interpolation to better estimate the function's values between grid points, reducing discretization error. For each interior point, the interpolated value is given by:

$$ \text{interpolated}[i] = a \cdot f[i-1] + b \cdot f[i] + c \cdot f[i+1] $$

where the coefficients depend on the fractional order $\alpha$:

• $a = \frac{\alpha(\alpha - 2)}{8}$
• $b = \frac{4 - \alpha^2}{4}$
• $c = \frac{\alpha(\alpha + 2)}{8}$

These interpolated values are then convolved with the Grünwald-Letnikov binomial coefficients and scaled by the step size raised to $-\alpha$.

---

Implementation Details

• Hybrid Convolution: For small arrays (num_points < 800), a direct NumPy convolution (np.convolve) is used. For larger arrays, a fast FFT-based convolution (scipy.signal.fftconvolve) is automatically selected for optimal performance.
• Input Flexibility: GLI accepts a callable function, NumPy array, or Python list. If a callable is provided, it is evaluated at num_points equally spaced points between domain_start and domain_end.

---

Function Signature

GLI(
    alpha: float,
    f_name: Callable | np.ndarray | list,
    domain_start: float = 0.0,
    domain_end: float = 1.0,
    num_points: int = 100,
) -> np.ndarray

---

Parameters

• alpha (float): The order of the fractional derivative.
• f_name (Callable, list, or np.ndarray): Function, lambda, or sequence of function values to be differintegrated.
• domain_start (float, optional): Left endpoint of the domain. Default is 0.0.
• domain_end (float, optional): Right endpoint of the domain. Default is 1.0.
• num_points (int, optional): Number of points in the domain. Default is 100.

---

Returns

• result (np.ndarray): Array containing the improved GL fractional derivative at each grid point.

---

Example Usage

import numpy as np
from differintP import GLI

# Example 1: Fractional derivative of a quadratic polynomial
result_poly = GLI(-0.5, lambda x: x**2 - 1)

# Example 2: Fractional derivative of sqrt(x) on [0, 1]
result_sqrt = GLI(0.5, np.sqrt, 0., 1., 100)

---

Notes

• The interpolation and convolution steps are fully vectorized for speed.
• Hybrid convolution ensures performance is optimal for both small and large arrays.
• This method is especially useful when higher accuracy is desired with moderate grid resolution.

---

References

• Oldham, K. & Spanier, J. (1974). The Fractional Calculus: Theory and Applications of Differentiation and Integration to Arbitrary Order. Academic Press.