function_solve_VIE_1

Solve the Volterra integral equation of the first kind

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

with callable kernel and right-hand side on an arbitrary mesh. Zero is not a permitted collocation node (both sides of the equation vanish at t=0).

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
g	callable	g(t) returns the right-hand side: scalar, \((d,)\), or \((d, m)\) for the matrix-valued case (\(m\) right-hand sides). Defaults to zero (trivial y=0). A 2-D return signals the matrix case.	None
soln_init_value	float or array_like	\(y(0)\). Required only when force_continuous=True; ignored otherwise. A warning is emitted if a value is passed when it has no effect. For a matrix-valued problem with force_continuous=True it must have shape \((d, m)\).	None
mesh_breakpoints	array_like	Strictly-increasing 1-D array starting at 0.	required
coll_divs	int, list of int	Collocation nodes lie at coll_choices[k] / coll_divs in (0, 1]. Zero is excluded from coll_choices. Defaults to coll_divs=3, coll_choices=[1, 2, 3]. Mutually exclusive with coll_nodes.	None
coll_choices	int, list of int	Collocation nodes lie at coll_choices[k] / coll_divs in (0, 1]. Zero is excluded from coll_choices. Defaults to coll_divs=3, coll_choices=[1, 2, 3]. Mutually exclusive with coll_nodes.	None
coll_nodes	array_like	Collocation node positions given directly as floats in \((0, 1]\) (zero excluded); see function_solve_VIE_2. Mutually exclusive with coll_divs/coll_choices. The convergence check (see Notes) is applied to both the integer and the coll_nodes paths.	None
force_continuous	bool	If True, use the continuous collocation method (Brunner's \(S_m^{(0)}\)): the solution is a globally \(C^0\) piecewise polynomial of degree $m = $ len(coll_choices), anchored at $y(0) = $ soln_init_value. This requires the last node to be the right endpoint, \(c_m = 1\) (see Notes for the convergence condition). The default discontinuous method (\(S_{m-1}^{(-1)}\)) imposes no such restriction and is generally the better default.	False
kernel_singularity		See function_solve_VIE_2.	None
return_function		See function_solve_VIE_2.	None
show_warnings		See function_solve_VIE_2.	None

Notes

Whether a given node set yields a convergent method depends on the nodes (Brunner 2004, smooth-kernel chapter). For a smooth kernel with \(|K(t, t)| \ge k_0 > 0\), writing \(c_1 < \dots < c_m\) for the nodes:

• Discontinuous (default): converges iff \(|\rho_m| := \prod_{i=1}^{m} (1 - c_i)/c_i \le 1\) (Thm 2.4.2), with global order \(m\) (reduced to \(m - 1\) at \(\rho_m = 1\)).
• Continuous (force_continuous): requires \(c_m = 1\) and converges iff \(|\rho_{m-1}| := \prod_{i=1}^{m-1} (1 - c_i)/c_i \le 1\) (Thm 2.4.5), with global order \(m + 1\) (reduced to \(m\) at \(\rho_{m-1} = 1\)).

Node sets violating these are rejected with a ValueError. The criteria do not apply to weakly-singular kernels (where \(|K(t, t)|\) is unbounded), so the check is skipped when kernel_singularity is given and convergence is then the caller's responsibility.

These theorems are proved for a uniform mesh. The amplification factors \(\rho\) depend only on the nodes, so the conditions carry over to quasi-uniform meshes (all interval widths within fixed constant multiples of each other). They are not guaranteed on strongly graded meshes -- e.g. the geometric refinement of optimal_graded_mesh, where the smallest interval shrinks far faster than the largest -- but such meshes normally accompany a declared kernel_singularity, for which the check is already skipped.

Source code in src/voles/_callable_solvers.py

@_escalate_complex_warning
def function_solve_VIE_1(*, kernel, g=None, soln_init_value=None,
                          mesh_breakpoints,
                          coll_divs: int = None, coll_choices: list[int] = None,
                          coll_nodes=None,
                          kernel_singularity=None,
                          return_function: bool = False,
                          force_continuous: bool = False,
                          show_warnings: bool = True,
                          _smooth_gl_order: int = 6):
    r"""Solve the Volterra integral equation of the first kind

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

    with callable kernel and right-hand side on an arbitrary mesh. Zero is not
    a permitted collocation node (both sides of the equation vanish at t=0).

    Parameters
    ----------
    kernel : callable
        ``kernel(u)`` returns $K(u)$: a scalar, or a $(d, d)$ matrix for vector
        and matrix-valued equations.
    g : callable, optional
        ``g(t)`` returns the right-hand side: scalar, $(d,)$, or $(d, m)$ for
        the matrix-valued case ($m$ right-hand sides). Defaults to zero
        (trivial y=0). A 2-D return signals the matrix case.
    soln_init_value : float or array_like, optional
        $y(0)$. Required only when ``force_continuous=True``; ignored
        otherwise. A warning is emitted if a value is passed when it has no
        effect. For a matrix-valued problem with ``force_continuous=True`` it
        must have shape $(d, m)$.
    mesh_breakpoints : array_like
        Strictly-increasing 1-D array starting at 0.
    coll_divs, coll_choices : int, list of int, optional
        Collocation nodes lie at ``coll_choices[k] / coll_divs`` in (0, 1].
        Zero is excluded from ``coll_choices``. Defaults to
        ``coll_divs=3, coll_choices=[1, 2, 3]``. Mutually exclusive with
        ``coll_nodes``.
    coll_nodes : array_like, optional
        Collocation node positions given directly as floats in $(0, 1]$ (zero
        excluded); see ``function_solve_VIE_2``. Mutually exclusive with
        ``coll_divs``/``coll_choices``. The convergence check (see Notes) is
        applied to both the integer and the ``coll_nodes`` paths.
    force_continuous : bool, optional
        If ``True``, use the continuous collocation method (Brunner's
        $S_m^{(0)}$): the solution is a globally $C^0$ piecewise polynomial of
        degree $m = $ ``len(coll_choices)``, anchored at $y(0) = $
        ``soln_init_value``. This requires the last node to be the right
        endpoint, $c_m = 1$ (see Notes for the convergence condition). The
        default discontinuous method ($S_{m-1}^{(-1)}$) imposes no such
        restriction and is generally the better default.
    kernel_singularity, return_function, show_warnings :
        See ``function_solve_VIE_2``.

    Notes
    -----
    Whether a given node set yields a convergent method depends on the nodes
    (Brunner 2004, smooth-kernel chapter). For a smooth kernel with
    $|K(t, t)| \ge k_0 > 0$, writing $c_1 < \dots < c_m$ for the nodes:

    - Discontinuous (default): converges iff
      $|\rho_m| := \prod_{i=1}^{m} (1 - c_i)/c_i \le 1$ (Thm 2.4.2), with global
      order $m$ (reduced to $m - 1$ at $\rho_m = 1$).
    - Continuous (``force_continuous``): requires $c_m = 1$ and converges iff
      $|\rho_{m-1}| := \prod_{i=1}^{m-1} (1 - c_i)/c_i \le 1$ (Thm 2.4.5), with
      global order $m + 1$ (reduced to $m$ at $\rho_{m-1} = 1$).

    Node sets violating these are rejected with a ``ValueError``. The criteria
    do not apply to weakly-singular kernels (where $|K(t, t)|$ is unbounded), so
    the check is skipped when ``kernel_singularity`` is given and convergence is
    then the caller's responsibility.

    These theorems are proved for a *uniform* mesh. The amplification factors
    $\rho$ depend only on the nodes, so the conditions carry over to
    *quasi-uniform* meshes (all interval widths within fixed constant multiples
    of each other). They are not guaranteed on *strongly graded* meshes -- e.g.
    the geometric refinement of ``optimal_graded_mesh``, where the smallest
    interval shrinks far faster than the largest -- but such meshes normally
    accompany a declared ``kernel_singularity``, for which the check is already
    skipped.
    """
    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, _divs, _choices = _resolve_node_pos(
        coll_nodes, coll_divs, coll_choices,
        default_divs=3, default_choices=[1, 2, 3],
        exclude_zero=True, fname="function_solve_VIE_1")

    # Reject node sets for which the chosen VIE-1 method does not converge.
    # The criteria assume a smooth kernel (Brunner Thm 2.4.2(b): |K(t,t)|>=k0>0),
    # so they are skipped when a weakly-singular kernel is declared.
    if kernel_singularity is None:
        if force_continuous:
            # Continuous S_m^(0): converges iff c_m = 1 and |rho_{m-1}| <= 1.
            if abs(node_pos[-1] - 1.0) > 1e-12:
                raise ValueError(
                    "force_continuous (continuous S_m^(0) collocation) requires the "
                    "last collocation node to be the right endpoint c_m = 1; got "
                    f"c_m = {node_pos[-1]:.6g}. Append 1.0 to your nodes.")
            rho = _vie1_cont_amplification(node_pos)
            if rho > 1.0 + 1e-9:
                raise ValueError(
                    "Collocation nodes do not yield a convergent continuous (S_m^(0)) "
                    f"VIE-1 solver: the leading-node amplification |rho_{{m-1}}| = "
                    f"prod_(i<m) (1 - c_i)/c_i = {rho:.4g} exceeds 1 (Brunner Thm 2.4.5).")
        elif not _vie1_convergent(node_pos):
            # Discontinuous S_{m-1}^{(-1)}: converges iff |rho_m| <= 1.
            rho = _vie1_amplification(node_pos)
            if _choices is not None:
                raise ValueError(
                    f"Collocation setting (coll_divs={_divs}, coll_choices={_choices}) "
                    f"does not produce a convergent VIE-1 solver and is not supported "
                    f"(amplification |rho_m| = {rho:.4g} > 1).")
            raise ValueError(
                f"Collocation nodes {np.array2string(node_pos, precision=6)} do not "
                f"produce a convergent VIE-1 solver: the amplification factor "
                f"|rho_m| = prod (1 - c_i)/c_i = {rho:.4g} exceeds 1. Choose nodes "
                f"with a smaller amplification (e.g. include c = 1, as Radau IIA does).")

    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 force_continuous and soln_init_value is None:
        raise ValueError("must specify soln_init_value when force_continuous=True")
    if not force_continuous and soln_init_value is not None and show_warnings:
        print("warning: setting soln_init_value has no effect when force_continuous=False.")

    # Complex dispatch. soln_init_value only contributes to the complex check
    # when force_continuous=True (otherwise it has no effect on the result).
    init_for_complex = soln_init_value if force_continuous else None
    if _samples_indicate_complex([kernel, g], mesh_breakpoints, init_for_complex):
        sample_u = float(np.diff(mesh_breakpoints)[0]) * 0.5
        d_orig = _detect_complex_d_orig(kernel, sample_u)
        init_wrapped = (_block_expand_init(soln_init_value, d_orig)
                        if force_continuous else None)
        result = function_solve_VIE_1(
            kernel=_block_wrap_kernel(kernel, d_orig),
            g=_block_wrap_g(g, d_orig),
            soln_init_value=init_wrapped,
            mesh_breakpoints=mesh_breakpoints,
            coll_nodes=node_pos,
            kernel_singularity=kernel_singularity,
            return_function=return_function,
            force_continuous=force_continuous,
            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}).")

    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)

    if not is_vector:
        g_arr = np.zeros((M, p), dtype=np.float64)
        if g is not None:
            for n in range(M):
                t_n = mesh_breakpoints[n]
                h_n = widths[n]
                for i in range(p):
                    g_arr[n, i] = float(g(t_n + node_pos[i] * h_n))
        if force_continuous:
            # Brunner S_m^(0): degree-m polynomial on {0} ∪ node_pos, with the
            # boundary value carried forward for continuity. Extended weight
            # tensor against the augmented basis (value cols + boundary col).
            cont_basis = _vie1_cont_basis_coefs(node_pos)
            W = _build_W_with_basis_scalar(
                kernel, mesh_breakpoints, node_pos,
                kernel_singularity, _smooth_gl_order, cont_basis)
            adv_U, adv_0 = _vie1_cont_advance(node_pos)
            y, boundary = _dlang_module.function_solve_vie1_cont_d(
                W, g_arr, adv_U, adv_0, float(soln_init_value))
            if return_function:
                polys = _build_vie1_cont_polynomials_scalar(
                    y, boundary, mesh_breakpoints, node_pos)
                return y, _SolutionFunction(polys, mesh_breakpoints, d=0)
            return y

        W = _build_W_scalar(kernel, mesh_breakpoints, node_pos,
                            kernel_singularity, _smooth_gl_order)
        y = _dlang_module.function_solve_vie1_d(W, g_arr)
        if return_function:
            polys = _build_polynomials(y, mesh_breakpoints, node_pos)
            y_func = _SolutionFunction(polys, mesh_breakpoints, d=0)
            return y, y_func
        return y

    # ----- 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 depends only on the kernel; built once and shared across columns. The
    # continuous (S_m^(0)) method uses the augmented value+boundary basis.
    if force_continuous:
        cont_basis = _vie1_cont_basis_coefs(node_pos)
        W = _build_W_with_basis_vector(kernel, mesh_breakpoints, node_pos,
                                       kernel_singularity, _smooth_gl_order, d,
                                       cont_basis)
        adv_U, adv_0 = _vie1_cont_advance(node_pos)
    else:
        W = _build_W_vector(kernel, mesh_breakpoints, node_pos,
                            kernel_singularity, _smooth_gl_order, d)

    sample_t = float(mesh_breakpoints[0] + node_pos[0] * widths[0])
    m = _detect_g_matrix_cols(g, d, sample_t)

    if m is not None:
        # Matrix-valued problem: m simultaneous right-hand sides sharing W.
        if force_continuous:
            init_mat = np.asarray(soln_init_value, dtype=np.float64)
            if init_mat.shape != (d, m):
                raise ValueError(
                    f"soln_init_value must have shape ({d}, {m}) for a "
                    f"matrix-valued problem with force_continuous=True; "
                    f"got shape {tuple(init_mat.shape)}")

        G = _sample_g_at_coll_matrix(g, mesh_breakpoints, node_pos, widths,
                                     M, p, d, m)

        if force_continuous:
            def _col_solve(j):
                g_arr_j = np.ascontiguousarray(G[:, :, :, j])
                return _dlang_module.function_solve_vie1_cont_vec_d(
                    W, g_arr_j, adv_U, adv_0,
                    np.ascontiguousarray(init_mat[:, j]))
            with ThreadPoolExecutor(max_workers=m) as ex:
                results = list(ex.map(_col_solve, range(m)))
            y = np.stack([r[0] for r in results], axis=3)        # (M, p, d, m)
            if return_function:
                boundary = np.stack([r[1] for r in results], axis=2)  # (M+1, d, m)
                polys = _build_vie1_cont_polynomials_matrix(
                    y, boundary, mesh_breakpoints, node_pos, d, m)
                return y, _SolutionFunction(polys, mesh_breakpoints, d=d, m=m)
            return y

        def _col_solve(j):
            g_arr_j = np.ascontiguousarray(G[:, :, :, j])
            return _dlang_module.function_solve_vie1_vec_d(W, g_arr_j)

        with ThreadPoolExecutor(max_workers=m) as ex:
            cols = list(ex.map(_col_solve, range(m)))
        y = np.stack(cols, axis=3)  # (M, p, d, m)

        if return_function:
            polys = _build_polynomials_matrix(y, mesh_breakpoints, node_pos,
                                              d, m)
            return y, _SolutionFunction(polys, mesh_breakpoints, d=d, m=m)
        return y

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

    if force_continuous:
        init_vec = np.asarray(soln_init_value, dtype=np.float64)
        if init_vec.shape != (d,):
            raise ValueError(
                f"soln_init_value must have shape ({d},) for vector kernel; "
                f"got shape {tuple(init_vec.shape)}")
        y, boundary = _dlang_module.function_solve_vie1_cont_vec_d(
            W, g_arr, adv_U, adv_0, init_vec)
        if return_function:
            polys = _build_vie1_cont_polynomials_vector(
                y, boundary, mesh_breakpoints, node_pos, d)
            return y, _SolutionFunction(polys, mesh_breakpoints, d=d)
        return y

    y = _dlang_module.function_solve_vie1_vec_d(W, g_arr)

    if return_function:
        polys = _build_polynomials_vector(y, mesh_breakpoints, node_pos, d)
        y_func = _SolutionFunction(polys, mesh_breakpoints, d=d)
        return y, y_func
    return y