RLpoint - iparsw/differintP GitHub Wiki

`RLpoint` – Riemann-Liouville Fractional Integral/Derivative at a Point

Overview

RLpoint efficiently computes the Riemann-Liouville (RL) fractional integral or derivative at the endpoint of a domain, using an optimized vectorized approach with the trapezoidal rule. This is ideal when you only require the RL operator at a single point (typically the rightmost grid point).

This method provides an 8x–60x speedup compared to naïve implementations, thanks to its fully vectorized coefficient calculations and use of NumPy’s optimized routines.

Mathematical Background

The Riemann-Liouville (RL) fractional differintegral of order $\alpha$ for a function $f(x)$ is classically defined as:

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

Discretizing this with the trapezoidal rule leads to a sum of the form:

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

where $c_{j,k}$ are precomputed weights that depend on $\alpha$, $k$, and $j$.

This implementation calculates these weights efficiently for the endpoint and applies them in a single dot product.

Implementation Details

Input Flexibility: Accepts a callable, NumPy array, or list of function values. If a function is given, it is evaluated at num_points equally spaced points between domain_start and domain_end.
Fully Vectorized: All coefficient and sum calculations are vectorized for maximum efficiency—no Python loops are used in the core computation.
Trapezoidal Rule: Uses the trapezoidal rule to improve the accuracy of the RL integral approximation.

Function Signature

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

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 (evaluation point). Default is 1.0.
num_points (int, optional): Number of discretization points. Default is 100.

Returns

result (float): The RL fractional integral/derivative at the endpoint.

Example Usage

import numpy as np
from differintP import RLpoint

# Fractional integral/derivative of sqrt(x) at x = 1
rl_value = RLpoint(0.5, np.sqrt, 0., 1., 100)

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

Notes

For the RL operator at all points, see the RL function.
This function is especially efficient for large num_points and high-precision needs.
Coefficient calculation is based on the discretized RL formula using the trapezoidal rule.