A numerical Bernstein splines approach for nonlinear initial value problems with Hilfer fractional derivative

Technical Summary: A Numerical Bernstein Splines Approach for Nonlinear Hilfer Fractional Initial Value Problems

Problem Statement
The paper addresses the numerical approximation of nonlinear fractional initial value problems (FIVPs) involving the Hilfer fractional derivative. The Hilfer derivative, parameterized by order $\alpha \in (0,1)$ and type $\beta \in [0,1]$, interpolates between the Riemann-Liouville ($\beta=0$) and Caputo ($\beta=1$) derivatives. A primary challenge in solving these problems is that for $\beta \neq 1$ with non-homogeneous initial conditions, solutions exhibit singular behavior as $t \downarrow 0$ (specifically, they are non-continuous and behave like $t^{\gamma-1}$ where $\gamma = \alpha + \beta(1-\alpha)$). Existing numerical methods, such as the widely used Adams–Bashforth–Moulton predictor-corrector scheme, are often tailored for the Caputo case or require solutions to be continuously differentiable, limiting their applicability to the singular solutions arising in the general Hilfer case. Furthermore, the literature lacks efficient methods with theoretical error bounds capable of handling the full range of singular solutions for nonlinear Hilfer FIVPs.

Methodology
The authors propose a numerical method based on Bernstein splines to approximate solutions to the nonlinear FIVP:
$$\begin{cases} {}_{0}D^{\alpha,\beta}_{t} y(t) = f(t, y(t)), & 0 < t \leq T \\ {}_{0}I^{1-\gamma}_{t} y(0) = \tilde{y}_0 \end{cases}$$
The methodology involves several key components:

1. Weighted Space Transformation: To handle the singularity at $t=0$, the problem is transformed into a weighted continuous function space $C^{1-\gamma}[0, T]$. Additionally, a time-domain shift $\varepsilon > 0$ is introduced, solving the problem on $[\varepsilon, T]$ where the solution is continuous, and then taking the limit as $\varepsilon \downarrow 0$.
2. Bernstein Spline Approximation: The solution is approximated using piecewise Bernstein polynomials (splines) defined over a knot collection. The method employs a local "knot-for-knot" iteration strategy. On each sub-interval, the solution is constructed via Picard iterations applied to the equivalent integral equation derived from the Hilfer derivative definition.
3. Operational Matrix Implementation: The authors derive closed-form expressions for the fractional integration of Bernstein splines using the incomplete beta function. This allows the fractional integral operator to be applied analytically to the spline basis, avoiding numerical quadrature errors associated with the history of the solution.
4. Vectorized and Parallelized Computation: The implementation utilizes tensor products and vectorization (specifically using numpy.einsum in Python) to compute the fractional integrals efficiently. This approach allows for parallelization with respect to the polynomial order $q$, mitigating the computational cost typically associated with the "memory effect" of fractional derivatives.

Key Contributions

• Theoretical Framework: The paper establishes existence and uniqueness results for the spline approximations under relaxed convergence assumptions (Theorem 3.1). It provides explicit convergence rates for the approximation error with respect to the polynomial order $q$ ($O(q^{-\alpha/2})$), the knot size $h$ ($O(h^\alpha)$), and the time-shift parameter $\varepsilon$ ($O(\varepsilon^\alpha)$) (Theorem 3.2).
• Handling Singularities: Unlike previous works restricted to continuous solutions or the Caputo case, this method is explicitly designed to capture the singular behavior of solutions for $\beta \in [0, 1)$ by utilizing the weighted space transformation and the $\varepsilon$-shift.
• Efficient Implementation: The authors present a vectorized algorithm that decouples the computational complexity from the knot size in terms of polynomial basis evaluation, offering a promising avenue for speedup compared to traditional step-by-step methods.
• Nonlinear Application: The method is successfully applied to a nonlinear system, the fractional Van der Pol oscillator, demonstrating its capability to handle nonlinear terms and varying $\beta$ values.

Numerical Results

• Convergence Validation: Numerical experiments on a polynomial test problem confirm the theoretical convergence rates. The error decreases as the knot size $h$ and the shift $\varepsilon$ decrease, and as the polynomial order $q$ increases (up to a point where numerical stability limits further gains).
• Comparison with Predictor-Corrector: When compared to the Adams–Bashforth–Moulton predictor-corrector method (Diethelm et al.) on a linear Caputo case ($\beta=1$), the Bernstein spline method demonstrates significantly higher accuracy (lower mean and supremum errors) for the same knot sizes. The predictor-corrector method showed a persistent bias and did not improve its supremum error with finer grids in the tested setup.
• Van der Pol Oscillator: Simulations of the Hilfer-fractional Van der Pol oscillator showed that the method converges to stable limit cycles regardless of the Hilfer type $\beta$. The results indicated that the system exhibits strong attraction to a unique limit cycle, consistent with semi-analytical findings for the Caputo case.

Significance and Claims
The paper claims to provide the first efficient numerical method with theoretical error bounds capable of obtaining the full range of solutions (including singular ones) for nonlinear Hilfer fractional initial value problems. The authors emphasize that their approach removes the stringent constraints on the time domain interval size found in other methods by using local iterations. While acknowledging that the spline method is currently more computationally demanding than the predictor-corrector method, the authors highlight its superior accuracy and its ability to handle singular, nonlinear systems where standard methods fail or lack theoretical guarantees. The work lays the groundwork for future extensions to boundary value problems and fractional partial differential equations (FPDEs).