voles 0.5.0

Collocation-method solvers for Volterra integral and integro-differential equations

voles

The Volterra Equation Solvers (VOLES) package is a collection of collocation-method solvers for Volterra integral and integro-differential equations. The algorithms used come from the book

Brunner H. Collocation Methods for Volterra Integral and Related Functional Differential Equations. Cambridge University Press; 2004.

The solvers are implemented as a compiled extension written in the D language. Performance should be on par with optimized C or FORTRAN code. All solvers support real-valued and complex-valued data, and scalar-, vector-, and matrix-valued equations.

Solvers

Two solver families are provided.

• The array-based family (solve_VIE_1, solve_VIE_2, solve_VIDE) is the fastest path when the kernel and forcing are already sampled on a uniform time grid. Supports scalar, vector, and matrix-valued solutions, and complex-valued data.

• The callable-input family (function_solve_VIE_1, function_solve_VIE_2, function_solve_VIDE) accepts the kernel, forcing term, and (for VIDE) coefficient $a$ as Python callables, runs on an arbitrary mesh via mesh_breakpoints, and handles weakly singular convolution kernels — useful for Abel-type problems where a graded mesh recovers the optimal convergence order. Supports scalar, vector, and matrix-valued cases, and complex-valued data. Requires scipy (install via the [callable] extra). The helper optimal_graded_mesh builds a Brunner-graded mesh for an $\alpha$-singular kernel.

Type-1 Volterra integral equation (VIE-1)

Given $K$ and $g$, solve for $y(t)$ in:

$$g(t) = \int_0^t K(t-s)\, y(s)\, ds$$

API	Inputs	Solution shape	Reference
solve_VIE_1	sampled arrays, uniform grid	scalar / vector / matrix	api/vie1/
function_solve_VIE_1	callables, arbitrary mesh	scalar / vector / matrix	api/function_vie1/

Type-2 Volterra integral equation (VIE-2)

Given $K$ and $g$, solve for $y(t)$ in:

$$y(t) = g(t) + \int_0^t K(t-s)\, y(s)\, ds$$

API	Inputs	Solution shape	Reference
solve_VIE_2	sampled arrays, uniform grid	scalar / vector / matrix	api/vie2/
function_solve_VIE_2	callables, arbitrary mesh	scalar / vector / matrix	api/function_vie2/

Volterra integro-differential equation (VIDE)

Given $K$, $a$, $g$, and initial value $y(0)$, solve for $y(t)$ in:

$$y'(t) = a(t)\, y(t) + g(t) + \int_0^t K(t-s)\, y(s)\, ds$$

API	Inputs	Solution shape	Reference
solve_VIDE	sampled arrays, uniform grid	scalar / vector / matrix	api/vide/
function_solve_VIDE	callables, arbitrary mesh	scalar / vector / matrix	api/function_vide/

Mesh helper: optimal_graded_mesh

Returns a Brunner-graded mesh $t_n = T (n/M)^r$ with grading exponent $r = p / (1 - \alpha)$, where $p$ is the number of collocation nodes per interval. Designed for weakly singular convolution kernels $K(u) \sim u^{-\alpha}$ with $\alpha \in [0, 1)$. Feed the result to a callable-input solver via mesh_breakpoints to recover the optimal convergence order on Abel-type problems.

API reference: api/optimal_graded_mesh/

Installation

pip install voles

Pre-built wheels are provided for Linux x86_64, macOS arm64 (Apple Silicon), and Windows x64. The D extension is bundled in the wheel and requires no extra tooling. Intel Macs are no longer supported as of 0.3.2; users can pin to volterra-equation-solvers==0.3.1 or build from source (see CONTRIBUTING.md).

Requirements: Python ≥ 3.10, numpy Optional extras:

• [callable] (scipy) — required for the function_solve_* family.
• [numba] (numba, scipy) — only needed for the array-based solvers when using non-standard collocation settings not compiled into the D extension.

To build from source (e.g. on an unsupported platform), see CONTRIBUTING.md.

Quick start

import numpy as np
from voles import solve_VIE_2

# y(t) = sin(t) satisfies this VIE-2 with K(s) = exp(-s)
time_step = 0.05
times = np.arange(0, 2.1, time_step)   # 42 points
kernel = np.exp(-times)
g = np.sin(times) - 0.5*(np.exp(-times) + np.sin(times) - np.cos(times))

# Default solver settings require length of form 4k+1; input will be
# truncated from 42 to 41, so soln has 41 elements, not 42.
soln = solve_VIE_2(
    kernel_values=kernel,
    g_values=g,
    time_step=time_step,
)
print(f"Max error: {max(abs(soln - np.sin(times[:len(soln)]))):.2e}")

All solvers accept return_polys=True to also return the piecewise polynomial solution as a list of numpy.polynomial.Polynomial objects.

The solvers require input arrays to satisfy an internal size constraint. Any length can be passed; if the length doesn't meet the constraint, the arrays are automatically truncated to the nearest valid length and a warning is printed. See the API reference for each solver for details.

Vector and Matrix Valued Equations

