volterra-equation-solvers 0.3.2

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

volterra-equation-solvers

Collocation-method solvers for Volterra integral and integro-differential equations, based on:

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.

Solvers

solve_VIE_1

Given functions $K$ and $g$, solves for $y(t)$ in the Type-1 Volterra integral equation (VIE-1):

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

The solver handles the following cases:

• All functions are scalar-valued.
• $y(t)$ and $g(t)$ are $d$-dimensional vectors and $K(t)$ is a $d \times d$ matrix.
• $y(t)$ and $g(t)$ are $d \times m$ matrices and $K(t)$ is a $d \times d$ matrix.

For a full description of the solver inputs/outputs and available settings, see: API reference

solve_VIE_2

Given functions $K$ and $g$, solves for $y(t)$ in the Type-2 Volterra integral equation (VIE-2):

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

The solver handles the following cases:

• All functions are scalar-valued.
• $y(t)$ and $g(t)$ are $d$-dimensional vectors and $K(t)$ is a $d \times d$ matrix.
• $y(t)$ and $g(t)$ are $d \times m$ matrices and $K(t)$ is a $d \times d$ matrix.

For a full description of the solver inputs/outputs and available settings, see: API reference

solve_VIDE

Given functions $K$, $a$, and $g$ and an initial value $y(0)$, solves for $y(t)$ in the Volterra integro-differential equation (VIDE):

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

The solver handles the following cases:

• All functions are scalar-valued.
• $y(t)$ and $g(t)$ are $d$-dimensional vectors and $K(t)$ and $a(t)$ are $d \times d$ matrices.
• $y(t)$ and $g(t)$ are $d \times m$ matrices and $K(t)$ and $a(t)$ are $d \times d$ matrices.

For a full description of the solver inputs/outputs and available settings, see: API reference

Installation

pip install volterra-equation-solvers

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: numba, scipy (only needed for non-standard collocation settings not supported by the D extension)

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

Quick start

import numpy as np
from volterra_equation_solvers 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 three 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 solver detects which case applies automatically from the shapes of the input arrays.

import numpy as np
from volterra_equation_solvers 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 volterra_equation_solvers 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 volterra_equation_solvers 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.