Fractional time derivatives¶

Caputo relaxation¶

\[ D_C^{0.6}u+u=0,\qquad u(0)=1, \]

with 48 Birk-Song modes and the default eliminated recurrence:

F = (inner(CaputoDerivative(u, 0.6), v) + inner(u, v)) * dx
stepper = FractionalTimeStepper(F, BirkSong(48), t, dt, u)

Increment t after each successful advance(). Representation error is controlled by the mode count, temporal error by dt, and algebraic error by the Firedrake solver parameters. See Separating error sources.

Changing the representation¶

Only one argument changes:

stepper = FractionalTimeStepper(F, Diethelm(48), t, dt, u)

Equal mode counts do not imply equal error. To compare fairly, fix the timestep and problem, refine each mode count independently, and measure both against the same analytic or high-accuracy reference. python -m benchmarks.diffusive_representations does this across every available spectrum and preserves the numerical results as CSV.

Use direct history instead, with no quadrature spectrum to tune:

stepper = FractionalTimeStepper(F, FullHistory(), t, dt, u)

This stores one solution increment per accepted step, so cost grows with the number of steps. It is the natural reference method (How the memory is advanced).

Changing the formulation¶

AuxiliaryODE solves the physical field and all modes together in \(m+1\) mixed fields, which makes the memory available to PETSc:

stepper = FractionalTimeStepper(
    F,
    BirkSong(8),
    t,
    dt,
    u,
    formulation=AuxiliaryODE(scheme="backward_euler"),
)

"trapezoidal" is the other supported scheme. The default recurrence instead solves one physical field and stores \(m\) histories, and is usually cheaper. caputo_auxiliary_ode.py shows the field split. See Solvers, field splits, and MPI for the PETSc side.

Riemann-Liouville¶

Switch only the marker:

F = (
    inner(RiemannLiouvilleDerivative(u, alpha), v)
    + inner(u, v)
) * dx
stepper = FractionalTimeStepper(F, BirkSong(48), t, dt, u)

The stepper evaluates

\[ D_{RL}^{\alpha}u(t) =D_C^\alpha u(t) +\frac{u(t_0)}{\Gamma(1-\alpha)}(t-t_0)^{-\alpha}, \]

where \(t_0\) is the value of t at construction and the trace term is exact. stepper.reset(u0, t0=new_time) establishes a new lower limit and trace. Both formulations and all representations support it, and FullHistory() uses the default eliminated formulation.

Remember that the Riemann-Liouville derivative of a nonzero constant is not zero. The interface assumes a classical trace \(u(t_0)\). Arbitrary prehistory and fractional-integral initial data are not represented (Caputo and Riemann-Liouville derivatives).

Exponential memory¶

ExponentialMemory and CaputoFabrizioOperator use TimeMemoryStepper instead, and are not fractional derivatives. See Exponential (fading) memory.