All solvers can solve for vector-valued and matrix-valued functions $y(t)$. When $y(t)$ is a $d$-dimensional vector, $g(t)$ is also a $d$-dimensional vector and $K(t)$ and $a(t)$ are $d \times d$ matrices. When $y(t)$ is a $d \times m$ matrix, $g(t)$ is also a $d \times m$ matrix and $K(t)$ and $a(t)$ are $d \times d$ matrices. The case is detected automatically: for the array-based family from the shapes of the input arrays, and for the callable family from the shape returned by g(t) (a (d, m) return — or a (d, m) soln_init_value for VIDE — selects the matrix case). The callable family builds the kernel weight tensor once and shares it across the $m$ columns, so a matrix solve is much cheaper than $m$ separate calls; see the callable-solver examples for a worked case.

import numpy as np
from voles import solve_VIE_1

# 2×2 VIE-1 with constant kernel K = [[3/2, -1/2], [-1/2, 3/2]],
# g(t) = [t + (3/2)t², t - (1/2)t²], and exact solution y(t) = [1+2t, 1]
time_step = 0.1
times = np.arange(0, 9.1, time_step)   # 91 pts = 10×3² + 1
N = len(times)

kernel = np.full((N, 2, 2), [[1.5, -0.5], [-0.5, 1.5]])

g = np.zeros((N, 2))
g[:, 0] = times + 1.5 * times**2
g[:, 1] = times - 0.5 * times**2

soln = solve_VIE_1(kernel_values=kernel, g_values=g, time_step=time_step)
# soln shape: (N, 2)
exact = np.column_stack([1 + 2*times, np.ones(N)])
print(f"Max error: {np.max(np.abs(soln - exact)):.2e}")

Complex-Valued Equations

All three solvers accept complex-valued inputs. Pass complex NumPy arrays for the kernel, forcing function, and (for VIDE) initial value, and the solver returns a complex-valued solution. This works for scalar, vector, and matrix cases alike.

import numpy as np
from voles import solve_VIE_2

time_step = 0.05
times = np.arange(0, 2.1, time_step)
kernel = np.exp(-1j * times)               # complex kernel
g = np.ones_like(times, dtype=complex)

soln = solve_VIE_2(kernel_values=kernel, g_values=g, time_step=time_step)
# soln is a complex-valued array

How the Collocation Method Works

The solvers approximate each component of $y(t)$ as a piecewise polynomial. The time axis is divided into mesh intervals, and on each interval the solution is represented by a polynomial whose coefficients are determined by requiring the Volterra equation to hold exactly at a set of collocation points within that interval.

Because the solution on each mesh interval is an explicit polynomial, the solver can optionally return it (see Polynomial Solutions below). This is useful for evaluating the solution at arbitrary times, differentiating, integrating, and so on.

Piecewise Polynomial Illustration

The figure below illustrates the idea. The time axis is split into mesh intervals (dashed vertical lines), and on each interval the solver finds a polynomial $p_i(t)$ that satisfies the equation at a set of collocation points (dots). Stitching the pieces together produces a piecewise polynomial that closely tracks the exact solution.

Polynomial Solutions

Passing return_polys=True to any solver returns a (soln_values, polys) tuple, where polys is a list of numpy.polynomial.Polynomial objects covering successive mesh intervals. These can be evaluated at any point, differentiated, integrated, and so on. The following example uses solve_VIDE to solve for $y(t) = \sin(t)$, then evaluates the solution and its derivative at a point not on the time grid:

import numpy as np
from voles import solve_VIDE

# y(t) = sin(t) satisfies this VIDE with K(s) = exp(-s), a(t) = -1
time_step = 0.1
times = np.arange(0, 9.1, time_step)   # 91 points
kernel = np.exp(-times)
a = np.full(len(times), -1.0)
g = 1.5*np.cos(times) + 0.5*np.sin(times) - 0.5*np.exp(-times)

soln_vals, polys = solve_VIDE(
    kernel_values=kernel,
    a_values=a,
    g_values=g,
    soln_init_value=0.0,
    time_step=time_step,
    return_polys=True,
)

# polys[i] is a numpy.polynomial.Polynomial covering mesh interval i.
# Evaluate at a point not on the time grid:
p = polys[0]                            # covers t ∈ [0, 0.4]
print(f"y(0.2)  ≈ {p(0.2):.6f},  exact = {np.sin(0.2):.6f}")

# Differentiate to recover y'(t):
print(f"y'(0.2) ≈ {p.deriv()(0.2):.6f},  exact = {np.cos(0.2):.6f}")

Benchmarks

All three solvers have the same expected asymptotic complexity in N, d, and m, where N is the number of input points, d is the number of rows in the solution, and m is the number of columns:

	Scalar	Vector (d×1)	Matrix (d×m)*
Time	O(N²)	O(N²d²)	O(N²d²m)
Memory	O(N)	O(Nd²)	O(Nd²m)

* The m columns of the solution are independent and the code runs them in parallel.

The quadratic time scaling arises because each new mesh step requires a history sum over all previous steps. The coll_divs and coll_choices parameters affect the constant factor but not the asymptotic scaling in N, d, and m.

Run on a GitHub Actions ubuntu-22.04 runner (2-core x86_64 VM on an Intel Xeon 8370C, 2.8 GHz base / 3.5 GHz boost). Mean time is averaged over a variable number of calibrated rounds (from ~9 for large inputs up to ~6000 for small inputs).

See the Getting Started page for complete examples.

Worked derivations of the analytic solutions used in the test suite are in docs/scalar_solutions.pdf and docs/coupled_vector_solutions.pdf.