RL – Riemann-Liouville Fractional Integral/Derivative (All Points)

Overview

RL computes the Riemann-Liouville (RL) fractional integral or derivative of a function at all points in a specified domain, using an optimized vectorized approach with the trapezoidal rule. This makes it possible to evaluate the RL operator efficiently on large arrays, providing up to a 14x speedup over naïve implementations.

---

Mathematical Background

The Riemann-Liouville (RL) fractional differintegral of order $\alpha$ is defined as:

$$ I^\alpha f(x) = \frac{1}{\Gamma(\alpha)} \int_{a}^{x} (x-t)^{\alpha-1} f(t) , dt $$

Discretizing this integral on a uniform grid using the trapezoidal rule yields:

$$ RL^\alpha f(x_k) \approx h^{-\alpha} \sum_{j=0}^{k} D_{k,j} , f(x_j) $$

where $D_{k,j}$ are weights precomputed for each $(k, j)$ pair, depending on $\alpha$ and the grid size. This implementation builds the RL matrix efficiently, then applies it to the function values in a single matrix multiplication.

---

Implementation Details

• Input Flexibility: Accepts a callable, NumPy array, or list of function values. If a function is provided, it is evaluated at num_points equally spaced points between domain_start and domain_end.

• Fully Vectorized: All computations are vectorized, leveraging NumPy’s fast matrix operations.

• Trapezoidal Rule: The RL operator is discretized using the trapezoidal rule for improved accuracy compared to simple Riemann sums.

• Optimized RL Matrix: The RL weights are precomputed in a matrix (RLmatrix), and the operation is performed via fast matrix multiplication.

---

Function Signature

RL(
    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): Order of the fractional differintegral to compute.
• f_name (Callable, list, or np.ndarray): The function or sequence to be differintegrated.
• domain_start (float, optional): Start of the domain. Default is 0.0.
• domain_end (float, optional): End of the domain. Default is 1.0.
• num_points (int, optional): Number of discretization points. Default is 100.

---

Returns

• result (np.ndarray): Array containing the RL fractional integral/derivative at each grid point.

---

Example Usage

import numpy as np
from differintP import RL

# Fractional integral/derivative of sqrt(x) at all points in [0, 1]
RL_sqrt = RL(0.5, np.sqrt, 0., 1., 100)

# Fractional integral/derivative of a polynomial
RL_poly = RL(0.5, lambda x: x**2 - 1, 0., 1., 100)

---

Notes

• For the RL operator at a single endpoint, see RLpoint.
• The underlying RL matrix (RLmatrix) is fully vectorized for high performance.
• Suitable for large arrays and high-precision computations.

---

See Also

• RLpoint: RL operator at the endpoint.
• GL: Grünwald-Letnikov derivative (all points).
• GLpoint: GL endpoint derivative.

---

References

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