PCsolver – Fractional Predictor-Corrector Solver

Overview

PCsolver numerically solves fractional initial value problems (IVPs) of the form $ D^\alpha y(x) = f(x, y(x)) $ using a predictor-corrector approach, adapted to fractional orders. The method generalizes classical Adams-Bashforth-Moulton predictor-corrector techniques to support non-integer (fractional) derivatives. The implementation supports multi-term initial conditions and is suitable for orders $\alpha > 1$.

---

Mathematical Background

This solver implements the fractional Adams-Bashforth-Moulton predictor-corrector scheme, a widely-used method for numerically integrating fractional differential equations (FDEs). For an FDE of order $\alpha$:

$$ D^\alpha y(x) = f(x, y(x)), \quad y^{(k)}(0) = y_k, \quad k = 0, 1, \ldots, \lceil \alpha \rceil - 1 $$

the algorithm computes a discrete approximation on a uniform grid. At each step, the next value is predicted using the previous solution and then corrected by evaluating both the right-hand side at the new point and a weighted sum (with PCcoeffs) over previous values.

The PCcoeffs function generates the necessary weights for the convolution sum, generalizing the classic coefficients to fractional orders.

For further theoretical background, see:

• Deng, W. (2007). "Short memory principle and a predictor-corrector approach for fractional differential equations." J. Comput. Appl. Math.
• Baskonus, H.M., Bulut, H. (2015). "On the numerical solutions of some fractional ordinary differential equations..." De Gruyter.

---

Implementation Details

• Handles Arbitrary Fractional Orders ($1 < \alpha$): Supports IVPs for any order $\alpha > 1$, automatically handling the required number of initial values.
• Supports Multi-Term Initial Conditions: Accepts a list of initial values matching the integer part of the order.
• Flexible Right-Hand Side: The function $f(x, y)$ can be supplied as a lambda or function handle.
• Discrete Grid: Solution is computed at num_points equally spaced points between domain_start and domain_end.
• Predictor and Corrector: Each step first predicts and then corrects the solution, using fractional weights computed by PCcoeffs.
• Dependency: Relies on scipy.special.factorial and Gamma (from scipy.special).

---

Function Signature

PCsolver(
    initial_values: list[float],
    alpha: float,
    f_name: Callable,
    domain_start: float = 0.0,
    domain_end: float = 1.0,
    num_points: int = 100,
) -> np.ndarray

---

Parameters

• initial_values (list[float]): List of initial values for the IVP. The number should be at least ceil(alpha).

• alpha (float): The order of the fractional derivative/integral ($ \alpha > 1 $).

• f_name (Callable): Function or lambda specifying the right side of the differential equation. Must accept two arguments: x and y.

• 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 grid points. Default is 100.

---

Returns

• y_correction (np.ndarray): Array containing the solution to the fractional IVP at each point in the grid.

---

Example Usage

import numpy as np
from differintP.solvers import PCsolver

# Example: Solve D^1.5 y(x) = y - x - 1 with y(0)=1, y'(0)=1
f = lambda x, y: y - x - 1
initial_values = [1, 1]
y_approx = PCsolver(initial_values, 1.5, f)

# Compare to analytical solution for validation (if known)
x = np.linspace(0, 1, 100)
theoretical = x + 1
assert np.allclose(y_approx, theoretical)

---

Notes

• Only supports $\alpha > 1$ (fractional orders less than 1 are not yet implemented).
• The coefficients are computed via the PCcoeffs helper function, which generates the fractional Adams-Bashforth-Moulton weights.
• For more advanced or higher-order FDE solvers, see related pages or references.