A Residual-Certified Product-Integration Method for Nonlinear Caputo–Volterra Equations

Technical Summary: A Residual-Certified Product-Integration Method for Nonlinear Caputo–Volterra Equations

Problem Statement
The paper addresses the numerical solution of nonlinear Caputo initial-value problems, which can be reformulated as weakly singular Volterra integral equations with Abel-type kernels. Specifically, the study focuses on the equation:
$$C D^\alpha_{0+} u(t) = F(t, u(t)), \quad u(0) = u_0, \quad 0 < \alpha < 1$$
where $F$ satisfies a local Lipschitz condition. The primary challenge addressed is solving the resulting nonlinear algebraic system arising from discretization, particularly when the global Lipschitz constant is large (exceeding the threshold required for standard Banach fixed-point theorems) and when the solution exhibits initial singularities. The paper seeks to determine if a higher-order product-integration rule can be implemented causally while providing a computable, rigorous bound for inexact nonlinear solves.

Methodology
The authors propose an implicit product-integration scheme using piecewise-linear interpolation of the nonlinear term $F(t, u(t))$ on a nonuniform mesh.

1. Discretization: The Abel kernel is integrated exactly against the linear interpolant. This yields a lower-triangular nonlinear system where the value at the current node $U_i$ appears implicitly.
2. Solvability Condition: Unlike strictly causal left-endpoint rules, the current-node term introduces a diagonal coefficient $\kappa_i L_i$. The method relies on a mesh-local solvability test: if $\mu_h = \max_i (\kappa_i L_i) < 1$, the discrete system has a unique solution. This condition is weaker than the standard global condition $S = L^* T^\alpha / \Gamma(\alpha+1) < 1$.
3. Iterative Solvers: Both global Picard iterations and nodewise (sequential) Picard iterations are proven to converge under the local condition $\mu_h < 1$, even when the global constant $S > 1$.
4. Residual Certification: A key component is the derivation of a componentwise algebraic error bound. For any inexact discrete vector $W$, the error $|W - U_h|$ is bounded by $(I - C_h)^{-1} |W - T_h W|$, where $C_h$ is a nonnegative lower-triangular matrix derived from the Lipschitz constants and quadrature weights. This bound is computed via a single forward substitution, separating the error due to incomplete nonlinear iteration from the discretization error.
5. Convergence Analysis:
   • For smooth solutions, the method achieves second-order convergence on sufficiently fine meshes.
   • For solutions with power-type initial singularities ($|g''(t)| \le C t^{\sigma-2}$), the error on a graded mesh ($t_j = T(j/M)^r$) is of order $O(M^{-\min\{r\sigma, 2\}})$.

Key Contributions
The paper claims novelty not in the quadrature rule itself (which is classical) or the use of graded meshes, but in the combination of the following elements:

• Exact Diagonal Test: An explicit, mesh-local condition for the existence and uniqueness of the discrete solution that does not require the global Lipschitz constant to be small.
• Weighted Contraction Argument: A proof that global and nodewise Picard iterations converge under the local condition $\mu_h < 1$.
• Single-Grid Residual Certificate: A computable, componentwise bound for the algebraic error of any inexact iterate. Unlike extrapolation-based estimators, this does not require a refined mesh or asymptotic error expansions.
• Error Separation: A two-level estimate that explicitly separates the error introduced by terminating the nonlinear iteration from the discretization error.
• Mesh-Independent Resolvent Bound: A blockwise Volterra estimate proving that the resolvent bound remains independent of the mesh size, closing the convergence proof.

Results
Numerical experiments validate the theoretical findings across four models:

1. Certificate Reliability: The certified bound consistently exceeds the observed algebraic error for global Picard iterates, even in cases where the global Lipschitz constant $S$ is significantly larger than 1 (e.g., $S \approx 7.7$). The "effectivity" of the estimator (ratio of bound to error) improves as iterations proceed.
2. Convergence Rates:
   • On smooth benchmarks, the proposed method achieves second-order convergence, outperforming the fractional Adams PECE method (order $\approx 1+\alpha$) and the explicit left-endpoint rule (order 1).
   • On problems with initial singularities, uniform meshes yield reduced convergence orders (e.g., $\approx 1.2$), while graded meshes successfully recover the second-order rate.
3. Non-Manufactured Problems: The method performs robustly on a nonlinear problem without a known exact solution, with iteration counts decreasing as the mesh is refined.
4. Computational Cost: The implicit method is computationally more expensive than explicit baselines (Adams PECE, left-endpoint) due to the nonlinear solves and weight construction. However, it provides a rigorous stopping criterion for the nonlinear solver and superior accuracy compared to the left-endpoint rule.

Significance and Claims
The paper modestly positions its contribution as a theoretical and algorithmic refinement of existing product-integration techniques rather than a discovery of new quadrature rules. Its significance lies in:

• Extending the applicability of implicit product-integration to cases where global contraction fails ($S > 1$) by utilizing local mesh properties.
• Providing a rigorous, computable stopping criterion for nonlinear solvers that does not rely on asymptotic assumptions or additional computations on refined meshes.
• Offering a unified framework that separates algebraic and discretization errors, allowing for reliable error estimation in nonlinear fractional Volterra equations.

The authors explicitly state they do not claim the first use of graded meshes, smoothing transformations, or product-integration rules for these equations, but rather the specific integration of a diagonal solvability test, weighted contraction proofs, and a single-grid residual certificate.