The Hilfer fractional derivative generalizes and interpolates between the commonly used Riemann–Liouville and Caputo fractional derivative. In general, solutions to Hilfer fractional derivative initial value problems are singular for t↓t₀, providing a challenge in numerical approximation. In this work, we present a Bernstein splines method to obtain accurate approximations to solutions of nonlinear Hilfer fractional derivative initial value problems, including solutions with singular behavior. We provide explicit convergence requirements and asymptotic convergence rates to analytical solutions. Moreover, the method is efficiently implemented, using a vectorized approach and parallelization in polynomial order. Numerical experiments show empirical convergence corresponding to analytical results and a significantly higher accuracy compared to the commonly used fractional Adams–Bashforth–Moulton predictor–corrector method. Furthermore, our method is able to simulate nonlinear systems with Hilfer fractional derivatives, as demonstrated by the application to the fractional Van der Pol oscillator.

We prove conditions for existence of analytical solutions for boundary value problems with the Hilfer fractional derivative, generalizing the commonly used Riemann-Liouville and Caputo operators. The boundary values, referred to in this paper as fractional-periodic, are fractional integral conditions generalizing recurrent solution values for the non-Caputo case of the Hilfer fractional derivative. Analytical solutions to the studied problem are obtained using a perturbation of the corresponding initial value problem with enforced boundary conditions. In general, solutions to the boundary value problem are singular for t ↓0. To overcome this singularity, we construct a sequence of converging solutions in a weighted continuous function space. We present a Bernstein spline-based implementation to numerically approximate solutions. We prove convergence of the numerical method, providing convergence criteria and asymptotic convergence rates. Numerical examples show empirical convergence results corresponding with the theoretical bounds. Moreover, the method is able to approximate the singular behavior of solutions and is demonstrated to converge for nonlinear problems. Finally, we apply a grid search to obtain correspondence to the original, non-perturbed system.