Skip to content

Example: Polynomial Solutions (return_function)

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 an explicit polynomial fixed by requiring the Volterra equation to hold at a set of collocation points within that interval.

Because the solution on each interval is an explicit polynomial, the solver can return it. Passing return_function=True to any solver returns a (soln_values, solution) tuple (return_polys=True is a deprecated alias). solution(t) evaluates the piecewise polynomial at any time, and solution also indexes/iterates like 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, solution = solve_VIDE(
    kernel_values=kernel,
    a_values=a,
    g_values=g,
    soln_init_value=0.0,
    time_step=time_step,
    return_function=True,
)

# solution(t) evaluates the piecewise polynomial directly:
print(f"y(0.2)  ≈ {solution(0.2):.6f},  exact = {np.sin(0.2):.6f}")

# solution also indexes like the per-interval polynomials:
p = solution[0]                         # numpy.polynomial.Polynomial on t ∈ [0, 0.4]

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