function_solve_VIDE

Solve the scalar Volterra integro-differential equation

\[y'(t) = a(t)\,y(t) + g(t) + \int_0^t K(t-s)\,y(s)\,ds,\quad y(0) = y_0\]

with callable inputs on an arbitrary mesh.

Parameters:

Name	Type	Description	Default
`kernel`	`callable`	`kernel(u)` returns \(K(u)\): a scalar, or a \((d, d)\) matrix for vector and matrix-valued equations.	required
`a`	`callable`	`a(t)` returns the coefficient \(a(t)\) (a scalar, or a \((d, d)\) matrix for vector/matrix equations). Defaults to zero. `a` does not depend on the right-hand side, so it is sampled once in the matrix case.	`None`
`g`	`callable`	`g(t)` returns the forcing term: scalar, \((d,)\), or \((d, m)\) for the matrix-valued case. Defaults to zero.	`None`
`soln_init_value`	`float or array_like`	\(y(0)\). Required. A scalar/`(d,)` value gives the scalar/vector case; a \((d, m)\) array selects the matrix-valued case (\(m\) right-hand sides).	required
`mesh_breakpoints`	`array_like`	Strictly-increasing 1-D array starting at 0.	required
`coll_divs`	`int, list of int`	Collocation node positions; see `function_solve_VIE_2` for the convention (differs from the array-based solvers). Defaults to `coll_divs=2, coll_choices=[0, 1, 2]`. Mutually exclusive with `coll_nodes`.	`None`
`coll_choices`	`int, list of int`	Collocation node positions; see `function_solve_VIE_2` for the convention (differs from the array-based solvers). Defaults to `coll_divs=2, coll_choices=[0, 1, 2]`. Mutually exclusive with `coll_nodes`.	`None`
`coll_nodes`	`array_like`	Collocation node positions given directly as floats in \([0, 1]\); see `function_solve_VIE_2`. Mutually exclusive with `coll_divs`/`coll_choices`.	`None`
`kernel_singularity`	`None, float, list of float, or callable`	Declare integrable singularities; see `function_solve_VIE_2`.	`None`
`return_function`	`bool`	If True, also return a callable solution wrapper.	`False`

Returns:

Name	Type	Description
`soln_values`	`ndarray`	\(y\) values at collocation nodes. Shape `(M, p)` (scalar), `(M, p, d)` (vector), or `(M, p, d, m)` (matrix-valued).
	`(soln_values, y_callable) : tuple`	When `return_function=True`.

Notes

Internally the unknowns are \(y'\) at the collocation nodes; the returned values are \(y\) itself (reconstructed via the antiderivative basis and the tracked boundary values).

Source code in src/voles/_callable_solvers.py

@_escalate_complex_warning
def function_solve_VIDE(*, kernel, a=None, g=None, soln_init_value,
                         mesh_breakpoints,
                         coll_divs: int = None, coll_choices: list[int] = None,
                         coll_nodes=None,
                         kernel_singularity=None,
                         return_function: bool = False,
                         show_warnings: bool = True,
                         _smooth_gl_order: int = 6):
    r"""Solve the scalar Volterra integro-differential equation

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

    with callable inputs on an arbitrary mesh.

    Parameters
    ----------
    kernel : callable
        ``kernel(u)`` returns $K(u)$: a scalar, or a $(d, d)$ matrix for vector
        and matrix-valued equations.
    a : callable, optional
        ``a(t)`` returns the coefficient $a(t)$ (a scalar, or a $(d, d)$ matrix
        for vector/matrix equations). Defaults to zero. ``a`` does not depend on
        the right-hand side, so it is sampled once in the matrix case.
    g : callable, optional
        ``g(t)`` returns the forcing term: scalar, $(d,)$, or $(d, m)$ for the
        matrix-valued case. Defaults to zero.
    soln_init_value : float or array_like
        $y(0)$. Required. A scalar/`(d,)` value gives the scalar/vector case; a
        $(d, m)$ array selects the matrix-valued case ($m$ right-hand sides).
    mesh_breakpoints : array_like
        Strictly-increasing 1-D array starting at 0.
    coll_divs, coll_choices : int, list of int, optional
        Collocation node positions; see ``function_solve_VIE_2`` for the
        convention (differs from the array-based solvers). Defaults to
        ``coll_divs=2, coll_choices=[0, 1, 2]``. Mutually exclusive with
        ``coll_nodes``.
    coll_nodes : array_like, optional
        Collocation node positions given directly as floats in $[0, 1]$; see
        ``function_solve_VIE_2``. Mutually exclusive with
        ``coll_divs``/``coll_choices``.
    kernel_singularity : None, float, list of float, or callable
        Declare integrable singularities; see ``function_solve_VIE_2``.
    return_function : bool, optional
        If True, also return a callable solution wrapper.

    Returns
    -------
    soln_values : ndarray
        $y$ values at collocation nodes. Shape ``(M, p)`` (scalar),
        ``(M, p, d)`` (vector), or ``(M, p, d, m)`` (matrix-valued).
    (soln_values, y_callable) : tuple
        When ``return_function=True``.

    Notes
    -----
    Internally the unknowns are $y'$ at the collocation nodes; the returned
    values are $y$ itself (reconstructed via the antiderivative basis and the
    tracked boundary values).
    """
    mesh_breakpoints = np.asarray(mesh_breakpoints, dtype=float)
    if mesh_breakpoints.ndim != 1 or len(mesh_breakpoints) < 2:
        raise ValueError("mesh_breakpoints must be 1-D with at least two entries")
    if not np.all(np.diff(mesh_breakpoints) > 0):
        raise ValueError("mesh_breakpoints must be strictly increasing")
    if mesh_breakpoints[0] != 0.0:
        raise ValueError("mesh_breakpoints[0] must be 0")

    node_pos, _, _ = _resolve_node_pos(
        coll_nodes, coll_divs, coll_choices,
        default_divs=2, default_choices=[0, 1, 2],
        exclude_zero=False, fname="function_solve_VIDE")

    if not callable(kernel):
        raise TypeError("kernel must be callable")
    if g is not None and not callable(g):
        raise TypeError("g must be callable or None")
    if a is not None and not callable(a):
        raise TypeError("a must be callable or None")

    # Complex dispatch.
    if _samples_indicate_complex([kernel, g, a], mesh_breakpoints, soln_init_value):
        sample_u = float(np.diff(mesh_breakpoints)[0]) * 0.5
        d_orig = _detect_complex_d_orig(kernel, sample_u)
        result = function_solve_VIDE(
            kernel=_block_wrap_kernel(kernel, d_orig),
            a=_block_wrap_a(a, d_orig),
            g=_block_wrap_g(g, d_orig),
            soln_init_value=_block_expand_init(soln_init_value, d_orig),
            mesh_breakpoints=mesh_breakpoints,
            coll_nodes=node_pos,
            kernel_singularity=kernel_singularity,
            return_function=return_function, show_warnings=show_warnings,
            _smooth_gl_order=_smooth_gl_order)
        if return_function:
            y_real, y_func_real = result
            return (_recombine_complex_y(y_real, d_orig),
                    _ComplexSolutionFunction(y_func_real, d_orig))
        return _recombine_complex_y(result, d_orig)

    M = len(mesh_breakpoints) - 1
    p = len(node_pos)
    widths = np.diff(mesh_breakpoints)

    from . import _dlang as _dlang_module
    max_p = _dlang_module.function_solve_max_p_d()
    if p > max_p:
        raise ValueError(
            f"number of collocation nodes p = {p} exceeds the maximum compiled "
            f"into the D extension ({max_p}).")

    # Detect scalar vs vector by sampling kernel at a non-singular point.
    sample_u = float(widths[0]) * 0.5
    is_vector, d = _detect_kernel_shape(kernel, sample_u)

    _maybe_warn_mesh_uniform_with_singularity(
        mesh_breakpoints, kernel_singularity, show_warnings)

    # Precomputed scalar tables: alpha[i,k] = I_k(c_i); w[k] = I_k(1)
    anti = _lagrange_antideriv_coefs(node_pos)  # (p, p+1)
    alpha = np.zeros((p, p), dtype=np.float64)
    for i in range(p):
        for k in range(p):
            alpha[i, k] = _eval_poly_at(anti[k], node_pos[i])
    w_vec = np.array([_eval_poly_at(anti[k], 1.0) for k in range(p)],
                     dtype=np.float64)

    # Sample g and a at collocation points
    def _sample_callable_scalar(f):
        out = np.zeros((M, p), dtype=np.float64)
        if f is None:
            return out
        for n in range(M):
            t_n = mesh_breakpoints[n]
            h_n = widths[n]
            for i in range(p):
                out[n, i] = float(f(t_n + node_pos[i] * h_n))
        return out

    if not is_vector:
        vide_basis = _vide_basis_coefs(node_pos)
        W = _build_W_with_basis_scalar(
            kernel, mesh_breakpoints, node_pos,
            kernel_singularity, _smooth_gl_order, vide_basis)
        g_arr = _sample_callable_scalar(g)
        a_arr = _sample_callable_scalar(a)
        y_prime, y_boundary = _dlang_module.function_solve_vide_d(
            W, g_arr, a_arr, alpha, w_vec, widths, float(soln_init_value))

        # Reconstruct y at collocation nodes: y_{n,i} = y_n + h_n * sum_k Y'_{n,k} alpha[i,k]
        y_at_coll = np.zeros((M, p), dtype=np.float64)
        for n in range(M):
            y_at_coll[n, :] = y_boundary[n] + widths[n] * (alpha @ y_prime[n])

        if return_function:
            polys = _build_vide_polynomials_scalar(
                y_prime, y_boundary, mesh_breakpoints, node_pos)
            y_func = _SolutionFunction(polys, mesh_breakpoints, d=0)
            return y_at_coll, y_func
        return y_at_coll

    # ----- Vector / matrix path -----
    max_d = _dlang_module.function_solve_max_d_d()
    if d > max_d:
        raise ValueError(
            f"kernel dimension d = {d} exceeds the maximum compiled into the "
            f"D extension ({max_d}).")

    # W and a(t) depend only on the kernel and a, not on the right-hand side,
    # so both are built once and shared across all columns in the matrix case.
    vide_basis = _vide_basis_coefs(node_pos)
    W = _build_W_with_basis_vector(
        kernel, mesh_breakpoints, node_pos,
        kernel_singularity, _smooth_gl_order, d, vide_basis)
    a_arr = _sample_a_at_coll(a, mesh_breakpoints, node_pos, widths, M, p, d)

    # Matrix-valued problem is signalled by a 2-D (d, m) initial value.
    init_arr = np.asarray(soln_init_value, dtype=np.float64)
    if init_arr.ndim == 2:
        if init_arr.shape[0] != d:
            raise ValueError(
                f"soln_init_value shape {tuple(init_arr.shape)} incompatible "
                f"with kernel dimension d = {d}.")
        m = init_arr.shape[1]
        G = _sample_g_at_coll_matrix(g, mesh_breakpoints, node_pos, widths,
                                     M, p, d, m)

        def _col_solve(j):
            g_arr_j = np.ascontiguousarray(G[:, :, :, j])
            return _dlang_module.function_solve_vide_vec_d(
                W, g_arr_j, a_arr, alpha, w_vec, widths,
                np.ascontiguousarray(init_arr[:, j]))

        with ThreadPoolExecutor(max_workers=m) as ex:
            results = list(ex.map(_col_solve, range(m)))

        # y_prime_j: (M, p, d); y_boundary_j: (M+1, d). Stack on a new last axis.
        y_prime = np.stack([r[0] for r in results], axis=3)      # (M, p, d, m)
        y_boundary = np.stack([r[1] for r in results], axis=2)   # (M+1, d, m)

        y_at_coll = np.zeros((M, p, d, m), dtype=np.float64)
        for n in range(M):
            # alpha (p,p) @ y_prime[n] (p,d,m) over the p axis -> (p,d,m)
            y_at_coll[n] = y_boundary[n] + widths[n] * np.tensordot(
                alpha, y_prime[n], axes=([1], [0]))

        if return_function:
            polys = _build_vide_polynomials_matrix(
                y_prime, y_boundary, mesh_breakpoints, node_pos,
                d, m)
            return y_at_coll, _SolutionFunction(polys, mesh_breakpoints,
                                                d=d, m=m)
        return y_at_coll

    # Vector path
    if init_arr.shape != (d,):
        raise ValueError(
            f"soln_init_value must have shape ({d},) for vector kernel; "
            f"got shape {tuple(init_arr.shape)}")

    g_arr = _sample_g_at_coll_vec(g, mesh_breakpoints, node_pos, widths, M, p, d)

    y_prime, y_boundary = _dlang_module.function_solve_vide_vec_d(
        W, g_arr, a_arr, alpha, w_vec, widths, init_arr)

    # Reconstruct y at collocation: y[n,i] = y_n + h_n * sum_k alpha[i,k] Y'_{n,k}
    y_at_coll = np.zeros((M, p, d), dtype=np.float64)
    for n in range(M):
        # alpha (p,p) @ y_prime[n] (p,d) -> (p,d)
        y_at_coll[n, :, :] = y_boundary[n] + widths[n] * (alpha @ y_prime[n])

    if return_function:
        polys = _build_vide_polynomials_vector(
            y_prime, y_boundary, mesh_breakpoints, node_pos, d)
        y_func = _SolutionFunction(polys, mesh_breakpoints, d=d)
        return y_at_coll, y_func
    return y_at_coll